% Given
Ma=0.44
deltaH=17940;
W_cp=0.97;
P_a=12.93252;
m_dot=26;
T_t4=2150;
T_a=506.21;
gamma=1.4;
R=53.35;
Cp=0.2399565;
rho_8= 7.952636E-04;
% Diffuser
a_a=sqrt(gamma*R*32.17*T_a)
u_a=Ma*a_a
Tta_Ta=(1+0.5*(gamma-1).*Ma^2)
T_ta=506.21 * 1.03872
T_t2=T_ta
Pta_Pa=(Tta_Ta)^((gamma)/(gamma-1))
Pt2_Pa=Pta_Pa
P_t2=Pta_Pa*P_a
pi_c=linspace(5,30);
%Compressor
P_t3=pi_c.*P_t2
tau_c=pi_c.^((gamma-1)./gamma)
T_t3=tau_c*T_t2
%Primary Combustor
P_t4=P_t3
m_dotf=(m_dot*Cp*(T_t4-T_t3))./ deltaH
%Turbine
T_t5=(((m_dotf.*Cp.*(T_t3-T_t2)+W_cp*m_dotf*Cp*T_a)./(m_dotf*Cp))-T_t4)*-1
tau_t=T_t5./T_t4
pi_t= 0.2157084
P_t5=pi_t.*P_t4
%Primary Nozzle
T_t6=T_t5
P_t6=P_t5
P_8= 12.93252 %exit pressure matches ambient pressure
T_8=T_t6.*((P_8./P_t6).^((gamma-1)/gamma))
u_8=[2*Cp.*(T_t6-T_8)*32.17*778.16].^0.5
a_8=(gamma*R*T_8.*32.17).^0.5
M_8=u_8./a_8
A_8=m_dotf.*144./(rho_8.*u_8*32.17)
% Function for the work coefficient
tau_5=(T_ta./T_a).*tau_c.*tau_t;
a=T_t4./T_a
b=T_ta./T_a
c=T_a./T_ta
d=[((a*c.*((tau_t-1)./tau_c))./b-1).^(0.5)-1]
C_we=(a*(1-(tau_5./(b*tau_c)))-(b*(tau_c-1)))+((gamma-1).*(Ma.^2).*((((a*c.*((tau_5-1)./tau_c))./(b-1)).^(0.5))-1))
%%Graph
% 1. Cwe vs pi_c
f1=figure(1);clf;
C_we=(a.*(1-(tau_5./(b*tau_c)))-(b*(tau_c-1)))+((gamma-1)*(Ma.^2).*((((a*c.*((tau_5-1)./tau_c))./(b-1)).^(0.5))-1))
plot (pi_c,C_we,'k')
title('Figure 1. C_we vs. \pi_c');
xlabel('\pi_c')
ylabel('C_w_e')
