%%This script will attempt to solve the Transonic Disturbance equations for
%%a thin airfoil.
%%
M_inf=0;
M_inf=input('Enter M_inf (Mach_numer)   ');
V_inf=M_inf;

%Ratio of Specific Heats
gamma=1.4;

%%GRID INDEX DEFINITIONS
I=2;                      J=2;               %Beginning Indices
IL=51;                    JL=51;             %Maximum Indices   

I_lead_edge=(IL-1)*.3;           I_trail_edge=(IL-1)*.7+1;

%SEPERATE X INDICES FOR 3 REGIONS (SUBSONIC--SONIC---SUBSONIC)
I_before_foil=[2:I_lead_edge];
I_over_foil=[I_lead_edge:I_trail_edge-1];
I_after_foil=[I_trail_edge:IL];

IO=[1:IL];              JO=[1:JL];        %Full Range Indices
IM=[2:IL-1];            JM=[2:JL-1];     %Interior Indices

%DISTANCES
c=1;                                         %Chord length
Dy=50*c;                  Dx=50*c;           %Max DIistance
del_y_min=.06*c/10;       del_x_min=c/20;    %The minimum separation in y near the foil

%The Grid was solved with the exponential mesh stretching formula discussed
%in the text on page 11 of Chapter 5
%The value of k was determined 
[kx,ky]=Stretching_Function(I,J,IL,JL,Dx,Dy,del_y_min,del_x_min);

J_fac=(JO-1)/(JL-1);             I_fac=(I_before_foil-1)/(I_lead_edge-1);


%GENERATE THE Y SPACING-----Spaced with the stretchin function
y=Dy*(exp(ky*J_fac)-1)/(exp(ky)-1);

%GENERATE THE X SPACING-----Spaced with the stretchin function

I_before_foil=[2:I_lead_edge];
I_over_foil=[I_lead_edge:I_trail_edge-1];
I_after_foil=[I_trail_edge:IL];

X_trail=Dx*(exp(kx*I_fac)-1)/(exp(kx)-1)+1;
X_over=(0:c/20:1);
X_lead=-fliplr(X_trail)+1;
x=[X_lead,X_over,X_trail];

%GENERATE THE MESH
[X,Y]=meshgrid(x,y);

%DRAW IT IF YOU"D LIKE
%figure;                 set(gca,'XLIM',[min(x),max(x)]);
%set(gca,'XTICK',x);     set(gca,'YTICK',y);                 grid;

%PUT THE AIRFOIL AT Y=0
tau=0.06
ys=zeros(1,IL);         ys(I_over_foil)=tau*.5*x(I_over_foil).*(1.0-x(I_over_foil));
dys=zeros(1,IL);        dys(I_over_foil)=tau*.5*(1-2*x(I_over_foil));

%SET UP THE INITIAL VALUES OF THE POTENTIAL
PHI=ones(JL,IL);
[PHI]=Murman_Cole_bc(PHI,IM,JM,IL,JL,M_inf,dys,del_y_min);

%Declare the Coefficient
Aij=zeros(JL,IL);       Ai_1j=zeros(JL,IL);

Aij(JM,IM)=1-M_inf^2-(gamma+1)*M_inf^2*(PHI(JM,IM+1)-PHI(JM,IM-1))./...
                                     (X(JM,IM+1)-X(JM,IM-1))./V_inf;
Ai_1j(JM,IM+1)=1-M_inf^2-(gamma+1)*M_inf^2*(PHI(JM,IM+1)-PHI(JM,IM-1))./...
                                           (X(JM,IM+1)-X(JM,IM-1))./V_inf;
mu=(Aij > 0);

%THE FOLLOWING USES THE COLUMN SWEEP FROM J=2 to J=JL

%%%SET UP DENOMINATOR TERMS FOR Y----ONLY NEED ONE COLUMN HERE
%%INDIVIDUAL
yj1=y(3:JL)';                                  yjm1=y(1:JL-2)';
yj0=y(2:JL-1)';
%%DENOMINATORS
y_j1_j0=(yj1-yj0).*(yj1-yjm1);                 y_j0_jm1=(yj0-yjm1).*(yj1-yjm1);

%%SET UP DENOMINATOR TERMS FOR X-----ONLY NEED ONE ROW HERE
%%INDIVIDUAL
xi1=X(2:JL-1,4:IL);                            xim2=X(2:JL-1,1:IL-3);
xi0=X(2:JL-1,3:IL-1);                          xim1=X(2:JL-1,2:IL-2);
%%DENOMINATORS
x_i1_i0_c=(xi1-xi0).*(xi1-xim1);               x_im1_im2_b=(xim1-xim2).*(xi0-xim2);
x_i0_im1_c=(xi0-xim1).*(xi1-xim1);             x_i0_im1_b= (xi0-xim1).*(xi0-xim2);

