clear; clc; close all;

%% Given
gamma = 1.4;
Rgas = 287;
T0 = 1000;
p0 = 1500e3;
pa = 101e3;
cp = gamma*Rgas/(gamma-1);
h0 = cp*T0;

%% Part A
pr = p0/pa;
Me = sqrt( (2/(gamma-1))*( pr^((gamma-1)/gamma) - 1 ) );
Te = T0/(1+(gamma-1)/2*Me^2);
ae = sqrt(gamma*Rgas*Te);
Ve = Me*ae;
fprintf('Me = %.5f,\nVe = %.3f m/s,\nKE/h0 = %.3f %%\n', Me, Ve, 0.5*Ve^2/h0*100);

%% Prandtl-Meyer
nu = @(M) sqrt((gamma+1)/(gamma-1)).*atan(sqrt((gamma-1)/(gamma+1).*(M.^2-1))) - atan(sqrt(M.^2-1));
nu_inv = @(nu_t) PMInverse(nu_t,gamma);
mach_angle = @(M) asin(1./M);

p_of_M = @(M) p0*(1+(gamma-1)/2*M.^2).^(-gamma/(gamma-1));
T_of_M = @(M) T0*(1+(gamma-1)/2*M.^2).^(-1);

%% Theta_w,max and n
nu_max = nu(Me);
theta_max = nu_max/2;
n  = 50;
rt = 0.02;

N_expected = n*(n+3)/2;
fprintf('theta_w,max = %.4f deg,\nn = %d,\ntotal points N = n(n+3)/2 = %d\n', rad2deg(theta_max), n, N_expected);

%% Increments
theta1 = theta_max/n;
dtheta = (theta_max - theta1)/(n-1);
theta  = theta1 + (0:n-1)'*dtheta;
nu_arr = theta;
Kminus = 2*theta;

corner = struct('x',0,'y',rt);
%% First Char. Line
M1  = nu_inv(nu_arr(1)); mu1 = mach_angle(M1);

s_minus1 = tan(theta(1)-mu1);
x1  = corner.x - corner.y/s_minus1;
Aaxis(1) = struct('x',x1,'y',0,'theta',theta(1),'nu',nu_arr(1),'M',M1,'mu',mu1);
Kplus_axis = zeros(n,1);
Kplus_axis(1) = 0;

% Interior Points along row 1
P = cell(n,n);
for j = 2:n
    a.x = corner.x; a.y = corner.y;
    a.theta = theta(j); a.mu = mach_angle(nu_inv(nu_arr(j)));
    if j == 2
        b = Aaxis(1);
    else
        b = P{1,j-1};
    end
    P{1,j} = unit_process(a,b,Kminus(j),Kplus_axis(1),nu_inv,mach_angle);
end

%% Axis Reflection
for i = 2:n
    th_i = theta(i); mu_i = mach_angle(nu_inv(nu_arr(i)));
    nu_c = Kminus(i); M_c = nu_inv(nu_c); mu_c = mach_angle(M_c);
    s_minus = tan(0.5*((th_i-mu_i)+(0-mu_c)));
    xc = corner.x - corner.y/s_minus;
    Aaxis(i) = struct('x',xc,'y',0,'theta',0,'nu',nu_c,'M',M_c,'mu',mu_c);
    Kplus_axis(i) = -Kminus(i);
    for j = i+1:n
        a2.x = corner.x; a2.y = corner.y;
        if i == 1
            a2.theta = theta(j); a2.mu = mach_angle(nu_inv(nu_arr(j)));
        else
            pa2 = P{i-1,j}; a2 = struct('x',pa2.x,'y',pa2.y,'theta',pa2.theta,'mu',pa2.mu);
        end
        if j-1 == i
            b2 = Aaxis(i);
        else
            b2 = P{i,j-1};
        end
        P{i,j} = unit_process(a2,b2,Kminus(j),Kplus_axis(i),nu_inv,mach_angle);
    end
end

%% Wall Points
W(1)=struct('x',0,'y',rt,'theta',theta_max,'nu',0,'M',1,'mu',pi/2);

