clearglobal;
clear;
xdel(winsid())
clc;
disp('AERODYNAMIC CHARACTERISTICS OF 2D SYMMETRICAL DOUBLE WEDGE AEROFOIL');
disp('*******************************************************************');
disp('------------------CODE DEVELOPED BY VIJAY GOPAL--------------------');
disp('*******************************************************************');
disp('** Assumptions **');
disp(' > Inviscid Flow ');
disp(' > Calorifically Perfect Gas ');
disp(' > No Detached Shock Occurance');
disp(' > The Free-Stream Flow is Supersonic and Directed Along Right-Side');
disp('*******************************************************************');
disp('> Enter the values of the quantities prompted ');
disp('');
eh=input('* Enter the value of half angle of symmetrical double-wedge aerofoil in degrees : ');
t_half=input('* Enter the value of half thickness in meters: ');
M=input('* Enter the value of free-stream Mach number (>1) : ');
alpha_max=input('* Enter the value of maximum angle of attack in degrees (<46): ');
P=input('* Enter the value of free-stream static pressure in pascals : ');
r=input('* Enter the value ratio of specific heats (>0): ');
R=input('* Enter the value of specific gas constant in J/(Kg.K) (>0) : ');
Te=input('* Enter the value of free-stream temperature in Kelvin (>0) : ');
Interpolation=input('* Enter the number of interpolations (>= 100): ');
Interpol=int(Interpolation);
disp('VALIDATING DATA....')
a=linspace(0,alpha_max,Interpol);
t=2*t_half;
c=t./tand(eh);
e=2.*eh;
rho=((P)./(R.*Te));
cs=sqrt(r.*R.*Te);
l=0.5*(sqrt(t.^2+c.^2));
Beta_check=linspace(0.0001,90,5*Interpol);
Theta_check=[];
for i=1:5*Interpol
    Theta_check(i)=atand((2).*cotd(Beta_check(i)).*((M.*sind(Beta_check(i))).^2-1)./((r+cosd(2.*Beta_check(i))).*(M.^2)+2));
end
Theta_max=max((Theta_check));
a_max=(abs(eh)+abs(max(a)));
if(a_max>Theta_max | M<1 | Interpol<100 | alpha_max>90 |Te<0 | R<0 | r<0)
    disp('.')
    disp('.')
    disp('Detached shock detected or invalid data entry: Program aborted !!!');
    abort
