format long e
clc
clear all

Cp=1004;
R=283;
gamma=1.4;
steps=100001;


Cf=0;
dq=0;
dshaft=0;
ductlength=1;
efficiency=1;
A(1)=.1;
A_exit=.4345;
Twall=0;
P(1)=101325;
T(1)=288;
rho(1)=P(1)/(R*T(1));
M(1)=0.5;
Pt(1)=P(1)*(1+(gamma-1)/2*M(1)^2)^(gamma/(gamma-1));
Tt(1)=T(1)*(1+(gamma-1)/2*M(1)^2);
    
for i=2:steps %Ctrl+c to stop code from running

    dx=ductlength/steps;

    rho(1)=P(1)/(R*T(1));
    a(1)=sqrt(gamma*R*T(1));
    U(1)=M(1)*a(1);
    

    dA=(A_exit-A(1))/(steps-1);
    A(i)=A(i-1)+dA;

    C(i)=2*pi*sqrt(A(i)/pi);

    U(i)=(efficiency*dshaft-(.5*Cf*C(i)*dx*U(i-1)^2+dA)/A(i-1)-(dq+dshaft)/(Cp*T(i-1)))/(U(i-1)-1/U(i-1)-1/(Cp*T(i-1)))+U(i-1);
    dU=U(i)-U(i-1);

    T(i)=(dq+dshaft-U(i-1)*dU)/Cp+T(i-1);
    dT=T(i)-T(i-1);

    rho(i)=(rho(i-1)*U(i-1)*A(i-1))/(U(i)*A(i));

    P(i)=-.5*Cf*C(i)*dx/A(i-1)+efficiency*dshaft*rho(i-1)-U(i-1)*dU*rho(i-1)+P(i-1);

    M(i)=U(i)/sqrt(gamma*R*T(i));

    Tt(i)=T(i)*(1+(gamma-1)/2)*M(i)^2;

    Pt(i)=P(i)*(1+(gamma-1)/2)^(gamma/(gamma-1));

   
end




disp(['Exit pressure = ' num2str(P(steps)),' Pa'])
disp(['Exit temperature = ' num2str(T(steps)),' K'])
disp(['Exit velocity = ' num2str(U(steps)),' m/s'])
disp(['Exit Mach = ' num2str(M(steps))])
disp(['Exit total pressure = ' num2str(Pt(steps)),' Pa'])
disp(['Exit total temperature = ' num2str(Tt(steps)),' K'])
disp(['Exit density = ' num2str(rho(steps)), ' kg/m^2'])
disp(['P/P*= ' num2str(P(1)/P(i))])
disp(['rho/rho*= ' num2str(rho(1)/rho(i))])
disp(['T/T*= ' num2str(T(1)/T(i))])
disp(['U/U*= ' num2str(U(1)/U(i))])
disp(['Pt/Pt*= ' num2str(Pt(1)/Pt(i))])
disp(['Tt/Tt*= ' num2str(Tt(1)/Tt(i))])