close all
clear
clc

% Gas properties and throat conditions
gamma = 5/3; % cp/cv

P0   = 3.0189E7; % [Pa] 
rho0 = 105.54;   % [kg/m^3]

% Compute throat velocity
c0 = sqrt(gamma*P0/rho0); % [m/s] speed of sound
u0 = c0;                  % [m/s] sonic

% Load nozzle data
dd = load('expanding_nozzle_geometry.dat');
x_vect = dd(:,1);
y_vect = dd(:,2);

A_vect = pi*y_vect.^2;
A_crit = A_vect(1);

% Create Mach-area relation
M_test = linspace(1, 200, 10000);  % max Mach: 200
A_test = A_crit./M_test.*( 2*(1+(gamma-1)/2*M_test.^2) / (gamma+1) ).^((gamma+1)/(2*(gamma-1)));

% Find Mach number in the loaded nozzle
M_vect = interp1(A_test, M_test, A_vect);

% Find velocity
u_vect = u0 * M_vect.*sqrt((gamma+1)/2./(1+(gamma-1)/2*M_vect.^2));

% Find pressure
P_vect = P0 * ((gamma+1)/2./(1+(gamma-1)/2*M_vect.^2)).^(gamma/(gamma-1));

% Find density
rho_vect = rho0 * (P_vect/P0).^(1/gamma);

% PLOT
figure
plot(x_vect, M_vect, 'b', 'linewidth', 2)
xlabel('x [m]', 'fontsize', 14)
ylabel('Mach number', 'fontsize', 14)

figure
semilogy(x_vect, u_vect, 'k', 'linewidth', 2)
xlabel('x [m]', 'fontsize', 14)
ylabel('velocity', 'fontsize', 14)

figure
semilogy(x_vect, P_vect, 'k', 'linewidth', 2)
xlabel('x [m]', 'fontsize', 14)
ylabel('pressure', 'fontsize', 14)

figure
semilogy(x_vect, rho_vect, 'k', 'linewidth', 2)
xlabel('x [m]', 'fontsize', 14)
ylabel('density', 'fontsize', 14)
