clear all   %Workspace leeren
clc         %Command Window leeren

R=287; %Gaskonstante
cp=1004.5; %Wärmekapazität
dt=10; %Schrittweite Temperatur festlegen
dp=1; %Schrittweite Druck festlegen
p0=0.3; %Startwert Druck
pmax=30; %Maximalwert Druck
T0=233; %Startwert Temperatur  
Tmax=1600; %Maximalwert Temperatur
k=cp/(cp-R); %Isentropenexponent
m0=750; %Massenstrom
my=5.3; % Nebenstromverhältnis
Fn=88000; %Schub
c0=220; %Fluggeschwindigkeit

T0=T0+c0^2/(2*cp);   %um negativen Schub am Nebenstrom zu vermeiden, Isobare nach oben verschieben

p=p0:dp:pmax; %Vektor Druck definieren
T=T0:dt:Tmax; %Vektor Temperatur definieren

[X,Y]=meshgrid(p,T);

L=(Y>=T0.*(X/p0).^(R/cp));  %Werte über Logik filtern

X_filtered=X(L);            %gefilterter Druckwert
Y_filtered=Y(L);            %gefilterter Temperaturwert

A=[X_filtered Y_filtered];  %Matrix mit Druck und Temperatur




%Kernstrom berechnen


p9=zeros(size(X_filtered));  
p9(:)=p0;                          %Annahme unkritisch durchströmt p9=p0
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;                          %Düse kritisch durchströmt falls Bedingung erfüllt
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);%Definition Querschnitt Düsenaustritt
A9(w9==0)=0;                 %Werte für unterste Isobare gleich Null setzen um nicht durch 0 zu teilen




%Nebenstrom berechnen

Ers=(Fn-((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
Ma19(100)=0.5;
Ma19(1000)=0.5;

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


p19=size(Ma19);   %Vektor der Größe Ma19 erzeugen
T19=size(Ma19);   %Vektor der Größe Ma19 erzeugen
p19=p0;           %Vektor p19 erzeugen, der im Falle Ma19>1 (Düse kritisch) überschrieben wird

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


%falls Düse kritisch : p19 neu berechnen

x0=1;   %Startwert für 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)m0*(my./(1+my))*(sqrt(k*R*T0*(p19/p0).^(R/cp))-c0)+(p19-p0)*(R*T0*(p19./p0).^(R/cp))./(p19*sqrt(k*R*T0*(p19./p0).^(R/cp)))-Ers(i);  %Bestimmungsgleichung für p19
        
        [p19(i), fval, exitflag]=fsolve(F,x0,options);  %lösen der Gleichung
           
        
        if exitflag <= 0
            p19(i) = NaN; % kein Ergebnis
        else
            x0 = p19(i);  %Startwert gleich Ergebnis von fslove
        end
        
        
    end
end
        

%for i=1:size(Ers)
    % F=@(p19)m0*(my./(1+my))*(sqrt(k*R*T0*(p19/p0).^(R/cp))-c0)+(p19-p0)*(R*T0*(p19./p0).^(R/cp))./(p19*sqrt(k*R*T0*(p19./p0).^(R/cp)))-Ers(i);
    %  x0=1;     %Startwert Iteration
    %[p19(i) fval exitflag] = fsolve(F,x0,options);
    %p19(i)=fsolve(F,x0,options);   %Gleichung nach p19 lösen
    % x0=p19(i);  %Neuer Startwert für Iteration ist Lösung der vorherigen Gleichung
%end
%fsolve('Ers(1)=m0*(my/(1+my))*(sqrt(k*R*T0*(p19/p0)^(R/cp))-c0)+(p19-p0)*(R*T0*(p19/p0)^(R/cp))/(p19*sqrt(k*R*T0*(p19/p0)^(R/cp)))','p19');

T19(Ma19>1)=T0.*(p19(Ma19>1)./p0).^(R/cp);