else
       disp('.')
       disp('VALID DATA ENTERED')
       disp('>START COMPUTATION')
       T2=[];
       Te2=[];
       rho2=[];
       B2=[];
       P2=[];
       M2=[];
       T3=[];
       rho3=[];
       Te3=[];
       B3=[];
       P3=[];
       M3=[];
       T4=[];
       Te4=[];
       rho4=[];
       M4=[];
       P4=[];
       T5=[];
       Te5=[];
       rho5=[];
       M5=[];
       P5=[];
       for i=1:Interpol
          T2(i)=abs(eh)-abs(a(i));
          T3(i)=abs(eh)+abs(a(i));
       end
          T4=2.*abs(eh);
          T5=2.*abs(eh);
          disp('FORE UPPER REGION COMPUTATION.........')
          disp('*');
          for i=1:Interpol
                    if (T2(i)>0) then
                        //disp('SHOCK WAVE FORMED ON FORE UPPER REGION ')
                         A=((1+((r-1)/2).*(M.^2)).*tand(T2(i)));
                         B=(1-(M.^2));
                         C=((1+((r+1)/2).*(M.^2)).*tand(T2(i)));
                         D=1;
                         betacubic=poly([D,C,B,A],'tanB','c');
                         tanB=roots(betacubic);
                         tanB_strong=max(real(tanB));
                         garbage=min(real(tanB));
                         for k=1:3
                                if(real(tanB(k))>garbage & real(tanB(k))<tanB_strong) then
                                     B2(i)=atand(real(tanB(k)));
                                     P2(i)=P.*(1+((2*r)./(r+1)).*((M.*sind(B2(i))).^2-1));
                                     M2(i)=(sqrt(((M.*sind(B2(i))).^2+((2)./(r-1)))./(((((2*r)./(r-1))).*(((M.*sind(B2(i)))).^2)-1))))./(sind(B2(i)-T2(i)));
                                     rho2(i)=rho.*(((r+1).*((M.*sind(B2(i))).^2))./(((r-1).*((M.*sind(B2(i))).^2))+2));
                                     Te2(i)=Te.*((P2(i)./P).*(rho./rho2(i)));
                                end
                         end
                         A=0;
                         B=0;
                         C=0;
                         D=0;
                         tanB=[0;0;0];
                     else
                      //disp('EXPANSION WAVE FORMED ON FORE UPPER REGION ')
                      v_M=sqrt((r+1)./(r-1)).*(atan(sqrt(((r-1)./(r+1)).*(M.^2-1))))-atan(sqrt(M.^2-1));
                      T_2exp(i)=abs(T2(i));
                      deff('[func]=f(M_2)','func=T_2exp(i).*(%pi./180)+v_M-(sqrt((r+1)./(r-1)).*(atan(sqrt(((r-1)./(r+1)).*(M_2.^2-1))))-atan(sqrt(M_2.^2-1)))');
                      [M_2]=fsolve(2,f);
                      M2(i)=M_2(1);
                      P2(i)=P.*((1+((r-1)/2).*(M.^2))./(1+((r-1)/2).*(M2(i).^2))).^(r./(r-1));
                      Te2(i)=Te.*((1+((r-1)/2).*(M.^2))./(1+((r-1)/2).*(M2(i).^2)));
                      rho2(i)=rho.*(((P2(i)./P).*(Te./Te2(i))));
                      M_2=[0];
                      v_M=[0]
              end
          end
          disp('FORE LOWER REGION COMPUTATION..........')
          disp('*');
          disp('SHOCK WAVE FORMED ON FORE LOWER REGION ')
          for i=1:Interpol
               A=((1+((r-1)/2).*(M.^2)).*tand(T3(i)));
               B=(1-(M.^2));
               C=((1+((r+1)/2).*(M.^2)).*tand(T3(i)));
               D=1;
               betacubic=poly([D,C,B,A],'tanB','c');
               tanB=roots(betacubic);
               tanB_strong=max(real(tanB));
               garbage=min(real(tanB));
               for k=1:3
                   if(real(tanB(k))>garbage & real(tanB(k))<tanB_strong) then
                        B3(i)=atand(real(tanB(k)));
                        P3(i)=P.*(1+((2*r)./(r+1)).*((M.*sind(B3(i))).^2-1));
                        M3(i)=(sqrt(((M.*sind(B3(i))).^2+((2)./(r-1)))./(((((2*r)./(r-1))).*(((M.*sind(B3(i)))).^2)-1))))./(sind(B3(i)-T3(i)));
                        rho3(i)=rho.*(((r+1).*((M.*sind(B3(i))).^2))./(((r-1).*((M.*sind(B3(i))).^2))+2));
                        Te3(i)=Te.*((P3(i)./P).*(rho./rho3(i)));
                   end
               end
                         A=0;
                         B=0;
                         C=0;
                         D=0;
                         tanB=[0;0;0];
          end
          disp('AFT UPPER REGION COMPUTATION...............')
          disp('*');
          disp('EXPANSION WAVES FORMED ON AFT UPPER REGION ')
          for i=1:Interpol
          funcprot(0)
          v_M4(i)=sqrt((r+1)./(r-1)).*(atan(sqrt(((r-1)./(r+1)).*(M2(i).^2-1))))-atan(sqrt(M2(i).^2-1));
          deff('[func]=f(M_4)','func=T4.*(%pi./180)+v_M4(i)-(sqrt((r+1)./(r-1)).*(atan(sqrt(((r-1)./(r+1)).*(M_4.^2-1))))-atan(sqrt(M_4.^2-1)))');
          [M_4]=fsolve(2,f);
          M4(i)=M_4(1);
          P4(i)=P2(i).*(((1+((r-1)/2).*(M2(i).^2))./(1+((r-1)/2).*(M4(i).^2))).^(r./(r-1)));
          Te4(i)=Te2(i).*((1+((r-1)/2).*(M2(i).^2))./(1+((r-1)/2).*(M4(i).^2)));
          rho4(i)=rho2(i).*((P4(i)./P2(i)).*(Te2(i)./Te4(i)));
          M_4=[0];
          v_M4=[0];
          end
          disp('AFT LOWER REGION COMPUTATION..............')
          disp('*');
          disp('EXPANSION WAVES FORMED ON AFT UPPER REGION ')
          for i=1:Interpol
            funcprot(0)
            v_M5(i)=sqrt((r+1)./(r-1)).*(atan(sqrt(((r-1)./(r+1)).*(M3(i).^2-1))))-atan(sqrt(M3(i).^2-1));
            deff('[func]=f(M_5)','func=T5.*(%pi./180)+v_M5(i)-(sqrt((r+1)./(r-1)).*(atan(sqrt(((r-1)./(r+1)).*(M_5.^2-1))))-atan(sqrt(M_5.^2-1)))');
            [M_5]=fsolve(2,f);
            M5(i)=M_5(1);
            P5(i)=P3(i).*(((1+((r-1)./2).*(M3(i).^2))./(1+((r-1)./2).*(M5(i).^2))).^(r./(r-1)));
            Te5(i)=Te3(i).*((1+((r-1)./2).*(M3(i).^2))./(1+((r-1)./2).*(M5(i).^2)));
            rho5(i)=rho3(i).*(((P5(i)./P3(i)).*(Te3(i)./Te5(i))));
            M_5=[0];
            v_M5=[0];
          end
          disp('COMPUTING LIFT AND DRAG...............')
          for i=1:Interpol
            Lift(i)=l.*((P5(i)-P2(i)).*cosd(abs(eh)-abs(a(i)))+(P3(i)-P4(i)).*cosd(abs(eh)+abs(a(i))));
            Drag(i)=l.*((-P5(i)+P2(i)).*sind(abs(eh)-abs(a(i)))+(P3(i)-P4(i)).*sind(abs(eh)+abs(a(i))));
            V=M.*cs;
            Moment(i)=(((c/2)^2)/2)*((P3(i)-P2(i))+(P4(i)-P5(i))+((P2(i)-P3(i))+(P5(i)-P4(i)))*((tand(eh))^2))
            dp=(0.5).*((P)./(R.*Te)).*(V.^2);
            CL(i)=(Lift(i))./(dp.*c);
            CD(i)=(Drag(i))./(dp.*c);
            CM(i)=(Moment(i))./(dp*c*c);
          end
