close all
clear
clc

% ============================================================================
% This plots the anaylitical result for the Sod shock problem.
% Equations can be found on Quartapelle & Auteri, Fluidodinamica Comprimibile.
%
% Subscripts are (assuming that PL > PR):
% L) left state
% 1) between rarefaction wave and contact discontinuity
% 2) between contact discontinuity and shock wave
% R) right state
% ============================================================================

% Time for plotting the solution
tPLOT = 0.228; % [s]

% Gas properties and physical constants
% gamma = 5/3;
gamma = 7/5;
m = 2.18e-25;
kB = 1.38e-23; % [J/K]
R = kB/m;

% Left and Right states
rhoL = 1; % [kg/m3]
uL   = 0; % [m/s] IMPOSED
PL   = 1; % [Pa]

rhoR = 1/8;  % [kg/m3]
uR   = 0; % [m/s] IMPOSED
PR   = 1/10; % [Pa]

% STEP 0) compute intermediate pressure, where post-rarefaction velocity is equal to post-shock velocity
cL = sqrt(gamma*PL/rhoL);
cR = sqrt(gamma*PR/rhoR);

Pfun = @(P) 2*cL/(gamma-1)*(1 - (P/PL).^((gamma-1)/(2*gamma))) - sqrt(2)*cR*(P/PR - 1)./sqrt( gamma*((gamma+1)*P/PR + (gamma-1)));

P1_guess = (PL + PR)/2;
P1 = fsolve(Pfun, P1_guess);

% STEP 1) Find conditions after rarefaction wave (conditions "1")
rho1 = rhoL*(P1/PL)^(1/gamma);
u1 = uL + 2*cL/(gamma-1)*(1 - (rho1/rhoL)^((gamma-1)/2));

% STEP 2) Compute results post-shock conditions (conditions "2")
P2 = P1; 
u2 = u1;

% STEP 3) Find the post-shock density from shock relations
MR = sqrt( 1/(2*gamma)*((gamma-1) + (gamma+1)*P2/PR) );
rho2 = rhoR*((gamma+1)*MR^2)/((gamma-1)*MR^2+2);

% +++++++++++++++++++ CREATE RESULT +++++++++++++++++++++++
xL = (uL-cL)*tPLOT;

c1 = sqrt(gamma*P1/rho1);
x_end_rar = (u1 - c1)*tPLOT;

x_dc = u1*tPLOT;

uShock = MR*cR;
xR = uShock*tPLOT;

xMAX = 1.5*xL;
xMIN = 1.5*xR;

xx = linspace(xMIN, xMAX, 1000);

rho = zeros(size(xx));
u   = zeros(size(xx));
P   = zeros(size(xx));

idL = find(xx <= xL);
idR = find(xx >= xR);
id_rar = find((xx > xL ) & (xx <= x_end_rar));
id_1   = find((xx > x_end_rar) & (xx <= x_dc));
id_2   = find((xx > x_dc)      & (xx <= xR));

rho(idL) = rhoL;
rho(idR) = rhoR;
rho(id_rar) = rhoL*(2/(gamma+1) + (gamma-1)/(gamma+1)*(uL-xx(id_rar)/tPLOT)/cL).^(2/(gamma-1));
rho(id_1)   = rho1;
rho(id_2)   = rho2;

u(idL) = uL;
u(idR) = uR;
u(id_rar) = uL*(gamma-1)/(gamma+1) + 2/(gamma+1)*(cL + xx(id_rar)/tPLOT);
u(id_1)   = u1;
u(id_2)   = u2;

P(idL) = PL;
P(idR) = PR;
P(id_rar) = PL*(rho(id_rar)/rhoL).^gamma;
P(id_1)   = P1;
P(id_2)   = P2;

% ---------------- PLOT ----------------
figure
%subplot(1,3,1)
plot(xx, rho, 'b', 'linewidth', 2);
xlim([min(xx), max(xx)])
ylabel('Density [kg/m3]')
xlabel('Position [m]')

%subplot(1,3,2)
figure
plot(xx, u, 'b', 'linewidth', 2);
xlim([min(xx), max(xx)])
ylabel('Velocity [m/s]')

%subplot(1,3,3)
figure
plot(xx, P, 'b', 'linewidth', 2);
xlim([min(xx), max(xx)])
ylabel('Pressure [Pa]')