psi=7;%循環靜溫比
f=0.0675;%燃料空氣當量比
fh_pr=2818000;%f*hpr[J/kg]
h_pr=fh_pr/f;%反應熱[J/kg]
T0=222;%參考溫度K
R=289.3;%熱力學常數
C_pc=1090;%壓縮過程定壓比熱容[J/（kg*K）]
C_pb=1510;%熱附加過程定壓比熱容[J/（kg*K）]
gamma_c=1.362;%壓縮過程比熱比
gamma_b=1.238;%熱附加過程比熱比
a0=sqrt（1.4*R*T0）;%入口聲速[m/s]
F=zeros（413）;%未裝機比推力[N*s/kg]
eta_0=zeros（413）;%總效率
eta_c=zeros（413）;%絕熱壓縮效率
eta_e=zeros（413）;%絕熱膨脹效率
I_sp=zeros（413）;%比沖[s]
for?i=1:41
????for?j=1:3
????Ma0=4+0.2*（i-1）;????
????V0=Ma0*a0;
????Sa0=V0*（1+R*T0/V0^2）;
????T3=psi*T0;
????V3=sqrt（V0^2-2*C_pc*T0*（psi-1））;
????%Sa3=V3*（1+R*T3/V3^2）;
????c=0.005+j*0.005;%取值0.01、0.015、0.02
????psi1=1+（（gamma_c-1）/2）*Ma0^2*（1-（1+c）^2*（1+（2-2*psi）/（（gamma_c-1）*Ma0^2）））;
????eta_c_=（1-0.05*j）*（（1-1/psi1）/（1-1/psi））;%eta1取值0.95、0.9、0.85
????eta_c（ij）=eta_c_;
????ratio_P_3_0=（psi/（psi*（1-eta_c_）+eta_c_））^（C_pc/R）;%壓縮部件壓強比
????ratio_A_3_0=psi*（V3/V0）/ratio_P_3_0;%壓縮部件面積比
????b=0.05+0.05*j;%取值0.1、0.15、0.2
????V4=V3*（（1+f*0.5）/（1+f）-b/（1+f））;