close all
clear
clc

graphics_toolkit('gnuplot');

% Gas parameters
% -----------
% Nitrogen --
% ------------
Rgas    = 296.80;
gamma   = 7/5;
cpOverR = 7/2;
% -----------

%% --------
%% Argon --
%% --------
%Rgas    = 208.13;
%gamma   = 5/3;
%cpOverR = 5/2;
%% --------

% Load nozzle geometry
XYZgeom = load('expanding_nozzle_geometry.dat');
Xgeom   = XYZgeom(:,1);
Ygeom   = XYZgeom(:,2);

% Mesh
Nx = 20000;
xx = linspace(min(Xgeom), max(Xgeom), Nx)';
yy = interp1(Xgeom, Ygeom, xx);

A_vect = pi*yy.^2; % Nozzle area

dAdx = diff(A_vect)./diff(xx); 
dAdx = [dAdx; dAdx(end)]; 

% Initial conditions (at sonic point)
P0   = 3.0189e7; % [Pa]
rho0 = 105.53;   % [kg/s]

% Estimate u0
c0 = sqrt(gamma*P0/rho0);
u0 = c0 + 1.1;

P_vect   = zeros(size(xx));
u_vect   = zeros(size(xx));
rho_vect = zeros(size(xx));

P_vect(1)   = P0;
u_vect(1)   = u0;
rho_vect(1) = rho0;

% Integrate
for id_x = 2:numel(xx)

  dx = xx(id_x) - xx(id_x-1); % mesh size

  % Previous position variables
  P = P_vect(id_x-1);
  u = u_vect(id_x-1);
  rho = rho_vect(id_x-1);
  A   = A_vect(id_x-1);

  % RHS
  dPdx   = rho*u^2/( A*(1+rho*u^2/P*(1/cpOverR - 1)) ) * dAdx(id_x - 1);
  dudx   = - 1/(rho*u) * dPdx;
  drhodx = u*rho^2/P/cpOverR*(1-cpOverR) * dudx;

  % Integrate
  P_vect(id_x)   = P_vect(id_x-1)   + dx*dPdx;
  u_vect(id_x)   = u_vect(id_x-1)   + dx*dudx; 
  rho_vect(id_x) = rho_vect(id_x-1) + dx*drhodx; 

end

% Post processing 
T_vect = P_vect./(rho_vect.*Rgas);
c_vect = sqrt(gamma*Rgas*T_vect);
M_vect = u_vect./c_vect;

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

figure
plot(xx, u_vect, 'linewidth', 2)
xlabel('x [m]', 'fontsize', 14)
ylabel('Velocity', 'fontsize', 14)