for i = 1:n
    if i < n
        ref = P{i,n};
    else
        ref = Aaxis(n);
    end
    pw = W(i);
    theta_w = ref.theta;
    mu_w    = ref.mu;
    s_wall = tan(0.5*(pw.theta + theta_w));
    s_plus = tan(0.5*((ref.theta + ref.mu) + (theta_w + mu_w)));
    xw = (pw.y - ref.y + s_plus*ref.x - s_wall*pw.x)/(s_plus - s_wall);
    yw = pw.y + s_wall*(xw - pw.x);
    W(i+1)=struct('x',xw,'y',yw,'theta',ref.theta,'nu',ref.nu,'M',ref.M,'mu',ref.mu);
end

Lnoz  = W(end).x;  Rexit = W(end).y;  AR = Rexit/rt;
fprintf('Nozzle length L = %.6f m,\nExit half-height = %.6f m,\nAe/At = %.5f\n', Lnoz, Rexit, AR);
fprintf('Exit wall Mach (MOC) = %.5f   (target Me = %.5f)\n', W(end).M, Me);

%% Sample Calculation Table
% Point 1
theta1 = dtheta;  nu1 = theta1;
M1  = nu_inv(nu1);  mu1 = mach_angle(M1);
Km1 = theta1+nu1;   Kp1 = theta1-nu1;
s_minus1 = tan(theta1-mu1);
x1 = 0 - rt/s_minus1;   y1 = 0;

fprintf('\nPoint 1: K-=%.4f  K+=%.4f  th=%.4f  nu=%.4f  M=%.4f  mu=%.4f  x=%.5f  y=%.5f\n', rad2deg(Km1), rad2deg(Kp1), rad2deg(theta1), rad2deg(nu1), M1, rad2deg(mu1), x1, y1);

% Point 2
theta2_corner = 2*dtheta;
mu2_corner    = mach_angle(nu_inv(theta2_corner));
Km2 = 2*theta2_corner;   Kp2 = Kp1;
theta2 = 0.5*(Km2+Kp2);  nu2 = 0.5*(Km2-Kp2);
M2 = nu_inv(nu2);        mu2 = mach_angle(M2);

s_minus2 = tan(0.5*((theta2_corner-mu2_corner)+(theta2-mu2)));
s_plus2  = tan(0.5*((theta1+mu1)+(theta2+mu2)));
x2 = (y1-rt+s_minus2*0-s_plus2*x1)/(s_minus2-s_plus2);
y2 = rt + s_minus2*(x2-0);

fprintf('Point 2: K-=%.4f  K+=%.4f  th=%.4f  nu=%.4f  M=%.4f  mu=%.4f  x=%.5f  y=%.5f\n', rad2deg(Km2), rad2deg(Kp2), rad2deg(theta2), rad2deg(nu2), M2, rad2deg(mu2), x2, y2);

% Point 3
theta3_corner = 3*dtheta;
mu3_corner    = mach_angle(nu_inv(theta3_corner));
Km3 = 2*theta3_corner;   Kp3 = Kp1;
theta3 = 0.5*(Km3+Kp3);  nu3 = 0.5*(Km3-Kp3);
M3 = nu_inv(nu3);        mu3 = mach_angle(M3);

s_minus3 = tan(0.5*((theta3_corner-mu3_corner)+(theta3-mu3)));
s_plus3  = tan(0.5*((theta2+mu2)+(theta3+mu3)));
x3 = (y2-rt+s_minus3*0-s_plus3*x2)/(s_minus3-s_plus3);
y3 = rt + s_minus3*(x3-0);

fprintf('Point 3: K-=%.4f  K+=%.4f  th=%.4f  nu=%.4f  M=%.4f  mu=%.4f  x=%.5f  y=%.5f\n', rad2deg(Km3), rad2deg(Kp3), rad2deg(theta3), rad2deg(nu3), M3, rad2deg(mu3), x3, y3);

