clc
clear 

//input data
P0=38*10^5 //combustion chamber pressure in Pa
T0=3500 //combustion chamber temperature in K
ma=41.67 //oxidizer flow rate in kg/s
MR=5 //mixture ratio
k=1.3 //adiabatic constant
R=287 //gas constant in J/kg-K
Pamb=0.0582*10^5 //ambient pressure in Pa
Pe=Pamb //exhaust pressure at sea level in Pa

//calculation 
mf=ma/MR //mass flow of fuel in kg/s 
mp=mf+ma //propellant mass flow in kg/s
Cp=(k*R)/(k-1) //specific heat at constant pressure in J/kg-k
p=P0/Pe //ratio of pressures at combustion chamber and exhaust 
Me=((((p^((k-1)/k))-1)*2)/(k-1))^0.5 //Mach number
t=1/(1+(((k-1)/2)*Me^2)) //ratio of exhaust temperature to combustion temperature
Te=t*T0 //exhaust temperature in Kelvin
a=(1/Me)*(((2/(k+1))+(((k-1)/(k+1))*Me^2))^((k+1)/(2*(k-1)))) //ratio of areas at exit section and throat section of the nozzle
Ce=(k*R*Te)^0.5*Me //exit velocity in the exhaust in m/s
Cj=Ce //average effective jet velocity in m/s, since Pe=Pamb
P1=P0*(2/(k+1))^(k/(k-1)) //pressure at throat section in Pa
T1=T0*(2/(k+1)) //temperature at throat section in K
d1=P1/(R*T1) //density of fuel at throat section in kg/m^3
C1=(k*R*T1)^0.5 //velocity at throat section in m/s
A1=(mp/(d1*C1))*10^4 //nozzle throat area in cm^2
Ae=a*A1 //exit area in cm^2
F=(mp*Ce)*10^-3 //thrust in kN
Cmax1=(2*Cp*T0)^0.5 //maximum possible velocity in m/s
Cf=(F*10^3)/(P0*A1*10^-4) //thrust coefficient, F in kN and A1 in m^2
Cch1=Cj/Cf //characteristic velocity in m/s

//output
printf('(A)nozzle throat area is %3.2f cm^2 \n (B)thrust is %3.1f kN \n (C)thrust coefficient is %3.2f \n (D)characteristic velocity is %3i m/s \n (E)exit velocity in exhaust is %3i m/s\n (F)maximum possible exhaust velocity is %3i m/s\n',A1,F,Cf,Cch1,Ce,Cmax1)