%%Y difference prefactors---These do not change-----b*phi_i_j-1--c*phi_i_j+1
b(1:JL-2)=(1./y_j0_jm1);                       b(JL-1)=0;     
c(2:JL-1)=(1./y_j1_j0);                        c(1)=-1;

%%Plotting
set (gcf,'doublebuffer','on');
im=imagesc(x,y,PHI);

num_steps=200    
%%Begin the loop
for n=1: num_steps
    for i=2: IL-1
        
        %%Case where i=2, we cannot use backward differencing.  No A_im1_j
        %%term
      if i==2
        Aij=1-M_inf^2-(gamma+1)*M_inf^2*(PHI(JM,3)-PHI(JM,1))./...
                                        (X(JM,3)-X(JM,1))./V_inf;
        mu=(Aij<0);
        x_i1_i0_c_first=(X(JM,3)-X(JM,2)).*(X(JM,3)-X(JM,1));
        x_i0_im1_c_first=(X(JM,2)-X(JM,1)).*(X(JM,3)-X(JM,1));
        
        %%LHS OF EQN---TERMS WITH PHI_i--bijPhi_i_j-1+aijPhi_i_j+cijPhi_i_j+1
        a(2:JL-1)=(1-mu).*Aij.*(-1./x_i1_i0_c_first-1./x_i0_im1_c_first)...
                  -1./y_j1_j0-1./y_j0_jm1;
        a(1)=1;             a(JL)=1;
        
        %%RHS OF EQUATION----TERMS INVOLVING PHI_i+1 and PHI_i-1
        f(2:JL-1)=-(1-mu).*Aij.*PHI(JM,3)./x_i1_i0_c_first...
                  -(1-mu).*Aij.*PHI(JM,1)./x_i0_im1_c_first;
        f(1)=-del_y_min*V_inf*dys(i);     f(JL)=1;
        
        %%Do the solver
        [L,U,alph,gamm]=LU_decomposition(a,b,c,JL);
        PHI(:,i)=LU_solve(alph,b,gamm,f,JL);
        
        %Plot
        set(im,'cdata',flipud(PHI),'xdata',x,'ydata',y);
        drawnow;

      else
        %%Now we can use backward differencing----Get on it
        Aij=1-M_inf^2-(gamma+1)*M_inf^2*(PHI(JM,i+1)-PHI(JM,i-1))./...
                                        (X(JM,i+1)-X(JM,i-1))./V_inf;
        Aim1j=1-M_inf^2-(gamma+1)*M_inf^2*(PHI(JM,i)-PHI(JM,i-2))./...
                                          (X(JM,i)-X(JM,i-2))./V_inf;
        muij=(Aij<0);
        muim1j=(Aim1j<0);
        
        %%LHS
        %%LHS OF EQN---TERMS WITH PHI_i--bijPhi_i_j-1+aijPhi_i_j+cijPhi_i_j+1
        a(2:JL-1)= (1-muij).*Aij.*(-1./x_i1_i0_c(:,i-2)-1./x_i0_im1_c(:,i-2))...
                  +(muim1j).*Aim1j.*(1./x_i0_im1_b(:,i-2))+...
                   (-1./y_j1_j0-1./y_j0_jm1);
        a(1)=1;                                         a(JL)=1;
        
         %%RHS OF EQUATION----TERMS INVOLVING PHI_i+1 and PHI_i-1
        f(2:JL-1)=-(1-muij).*Aij.*PHI(JM,i+1)./x_i1_i0_c(:,i-2)...
                  -(1-muij).*Aij.*PHI(JM,i-1)./x_i0_im1_c(:,i-2)...
                  +muim1j.*Aim1j.*(1./x_i0_im1_b(:,i-2)+1./x_im1_im2_b(:,i-2)).*PHI(JM,i-1)...
                  -muim1j.*Aim1j.*(1./x_im1_im2_b(:,i-2)).*PHI(JM,i-2);
        f(1)=-del_y_min*V_inf*dys(i);     f(JL)=1;
        
        %%Do the solver
        [L,U,alph,gamm]=LU_decomposition(a,b,c,JL);
        PHI(:,i)=LU_solve(alph,b,gamm,f,JL);
        
        %%PLOT
        set(im,'cdata',flipud(PHI),'xdata',x,'ydata',y);
        drawnow;
        
      end
      
    end
end