close all
clear all

k=1.1323;
pc=25*10^5;
pa=[101325, 54012, 12041, 2510]';
cstar=1550;
Me=(2:0.05:4);
e=((k+1)./2).^(-(k+1)./(2.*(k-1))).*(1+(k-1)./2*Me.^2).^((k+1)./(2.*(k-1)))./Me;
pe=pc.*(1+(k-1)./2*Me.^2).^(k./(1-k));

cF = k.*sqrt((2./(k-1)).*(2./(k+1)).^((k+1)./(k-1)).*(1-(pe/pc).^((k-1)./k)))+(pe-pa)./pc.*e;
Isp=cF.*cstar;

figure(1)
plot(Me,e)
grid
xlabel("$M_e$ [-]",'Interpreter','latex')
ylabel("$\epsilon$ [-]",'Interpreter','latex')

figure(2)
plot(Me,pe)
grid
xlabel("$M_e$ [-]",'Interpreter','latex')
ylabel("$p_e$ [Pa]",'Interpreter','latex')

figure(3)
plot(Me,Isp)
grid
xlabel("$M_e$ [-]",'Interpreter','latex')
ylabel("$I_{sp}$ [m/s]",'Interpreter','latex')
ylim([0 3500])