%% Coordinate Table
wall_x=[W.x]'; wall_y=[W.y]'; wall_M=[W.M]'; wall_th=[W.theta]';
wall_p=p_of_M(wall_M); wall_T=T_of_M(wall_M);
ContourTbl = table((0:n)', wall_x, wall_y, rad2deg(wall_th), wall_M, wall_p, wall_T, 'VariableNames', {'i','x_m','y_m','theta_deg','Mach','p_Pa','T_K'});
disp(ContourTbl(1:51,:));
writetable(ContourTbl, "coordinatetable.csv")
%% Plots (Contour, Area ratio, Axis & wall properties)
figure; hold on; box on;
for i = 1:n
    xr = Aaxis(i).x;  yr = Aaxis(i).y;
    for j = i+1:n
        xr(end+1) = P{i,j}.x;  yr(end+1) = P{i,j}.y;
    end
    xr(end+1) = W(i+1).x;  yr(end+1) = W(i+1).y;
    plot(xr*1000, yr*1000, '-', 'Color',[0.75 0.75 0.75], 'LineWidth',0.5);
end
for j = 1:n
    xc = corner.x;  yc = corner.y;
    for i = 1:j-1
        xc(end+1) = P{i,j}.x;  yc(end+1) = P{i,j}.y;
    end
    xc(end+1) = Aaxis(j).x;  yc(end+1) = Aaxis(j).y;
    plot(xc*1000, yc*1000, '-', 'Color',[0.75 0.75 0.75], 'LineWidth',0.5);
end
wall_x=[W.x]'; wall_y=[W.y]';
plot(wall_x*1000,  wall_y*1000, 'b-','LineWidth',2);
plot(wall_x*1000, -wall_y*1000, 'b-','LineWidth',2);
yline(0,'k--');

xlabel('x (mm)'); ylabel('y (mm)'); axis equal; grid on;
title(sprintf('Nozzle Contour with MOC Mesh (n = %d)', n));
saveas(gcf,'nozzlecontour.png');

figure; plot(wall_x/Lnoz, wall_y/rt,'r-o','MarkerSize',3); grid on;
xlabel('x / L'); ylabel('A / A_{throat}'); title('Area Ratio vs. Normalized Axial Distance');
saveas(gcf,'arearatio.png')

axis_x = [0; [Aaxis.x]']; axis_M = [1; [Aaxis.M]'];
[axis_x, idx] = sort(axis_x); axis_M = axis_M(idx);
axis_p = p_of_M(axis_M); axis_T = T_of_M(axis_M);

figure;
subplot(3,1,1); plot(axis_x*1000,axis_M,'b-','LineWidth',1.5); ylabel('M'); grid on; title('Centerline Properties');
subplot(3,1,2); plot(axis_x*1000,axis_p/1000,'g-','LineWidth',1.5); ylabel('p (kPa)'); grid on;
subplot(3,1,3); plot(axis_x*1000,axis_T,'r-','LineWidth',1.5); ylabel('T (K)'); xlabel('x (mm)'); grid on;
saveas(gcf,'centerlineproperties.png')

figure;
subplot(3,1,1); plot(wall_x*1000,wall_M,'b-','LineWidth',1.5); ylabel('M'); grid on; title('Wall Properties');
subplot(3,1,2); plot(wall_x*1000,wall_p/1000,'g-','LineWidth',1.5); ylabel('p (kPa)'); grid on;
subplot(3,1,3); plot(wall_x*1000,wall_T,'r-','LineWidth',1.5); ylabel('T (K)'); xlabel('x (mm)'); grid on;
saveas(gcf, "wallproperties.png")

%% MOC unit process
function pt = unit_process(a,b,Km,Kp,nu_inv,mach_angle)
    theta_c = 0.5*(Km+Kp);
    nu_c    = 0.5*(Km-Kp);
    M_c  = nu_inv(nu_c); mu_c = mach_angle(M_c);
    s_minus = tan(0.5*((a.theta-a.mu)+(theta_c-mu_c)));
    s_plus  = tan(0.5*((b.theta+b.mu)+(theta_c+mu_c)));
    xc = (b.y-a.y+s_minus*a.x-s_plus*b.x)/(s_minus-s_plus);
    yc = a.y + s_minus*(xc-a.x);
    pt = struct('x',xc,'y',yc,'theta',theta_c,'nu',nu_c,'M',M_c,'mu',mu_c);
end

%% Newton Raphson Inverse

function M = PMInverse(nu_target,gamma)
tol = 1e-12;
M = 2.0;
A = sqrt((gamma+1)/(gamma-1));
B = sqrt((gamma-1)/(gamma+1));

for k = 1:30
    f = A*atan(B*sqrt(M^2-1)) - atan(sqrt(M^2-1)) - nu_target;
    df = sqrt(M^2-1)/(M*(1+(gamma-1)/2*M^2));

    Mnew = M - f/df;
    if abs(Mnew-M) < tol
        M = Mnew;
        return
    end
    M = Mnew;
end
end