%% Shockwave Calcs
clc;clear;close all;
sympref('FloatingPointOutput',true);

%% Known info
M0 = 2.1;                  % Free stream velocity
M4_up = 1.273;             % Velocity upstream of shock 4 (normal shockwave)
gamma = 1.4;               % specific heat ratio, assumed to be 1.4

%M4 needs to be 0.8, vary M4_up until it is - DONE
M4 = sqrt((((gamma-1)*(M4_up^2))+2)/((2*gamma*(M4_up^2))-(gamma-1)));

%% Unknowns
syms M1;                         % Mach number after first oblique shockwave
syms M2;                         % Mach number after second oblique shockwave
syms M3;                         % Mach number after third oblique shockwave
syms theta1;                     % Shock wave angle of first shockwave
syms theta2;                     % Shock wave angle of second shockwave
syms theta3;                     % Shock wave angle of third shockwave
syms d1;                         % deflection angle (half ramp angle) of M1 downstream flow
syms d2;                         % deflection angle (half ramp angle) of M2 downstream flow
syms d3;                         % deflection angle (half ramp angle) of M3 downstream flow

%% Equations for solving Mach numbers and angles


eqn_M1 = sqrt(((36*M0^4*sind(theta1)^2)-5*(M0^2*sind(theta1)^2-1)*(7*M0^2*sind(theta1)^2+5))/((7*M0^2*sind(theta1)^2-1)*(M0^2*sind(theta1)^2+5)));
eqn_d1 = acotd(tand(theta1)*((6*M0^2)/(5*(M0^2*sind(theta1)^2-1))-1));
eqn_M2 = sqrt(((36*M1^4*sind(theta2)^2)-5*(M1^2*sind(theta2)^2-1)*(7*M1^2*sind(theta2)^2+5))/((7*M1^2*sind(theta2)^2-1)*(M1^2*sind(theta2)^2+5)));
eqn_d2 = acotd(tand(theta2)*((6*M1^2)/(5*(M1^2*sind(theta2)^2-1))-1));
eqn_M3 = sqrt(((36*M2^4*sind(theta3)^2)-5*(M2^2*sind(theta3)^2-1)*(7*M2^2*sind(theta3)^2+5))/((7*M2^2*sind(theta3)^2-1)*(M2^2*sind(theta3)^2+5)));
eqn_d3 = acotd(tand(theta3)*((6*M2^2)/(5*(M2^2*sind(theta3)^2-1))-1));

eqn_7 = (M0*sind(theta1)) == (M1*sind(theta2));
eqn_8 = (M1*sind(theta2)) == (M2*sind(theta3));
eqn_9 = M3 == M4_up;

eqns = [eqn_M1, eqn_d1, eqn_M2, eqn_d2, eqn_M3, eqn_d3, eqn_7,eqn_8,eqn_9];

Solution = solve(eqns,[M1, M2, M3, theta1, theta2, theta3, d1, d2, d3])