clear all
clc
Aneck = 0.15;
Aexit = 0.15+2*5*sind(13);
gamma = 1.4;
p0 = 7000;
syms M
pNeck = p0/(1+(gamma-1)/2)^(gamma/(gamma-1));

Mexit = vpasolve(Aexit/Aneck == 1/M*(2/(gamma+1)*(1+(gamma-1)/2*M^2))^((gamma+1)/(2*(gamma-1))), M)

pe = p0/(1+(gamma-1)/2*Mexit(2)^2)^(gamma/(gamma-1))

% h = 11486m

grad = (-0.075-(-2))/2;
n = 50;
x = 0:(2+5*cosd(13))/n:  2+ 5*cosd(13);
for i = 1:length(x)
    if x(i) <= 2
        A(i) = 4-2*grad*x(i);
        MTemp = vpasolve(A(i)/0.15 == 1/M*(2/(gamma+1)*(1+(gamma-1)/2*M^2))^((gamma+1)/(2*(gamma-1))), M);
        Mn(i) = MTemp(1);
    else
        A(i) = 0.15+2*tand(13)*(x(i)-2);
        MTemp = vpasolve(A(i)/0.15 == 1/M*(2/(gamma+1)*(1+(gamma-1)/2*M^2))^((gamma+1)/(2*(gamma-1))), M);
        Mn(i) = MTemp(2);
    end
    p(i) = p0/(1+(gamma-1)/2*Mn(i)^2)^(gamma/(gamma-1));
end
plot(x,Mn,x,p/p0)

rhoExit = rho
thrust = Mn(length(x))