end
disp('XXXXXXXXXXXXXXXXXXXXXXXXX   END COMPUTATION    XXXXXXXXXXXXXXXXXXXXXXXX');
disp('')
disp('----------------------- DEFAULT  POST PROCESSING  ---------------------');
disp('Plotting Results....!');
z=strcat(["Number of Interpolations = ",string(Interpol)])
figure(1);
plot2d(a,CD,[2],nax=[10,10,10,10],leg=z);
xlabel('Angle of attack');
ylabel('Section Drag Coefficient');
title('Plot of Section Drag Coefficient V/S Angle of Attack','color','blue');  
figure(2);
plot2d(a,CL,[2],nax=[10,10,10,10],leg=z);
xlabel('Angle of attack');
ylabel('Section Lift Coefficient');  
title('Plot of Section Lift Coefficient V/S Angle of Attack','color','blue');  
figure(3);
plot2d(a,Drag,[2],nax=[10,10,10,10],leg=z);
xlabel('Angle of attack');
ylabel('Drag force per unit span(N/m)');
title('Plot of Drag Force Per Unit Span (N/m) V/S Angle of Attack','color','blue');    
figure(4);
plot2d(a,Lift,[2],nax=[10,10,10,10],leg=z);
xlabel('Angle of attack');
ylabel('Lift Force per unit span (N/m)');
title('Plot of Lift Force Per Unit Span (N/m) V/S Angle of Attack','color','blue');  
figure(5);
plot2d(a,Moment,[2],nax=[10,10,10,10],leg=z);
xlabel('Angle of attack');
ylabel('Moment per unit span about geometrical center (N)');
title('Plot of Moment Per Unit Span About the Geometrical Center (N) V/S Angle of Attack','color','blue');  
figure(6);
plot2d(a,CM,[2],nax=[10,10,10,10],leg=z);
xlabel('Angle of attack');
ylabel('Coefficent of Moment per unit span');
title('Plot of Coefficient of Moment per Unit Span About the Geometrical Center V/S Angle of Attack','color','blue');
disp('------------------------ USER  POST PROCESSING  --------------------------');
decide=1;
while (decide==1)
for lx=1:10000
    alpha=0;
    alpha=input('* Enter the desired angle of attack in degrees (Expected less than maximum angle of attack) : ');
    Mach_Number_FU=0;
    Mach_Number_FL=0;
    Mach_Number_AU=0;
    Mach_Number_AL=0;
    Temperature_Kelvin_FU=0;
    Temperature_Kelvin_FL=0;
    Temperature_Kelvin_AU=0;
    Temperature_Kelvin_AU=0;
    Pressure_Pascals_FU=0;
    Pressure_Pascals_FL=0;
    Pressure_Pascals_AU=0;
    Pressure_Pascals_AL=0;
    fore=[0];
    aft=[0];
    ch=0;
    CP_l=0;
    CP_u=0;
    i1=0;
    z_alpha=string(alpha);
    z_AOA=strcat(["Angle of Attack = ",z_alpha," degrees"]);
    i1=int((Interpol/alpha_max).*alpha);
    if i1==0 then
        i1=1;
    end
    if alpha>alpha_max then 
        disp('Invalid Angle of attack, POST PROCESSING ABORTED')
        disp('To repeat evaluation for different angle of attack enter ''resume'' or to end program enter ''abort'' ');
        pause
        break
    else
    disp('****************Aerofoils Fore Upper Region Data****************** ')
    disp(z_AOA);
    Mach_Number_FU=M2(i1);
    disp(Mach_Number_FU,'Mach Number');
    Temperature_Kelvin_FU=Te2(i1);
    disp(Temperature_Kelvin_FU,'Temperature in Kelvin');
    Pressure_Pascals_FU=P2(i1);
    disp(Pressure_Pascals_FU,'Pressure in pascals')
    disp(((Pressure_Pascals_FU)./(R.*Temperature_Kelvin_FU)),'Density in Kg/m^3');
    disp('****************Aerofoils Fore Lower Region Data*******************')
    disp(z_AOA);
    Mach_Number_FL=M3(i1);
    disp(Mach_Number_FL,'Mach Number')
    Temperature_Kelvin_FL=Te3(i1);
    disp(Temperature_Kelvin_FL,'Temperature in Kelvin')
    Pressure_Pascals_FL=P3(i1);
    disp(Pressure_Pascals_FL,'Pressure in pascals');
    disp(((Pressure_Pascals_FL)./(R.*Temperature_Kelvin_FL)),'Density in Kg/m^3');
    disp('****************Aerofoils Aft Upper Region Data********************')
    disp(z_AOA);
    Mach_Number_AU=M4(i1);
    disp(Mach_Number_AU,'Mach Number')
    Temperature_Kelvin_AU=Te4(i1);
    disp(Temperature_Kelvin_AU,'Temperature in Kelvin')
    Pressure_Pascals_AU=P4(i1);
    disp(Pressure_Pascals_AU,'Pressure in pascals')
    disp(((Pressure_Pascals_AU)./(R.*Temperature_Kelvin_AU)),'Density in Kg/m^3');
    disp('****************Aerofoils AFT Lower Region Data********************')
    disp(z_AOA);
    Mach_Number_AL=M5(i1);
    disp(Mach_Number_AL,'Mach Number')
    Temperature_Kelvin_AL=Te5(i1);
    disp(Temperature_Kelvin_AL,'Temperature in Kelvin')
    Pressure_Pascals_AL=P5(i1);
    disp(Pressure_Pascals_AL,'Pressure in Pascals')
    disp(((Pressure_Pascals_AL)./(R.*Temperature_Kelvin_AL)),'Density in Kg/m^3');
    disp('****************General Aerfoil data******************************')
    disp(z_AOA);
    Coefficient_Lift=CL(i1);
    disp(Coefficient_Lift,'Section Lift Coefficient')
    Coefficient_Drag=CD(i1);
    disp(Coefficient_Drag,'Section Drag Coefficient')
    Coefficient_Moment=CM(i1);
    disp(Coefficient_Moment,'Section Moment Coefficient')
    Lift_Force_Newton=Lift(i1);
    disp(Lift_Force_Newton,'Lift Force per Unit Span [N/m]')
    Drag_Force_Newton=Drag(i1);
    disp(Drag_Force_Newton,'Drag Force Per Unit Span [N/m]')
    Moment_Newton_perMeter=Moment(i1);
    disp(Moment_Newton_perMeter,'Moment About the Geometrical Center Per Unit Span [N]')
    disp('Plotting Coefficient of pressure curves')
    fore=[ones(1,5),zeros(1,5)];
    aft=[zeros(1,5),ones(1,5)];
    CP_u=((P2(i1)-P).*fore+(P4(i1)-P).*aft)./dp;
    CP_l=((P3(i1)-P).*fore+(P5(i1)-P).*aft)./dp;
    ch=(c./10):(c./10):(c);
    zero=zeros(1,10);
    Z_A=strcat(["Plot of Coefficient of Pressure v/s Chord Length at Angle of Attack = ",z_alpha," degrees"]);
    figure
    plot2d(ch,CP_u,[1],nax=[10,10,10,10])
    plot2d(ch,CP_l,[2],nax=[10,10,10,10])
    plot2d(ch,zero,[5],nax=[10,10,10,10])
    legends(['Upper Surface';'Lower Surface';'Free-Stream'],[1,2,5],opt="ur")
    xlabel('Chord')
    ylabel('Coefficient of Pressure')  
    title([Z_A],'color','blue')
    end
    disp('To repeat evaluation for different angle of attack enter ''resume'' or to end program enter ''abort'' ');
    pause
    disp('******************REPEATING POST PROCESSING****************************')
  end
disp('******************REPEATING POST PROCESSING****************************')
decide=1;
end



