%clear all                                                                 %Workspace leeren
%clc                                                                       %Command Window leeren

R=287;                                                                     %Gaskonstante
cp=1004.5;                                                                 %Waermekapazitaet
tN=100;                                                                    %Anzahl Stuetzstellen Temperatur
pN=20;                                                                     %Anzahl Stuetzstellen Druck
p0=0.3;                                                                    %Startwert Druck in hPa
pMax=30;                                                                   %Maximalwert Druck in hPa
t0=233;                                                                    %Startwert Temperatur  
tMax=1600;                                                                 %Maximalwert Temperatur
k=cp/(cp-R);                                                               %Isentropenexponent
m0=750;                                                                    %Massenstrom
mY=5.3;                                                                    %Nebenstromverhaeltnis
F_n=88000;                                                                 %Schub
c0=220;                                                                    %Fluggeschwindigkeit







%% Rechengebiet erstellen


p=linspace(p0,pMax,pN);                                                    %Vektor Druck definieren
T=linspace(t0,tMax,tN);                                                    %Vektor Temperatur definieren

[X,Y]=meshgrid(p,T);                                                       %Wertekombinationen per Meshgrid erstellen

delta_S=cp*log(Y./t0)-R*log(X./p0);                                        %delta_S berechnen, Entropiedifferenz zweier Punkte im Rechengebiet 
                                                             
t_id=Y.*(p0./X).^(R/cp);                                                   %t_id berechnen

c_id=sqrt(2*cp.*(Y-t_id));                                                 %c_id berechnen

L1=delta_S<0; 
L2=c_id<c0 ;

X(L1|L2)=0;                                                                %falls L1 oder L2 erfuellt soll das entsprechende Element der Matrix X gleich 0 gesetzt werden
Y(L1|L2)=0;                                                                %falls L1 oder L2 erfuellt soll das entsprechende Element der Matrix Y gleich 0 gesetzt werden

%% Kernstrom
                                                                           %Annahme: unkritisch durchströmt
p9=p0*ones(size(X));                                                       %Matrix p9 mit Werten p0 erzeugen,gleiche Größe wie Matrix X                    
t9=Y.*(p9./X).^(R/cp);                                                     %Matrix t9 berechnen 
w9=sqrt(2*cp.*(Y-t9));                                                     %Matrix w9 berechnen 
Ma9=w9./sqrt(k*R.*t9);                                                     %Matrix Ma9 berechnen 


L3=Ma9>1;                                                                  %falls kritisch durchströmt Matrizen t9,p9,w9 neu berechnen
t9(L3) = (2./(k+1)).*Y(L3); 
p9(L3)=X(L3).*(t9(L3)./Y(L3)).^(cp/R);
w9(L3)=sqrt(2*cp.*(Y(L3)-t9(L3)));


A9=(m0/(1+mY))*R.*(t9./(10^5.*p9.*w9));                                    %Querschnitt A9 berechnen


%% Nebenstrom 

Ers=(F_n-((m0*(1/(1+mY)).*(w9-c0)+10^5.*A9.*(p9-p0)-m0*c0)))...
    /(m0*(mY/(1+mY)));                                                     %Matrix 'Ers' mit entsprechenden Konstanten definieren
Ers(isnan(Ers))=0;                                                         %falls in Matrix 'Ers' NaN erscheint: gleich Null setzen

p19=p0*ones(size(X));                                                      %Matrix p19 mit Werten p0 erzeugen,gleiche Groeße wie Matrix X                                                                   
t19=t0*ones(size(X));                                                      %Matrix t19 mit Werten t0 erzeugen,gleiche Groeße wie Matrix X               
w19=Ers+c0;                                                                %Matrix w19 erzeugen
w19(w19==c0)=0;                                                            %falls Element der Matrix w19=c0 ist, gleich Null setzen, da diese Werte aus Rechengebiet fallen
Ma19=w19./(sqrt(k*R*t19));                                                 %Matrix Ma19 berechnen




options = optimset('Display', 'off');                                      %Iterationsanzeige im Command-Window ausschalten

x0=51245;                                                                  %Startwert Iteration

for j=1:size(Ers,2)
    
for i=1:size(Ers,1)

    if Ma19(i,j) > 1                                                        
        F=@(p19)(sqrt(k*R*t0*(p19./p0).^(R/cp))-c0)+10^5.*(p19-p0)*(R*t0*...     %Bestimmungleichung für p19
        (p19./p0).^(R/cp))./(10^5.*p19.*sqrt(k*R*t0*(p19./p0).^(R/cp)))-Ers(i,j);                        
        
        [p19(i,j), fval, exitflag]=fsolve(F,x0,options);                    
            if exitflag <= 0
                p19(i,j) = NaN;                                             
              else
                x0 = p19(i,j);                                             %das berechnete p19 als neuen Startwert nutzen 
            end
    end
end
end
p19(p19==p0)=0;


t19(Ma19>1)=t0.*(p19(Ma19>1)./p0).^(R/cp);                                 %für alle Ma19 für die gilt Ma19>1, muss t19 neu berechnet werden
w19(Ma19>1)=sqrt(k*R.*t19(Ma19>1));                                        %für alle Ma19 für die gilt Ma19>1, muss w19 neu berechnet werden

A19=m0*(mY/(1+mY))*R.*(t19./(10^5.*p19.*w19));                                   %Berechnung Querschnitt Duese Nebenstrom


%% Exergieberechnung Kernstrom

ex_Kern=cp.*(t9-t0)-t0*(cp*log(t9./t0)-R*log(p9./p0))+0.5.*(w9-c0).^2;     %Berechnungsvorschrift Exergie Kernstrom

%% Exergieberechnung Nebenstrom

ex_Neben=cp.*(t19-t0)-t0*(cp*log(t19./t0)-R*log(p19./p0))+0.5*(w19-c0).^2; %Berechnungsvorschrift Exergie Nebenstrom

%% Exergetischer Wirkungsgrad

P_vor=F_n*c0;                                                              %Gesamte Vortriebsleistung
P_vor_Kern=m0*(1/(1+mY)).*(w9-c0)+10^5.*A9.*(p9-p0)-m0*(1/(1+mY))*c0;            %Vortriebsleistung Kernstrom
P_vor_Neben=m0*(mY/(1+mY)).*(w19-c0)+10^5.*A19.*(p19-p0)-m0*(mY/(1+mY)).*c0;     %Vortriebsleistung Nebenstrom

ex_Wirk_ges=P_vor./(m0.*(ex_Kern+ex_Neben));                               %Exergetischer Wirkungsgrad "gesamt"
ex_Wirk_ges(isnan(ex_Wirk_ges))=0;
ex_Wirk_Kern=P_vor_Kern./(m0*(1/(1+mY)).*ex_Kern);                         %Exergetischer Wirkungsgrad "Kern"
ex_Wirk_Kern(isnan(ex_Wirk_Kern))=0;
ex_Wirk_Neben=P_vor_Neben./(m0*(mY/(1+mY)).*ex_Neben);                     %Exergetischer Wirkungsgrad "Neben"
ex_Wirk_Neben(isnan(ex_Wirk_Neben))=0;