R=287; %Gaskonstante
cp=1004.5; %Wärmekapazität
dt=1; %Schrittweite Temperatur festlegen
dp=0.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

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



% AB HIER SCHEITERT ES....

%idx2=Ma19>1;         %Düse kritisch durchströmt falls Bedingung erfüllt


% solve('Ers(idx)=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')