clear all                                                                  %Workspace leeren
clc                                                                        %Command Window leeren

R=287;                                                                     %Gaskonstante
cp=1004.5;                                                                 %Wärmekapazität
tN=100;                                                                    %Schrittweite Temperatur festlegen
pN=20;                                                                     %Schrittweite Druck festlegen
p0=0.3*10^5;                                                               %Startwert Druck in hPa
pMax=30*10^5;                                                              %Maximalwert Druck in hPa
t0=233;                                                                    %Startwert Temperatur  
tMax=1600;                                                                 %Maximalwert Temperatur
k=cp/(cp-R);                                                               %Isentropenexponent
m0=750;                                                                    %Massenstrom
mY=5.3;                                                                    %Nebenstromverhältnis
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

A=zeros(size(T,2),2);                                                      %leere Matrix A erstellen

for i=size(X,2):-1:1;                                                      %Matrix A mit Werten aus Meshgrid beschreiben
B=[X(:,i) Y(:,i)];
A=vertcat(B,A);
end


A(find(A(:,1) == 0),:) = [];                                               %Nullwerte von A=zeros... herauslöschen




delta_S=cp*log(A(:,2)./t0)-R*log(A(:,1)./p0);                              %delta_S berechnen, Entropiedifferenz zweier Punkte im Rechengebiet                   
A(delta_S<0,:)=[] ;                                                        %für delta_S<0  Werte filtern

t_id=A(:,2).*(p0./A(:,1)).^(R/cp);                                         %T_id berechnen
c_id=sqrt(2*cp.*(A(:,2)-t_id(:)));                                         %c_id berechnen
A(c_id<c0,:)=[];                                                           %falls c_id<c0 soll das Wertepaar (Druck;Temperatur) in Matrix A gelöscht werden

[Zeile Spalte]=size(A);                                                    %für folgende Berechnungen Größe der Matrix A ausgeben





%% Kernstrom berechnen

p9=repmat(p0,Zeile,1);                                                     %Vektor p9 erzeugen,mit Anzahl Zeilen Matrix A und Werten p0
                                                                           %da Annahme: unkritisch durchströmt

                                                                           
t9=A(:,2).*(p9./A(:,1)).^(R/cp);                                           %statische Temperatur Düsenaustritt
w9=sqrt(2*cp.*(A(:,2)-t9));                                                %Geschwindigkeit Düsenaustritt
Ma9=w9./sqrt(k*R.*t9);                                                     %Machzahl berechnen

idx=Ma9>1;                                                                 %Prüfen ob Düse kritisch durchströmt ist und ggf. T9,p9 und w9 neu berechnen
t9(idx) = (2./(k+1)).*A(idx,2); 
p9(idx)=A(idx,1).*(t9(idx)./A(idx,2)).^(cp/R);
w9(idx)=sqrt(2*cp.*(A(idx,2)-t9(idx)));


A9=(m0/(1+mY))*R.*(t9./(p9.*w9));                                          %Berechnung Querschnitt Düse Kernstrom
                                                               






%% Nebenstrom berechnen



Ers=(F_n-((m0*(1/(1+mY)).*(w9-c0)+A9.*(p9-p0)-m0*c0)))/(m0*(mY/(1+mY)));   %zur Übersichtlichkeit Variable "Ers" definieren
p19=p0;                                                                    %Annahme unkritisch durchströmt p19=p0
t19=t0;                                                                    %dann Vorgabe laut Aufgabenstellung: T0=T19
w19=Ers+c0;                                                                %Geschwindigkeit Düsenautritt
      
Ma19=w19./(sqrt(k*R*t19));                                                 %Machzahl Düsenaustritt

tt19=t19+w19.^2./(2*cp);                                                   %Totaltemperatur Nebenstrom
pt19=p19*(tt19./t19).^(cp/R);                                              %Totaldruck Nebenstrom


                                                                           %Prüfen ob Düse kritisch durchströmt

                           

p19=repmat(p19,size(Ma19));                                                %Vektor p19 erzeugen,mit Größen Ma19 und Werten p0
t19=repmat(t19,size(Ma19));                                                %Vektor T19 erzeugen,mit Größen Ma19 und Werten T19



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

                                                                           

x0=51245;                                                                  %Startwert Iteration

for i=1:size(Ers)

    if Ma19(i) > 1                                                         %für alle Ma19 für die gilt Ma19>1, muss p19 neu berechnet werden
        F=@(p19)(sqrt(k*R*t0*(p19/p0).^(R/cp))...                          %Bestimmungsgleichung für p19
        -c0)+(p19-p0)*(R*t0*(p19./p0).^(R/cp))./...
        (p19*sqrt(k*R*t0*(p19./p0).^(R/cp)))-Ers(i);                        
        
        [p19(i), fval, exitflag]=fsolve(F,x0,options);                     %lösen der Bestimmungsgleichung nach p19
            if exitflag <= 0
                p19(i) = NaN;                                              %falls Bedingung erfüllt,gibt es keine Lösung die konvergiert
              else
                x0 = p19(i);                                               %Startwert gleich Ergebnis von fslove
            end
    end
end

 
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./(p19.*w19));                                   %Berechnung Querschnitt Düse 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)+A9.*(p9-p0)-m0*(1/(1+mY))*c0;            %Vortriebsleistung Kernstrom
P_vor_Neben=m0*(mY/(1+mY)).*(w19-c0)+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_Kern=P_vor_Kern./(m0*(1/(1+mY)).*ex_Kern);                         %Exergetischer Wirkungsgrad Kernstrom
ex_Wirk_Neben=P_vor_Neben./(m0*(mY/(1+mY)).*ex_Neben);                     %Exergetischer Wirkungsgrad Nebenstrom

mesh(p,T,ex_Wirk_ges);

%%Excel Dateien schreiben

xlswrite('Exergieberechnung',{'pt9','tt9','p9','t9','w9','Ma9','p19','t19','w19','Ma19'},1,'A1:J1')
xlswrite('Exergieberechnung',{'ex_Wirk_ges','ex_Wirk_Kern','ex_Wirk_Neben'},1,'L1:N1')

xlswrite('Exergieberechnung',A(:,1),1,'A2')
xlswrite('Exergieberechnung',A(:,2),1,'B2')

xlswrite('Exergieberechnung',p9,1,'C2')
xlswrite('Exergieberechnung',t9,1,'D2')
xlswrite('Exergieberechnung',w9,1,'E2')
xlswrite('Exergieberechnung',Ma9,1,'F2')

xlswrite('Exergieberechnung',p19,1,'G2')
xlswrite('Exergieberechnung',t19,1,'H2')
xlswrite('Exergieberechnung',w19,1,'I2')
xlswrite('Exergieberechnung',Ma19,1,'J2')

xlswrite('Exergieberechnung',ex_Wirk_ges,1,'L2')
xlswrite('Exergieberechnung',ex_Wirk_Kern,1,'M2')
xlswrite('Exergieberechnung',ex_Wirk_Neben,1,'N2')