clc;clear all;
%% Problem 4

%Nomenclature for numbering:
%No number means freestream
%1 means upper left panel
%2 means upper right panel
%3 means bottom panel

M = 3;
gamma = 1.4;
alpha = 5; %degrees
epsilon = 3; %Angle the panels of the triangle, degrees

%% Calculate P1/P using expansion wave relations
theta_R1 = alpha - epsilon;

%Using Appendix C:
[~,v] = flowprandtlmeyer(gamma,M); %Prandtl-Meyer for Mach 3
v1 = v + theta_R1; %New Prandtl - Meyer

%Mach in Region 2
M1 = flowprandtlmeyer(gamma,v1,'nu');

%P/P0
[~,~,stag_pressure_ratio] = flowisentropic(gamma,M);

%P1/P01
[~,~,stag_pressure_ratio_1] = flowisentropic(gamma,M1);

%Answer
pressure_ratio_R1 = stag_pressure_ratio_1/stag_pressure_ratio; %P1/P, region 1

%% Calculate P2/P using expansion wave relations
theta_R2 = 2*epsilon;
v2 = theta_R2 + v1; %degrees, turn once at tip, turn again for panel

%Mach in Region 2
M2 = flowprandtlmeyer(gamma,v2,'nu');

%P2/P02
[~,~,stag_pressure_ratio_2] = flowisentropic(gamma,M2);

%ANSWER
pressure_ratio_R2 = pressure_ratio_R1*(stag_pressure_ratio_2/stag_pressure_ratio_1); %P2/P, accounts for effects of Region 1.

%% Calculate P3/P1 using oblique shock relations
theta_R3 = alpha; %degrees, total deflection from zero aoa

lambda = sqrt((M^2-1)^2 - 3*(1+((gamma-1)/2)*M^2)*(1+((gamma+1)/2)*M^2)*tand(theta_R3)^2);
X = ((M^2-1)^3-9*((1+((gamma-1)/2)*M^2)*(1 + ((gamma-1)/2)*M^2 + ((gamma+1)/4)*M^4)*tand(theta_R3)^2))/lambda^3;
delta = 1;

%Calculate Shock Angles, degrees
beta = (atan((M^2-1+2*lambda*cos((4*pi*delta + acos(X))/3))/(3*(1 + ((gamma-1)/2)*M^2)*tand(theta_R3))))*180/pi;

Mn = M*sind(beta);

%Mach in Region 3
Mn3 = sqrt((1 + ((gamma-1)/2)*Mn^2)/(gamma*Mn^2 - (gamma-1)/2));
M3 = Mn3/sind(beta-theta_R3);


%ANSWER
pressure_ratio_R3 = 1 + ((2*gamma)/(gamma+1))*(Mn^2-1); %P3/P

%% Calculate Cp
Cp1 = (2/(gamma*M^2))*(pressure_ratio_R1 - 1);
Cp2 = (2/(gamma*M^2))*(pressure_ratio_R2 - 1);
Cp3 = (2/(gamma*M^2))*(pressure_ratio_R3 - 1);

disp(['Cp of Panel 1 is: ',num2str(Cp1)])
disp(['Cp of Panel 2 is: ',num2str(Cp2)])
disp(['Cp of Panel 3 is: ',num2str(Cp3) newline])

%% Calculate Cl and Cd
Cn = .5*((Cp3 - Cp1)+(Cp3 - Cp2));
Ca = .5*tand(epsilon)*(Cp1 - Cp2);

Cl = (Cn*cosd(alpha)) - (Ca*sind(alpha));
Cd = (Cn*sind(alpha)) + (Ca*cosd(alpha));

disp(['Cl is: ',num2str(Cl)])
disp(['Cd is: ',num2str(Cd)])