%%Lance Simms                                                       1/30/06

%%This script will attempt to solve the Transonic Disturbance equations for
%%a Circular Chord Thin Airfoil.

plot_density=0;
M_inf=0;
M_inf=input('Enter M_inf (Mach_numer)   ');
V_inf=M_inf;

num_steps=200    
num_steps=input('Enter Number of Steps       ');

%Ratio of Specific Heats and quantities
gamma=1.4;                          rho_inf=1;
a_inf=V_inf/M_inf;                  p_inf=rho_inf*a_inf^2/gamma;
gamma_exp=(gamma/(gamma-1));

%%GRID INDEX DEFINITIONS
I=2;                                J=2;               %Beginning Indices
IL=51;                              JL=51;             %Maximum Indices   
%Lead edge occurs at 3/10 of x and trailing edge at 7/10+1
I_lead_edge=(IL-1)*.3+1;            I_trail_edge=(IL-1)*.7+2;

%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,I_lead_edge,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;
%figure;

%PUT THE AIRFOIL AT Y=0
tau=.03;                c=1;
[xo,yo,R]=Circular_Chord(tau,c);
ys=zeros(1,IL);         ys(I_over_foil)=sqrt(R^2-(x(I_over_foil)-xo).^2)+yo;
dys=zeros(1,IL);        dys(I_over_foil)=-(x(I_over_foil)-xo)./...
                                           sqrt(R^2-(x(I_over_foil)-xo).^2);
                                       
%SET UP THE INITIAL VALUES OF THE POTENTIAL
PHI=ones(JL,IL);

%%%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)';
Yj1=Y(3:JL,3:IL-1);                             Yjm1=Y(1:JL-2,3:IL-1);
Yj0=Y(2:JL-1,3:IL-1);

%%DENOMINATORS
y_j1_j0=(yj1-yj0).*(yj1-yjm1);                  y_j0_jm1=(yj0-yjm1).*(yj1-yjm1);
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');
if plot_density==1
    subplot(3,1,3);
    im=imagesc(x,y,PHI);
end

cp_ind=[min(I_over_foil)-3:max(I_over_foil)+3];
subplot(3,1,1);
cpplot=plot(x(cp_ind),ones(size(cp_ind)))
subplot(3,1,2);
Residplot=semilogy([1],[1])
R_max_old=1;

%%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);
        
        if plot_density==1
          set(im,'cdata',flipud(PHI),'xdata',x,'ydata',y);
          drawnow;
        end
        
      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
        if plot_density==1
          set(im,'cdata',flipud(PHI),'xdata',x,'ydata',y);
          drawnow;
        end
      end
        
    end
    
    %The span in x has been covered.  Plot cp vs x
    u=V_inf+(PHI(1,cp_ind+1)-PHI(1,cp_ind-1))./...
            (x(cp_ind+1)-x(cp_ind-1));
    v=(PHI(2,cp_ind)-PHI(1,cp_ind))./...
      (del_y_min);
    vel=(u.^2+v.^2)./V_inf^2;
    p=p_inf*(1-(gamma-1)*(M_inf^2.*(vel-1))/2).^gamma_exp;
    neg_cp=-(p-p_inf)/(.5*rho_inf*V_inf^2);
    set(cpplot,'ydata',neg_cp);
    drawnow;
   
    %Get the maximum Residual Rij
    Aij=1-M_inf^2-(gamma+1)*(M_inf^2).*(PHI(JM,4:IL)-PHI(JM,2:IL-2))./...
                                       (X(JM,4:IL)-X(JM,2:JL-2))./V_inf;
    Aim1j=1-M_inf^2-(gamma+1)*(M_inf^2).*(PHI(JM,3:IL-1)-PHI(JM,1:IL-3))./...
                                       (X(JM,3:IL-1)-X(JM,1:IL-3))./V_inf;
    muij=(Aij<0);
    muim1j=(Aim1j<0);
    R=(1-muij).*Aij.*((PHI(JM,4:IL)-PHI(JM,3:IL-1))./(x_i1_i0_c)-...
                      (PHI(JM,3:IL-1)-PHI(JM,2:IL-2))./(x_i0_im1_c))...
     +muim1j.*Aim1j.*((PHI(JM,3:IL-1)-PHI(JM,2:IL-2))./(x_i0_im1_b)-...
                      (PHI(JM,2:IL-2)-PHI(JM,1:IL-3))./(x_im1_im2_b))...
                    +((PHI(3:JL,3:IL-1)-PHI(2:JL-1,3:IL-1))./(Y_j1_j0)-...
                      (PHI(2:JL-1,3:IL-1)-PHI(1:JL-2,3:IL-1))./(Y_j0_jm1));
   
    R_max=max(max(abs(R)))
    if n==1
       R_max_arr=[R_max];
       n_arr=[n]
    else
       R_max_arr=[R_max_arr,R_max];
       n_arr=[n_arr,n];
    end
    set(Residplot,'xdata',n_arr,'ydata',R_max_arr);
              
end
            
%%Generate Final Figures to turn in 
%%CP PLOT
figure;
subplot(2,1,1);
cpplot=plot(x(cp_ind),neg_cp)
xlabel('x/c'); ylabel('-Cp');
title(sprintf('Cp vs x for flow with M_{inf}=%2.3f',M_inf));
%%Residual PLOT
subplot(2,1,2);
semilogy(n_arr,R_max_arr);
xlabel('Step Number'); ylabel('Residual');
Title(sprintf('Maximum residual for M_{inf}=%2.3f',M_inf));

