%FINAL PROJECT #9--------------COMPUTATION OF LINEARIZED EULER EQUATIONS
%Lance Simms 12/2005
%
%u=velocity in x direction      v=velocity in y direction
%Flow is flowing in from left.  Assumed to be uniform rho=1, u=M_inf, v=0
%Linearization about p=1, u=M_inf, v=0
%
%Thin airfoil present having shape given by
%   ywall=Tx(1-x)/2     0<x<1
%   ywall=0             x<0 x>1
%
%Euler Equations are written as 
%
%   DQ       DQ     DQ
%   --  +  A -- + B --  = 0
%   Dt       Dx     Dy
%
%  where:     |     0         1      0   |      |  rho  |
%         A = | -M_inf^2   2 M_inf   0   |, Q = | rho u |
%             |     0         0    M_inf |      | rho v |
%
%             |     0         0      1    |
%         B = |     0         0     M_inf |
%             |     1         0      0    |
%  Inputs:
%         M_inf :  Free Stream Mach number
%         nx    :  Number of points on airfoil (x=0,1)
%         cfl   :  dt = min(dx,dy)*cfl/(M_inf+1 + 1)
%                     where Spectral radius of A is (M_inf+1) and B is 1
%         nmax  : Number of time steps to go
%         nout  : Output frequency of plot of cp  (0 for no plots)
%         meth  : 1: Explicit Euler  / 2: RK4 / 3: Implicit Euler
%         order : 0 : central 2nd / 1 : 1st order upwind / 2 : 2nd order upwind
%         ichar : 0 : simple bc / 1 : char bc / -1 simpler bc
%
%
clear;

%MACH NUMBER

M_inf=.8;
M_inf=input('Enter M_inf');

%Grid Generation
xmin=-1.0;          xmax=3.0;
ymin=0.0;           ymax=2.0;

nx=10;
nx=input('Enter number of points on body:');

dx=1.0/(nx-1);      dy=dx; 
x=[xmin:dx:xmax]';  y=[ymin:dy:ymax]';
jmax=size(x,1);     kmax=size(y,1);

%Generate meshgrid having x values as rows and y values as columns
[X,Y]=meshgrid(x,y);

jle=nx; jte=2*nx-1;
fprintf('Leading edge point= %g trailing edge point = %g \n', jle,jte);
fprintf('X of Leading edge point =%g    X of trailing edge point = %g \n', x(jle),x(jte));

%Surface Definition: tau is the thickness
tau=0.1;   
%tau=input ('tau);
s=[jle:jte];
ys=zeros(jmax,1);   ys(s)=tau*0.5*x(s).*(1.0-x(s));

%Try to figure out indices of airfoil
foil=0;
for xind=1:length(x)
    newfoil=find(y<ys(xind));
    foil=[foil,xind*length(y)+newfoil'];
end
foil=foil(2:length(foil));

%Initialze elements of Q
Q(:,:,1)=ones(jmax,kmax,1);         %rho
Q(:,:,2)=M_inf*ones(jmax,kmax,1);   %rho u
Q(:,:,3)=zeros(jmax,kmax,1);        %rho v

%Initialize storage for Runge-Kutta Methods
Q1=zeros(jmax,kmax,3);
Q2=zeros(jmax,kmax,3);
Q3=zeros(jmax,kmax,3);

%Jacobian Matrices
A=zeros(3,3);                                       B=zeros(3,3);

A(1,:)=[0,          1,           0    ];            B(1,:)=[0,   0,   1    ];
A(2,:)=[-M_inf^2+1  2.0*M_inf,   0    ];            B(2,:)=[0,   0,   M_inf];
A(3,:)=[0,          0,           M_inf];            B(3,:)=[1,   0,   0    ];

%FLUX SPLIT
[evecA,evalA]=eig(A);                 [evecB,evalB]=eig(B);
evalAp=(evalA+abs(evalA))/2;          evalBp=(evalB+abs(evalB))/2;
evalAm=(evalA-abs(evalA))/2;          evalBm=(evalB-abs(evalB))/2;

Ap=evecA*evalAp*inv(evecA);           Bp=evecB*evalBp*inv(evecB);
Am=evecA*evalAm*inv(evecA);           Bm=evecB*evalBm*inv(evecB);

%Thin airfoil BC v=M_inf D(ys)/Dx
Q(s,1,3)=M_inf*(ys(s+1)-ys(s-1))*0.5/dx;

%Indices
J=[1:jmax];         K=[1:kmax];
JM=[2:jmax-1];      KM=[2:kmax-1];

%INPUTS
nmax=300;
nmax=input('Enter nmax ');
nmax_tot=nmax;

%cfl=0.5;
cfl=input('Enter cfl');

meth=0;
meth=input('Enter 0 for RK4 \n      1 for Euler Explicit\n      2 for Euler implicit: ');

order=1;    
epse=0.0;
%epse=Enter('2nd order diss coef epse   ',epse);
Eig_t=2.0+M_inf;

dt=min(dx,dy)*cfl/Eig_t;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%Eventually placed in separate function

R = rhs_cen(Q,JM,KM,dx,dy,A,B,epse);

ncount=0;
cpu0=cputime;

figure;
set (gcf,'doublebuffer','on');
im=imagesc(x,y,Q(:,:,1));
maxp=text(.1,1.05,'maxp=','units','normalized');
maxpu=text(.3,1.05,'maxpu=','units','normalized');
maxpv=text(.5,1.05,'maxpv=','units','normalized');

if (meth==0)
for nstep=1:nmax_tot
    ncount=ncount+1;
    
    Q1=zeros(jmax,kmax,3);
    Q2=zeros(jmax,kmax,3);
    Q3=zeros(jmax,kmax,3);

    %First RK4 Step:        Q'n+1/2=Qn+.5hR(Q)n
    R1=rhs_cen(Q1,JM,KM,dx,dy,A,B,epse);
    Q1=Q+0.5*dt*R1;
    
    %Second RK4 Step:       Q-n+1/2=Qn+.5hR(Q')n
    R2=rhs_cen(Q2,JM,KM,dx,dy,A,B,epse);
    Q2=Q+0.5*dt*R2;
    
    %Third RK4 Step:        Q/n+1=Qn+hR(Q-)n+1/2
    R3=rhs_cen(Q3,JM,KM,dx,dy,A,B,epse);
    Q3=Q+dt*R3;
    
    %Last RK4 Step:         Qn=1/6*h*[R(Q)n+2(R(Q')n+1/2+R(Q-)n+1/2)+R(Q)/n+1]
    R4=rhs_cen(Q3,JM,KM,dx,dy,A,B,epse);

    Q(JM,KM,:) = Q(JM,KM,:) + (1/6)*dt*(R1(JM,KM,:)+2*(R2(JM,KM,:)+R3(JM,KM,:))+R4(JM,KM,:));
    
    %Apply Boundary conditions
    Q=bc8(Q,JM,KM,jmax,kmax,M_inf);
    
    PlotQ=Q(:,:,1)';
    PlotQ(foil)=1;
    %subplot(2,1,2);
   
    set(gcf,'CurrentObject',im);    
    %imagesc(x,y,flipud(PlotQ));
    %imagesc(x,y(kmax:-1:1),Q(:,:,1)');
    set(gcf,'units','normalized');
    set(im,'cdata',flipud(PlotQ),'xdata',x,'ydata',y);
    set(maxp,'string',sprintf('max_p=%2.2f',max(max(Q(:,:,1)))));
    set(maxpu,'string',sprintf('max_pu=%2.2f',max(max(Q(:,:,2)))));
    set(maxpv,'string',sprintf('max_pv=%2.2f',max(max(Q(:,:,3)))));

    refresh;
    drawnow;
    %[flag]=figflag(gcf,0);
    pause(0.1);
    end
end

%%EULER EXPLICIT
if (meth==1)
 
    for nstep=1:nmax_tot

    R=rhs_flux_split_so(Q,J,K,jmax,kmax,dx,dy,Ap,Am,Bp,Bm,epse);
    Q(JM,KM,:) = Q(JM,KM,:) + dt*R(JM,KM,:);
    
    %Boundary Conditions
    Q=bc8(Q,JM,KM,jmax,kmax,M_inf);
    
    %Position Airfoil in Plot
    PlotQ=Q(:,:,1)';
    PlotQ(foil)=1;
   
    imagesc(x,y,flipud(PlotQ));
    %imagesc(x,y(kmax:-1:1),Q(:,:,1)');
    pause(0.1);
    end
end

if (meth==1)
 
    for nstep=1:nmax_tot

    R=rhs_cen(Q,JM,KM,dx,dy,A,B,epse);
    Q(JM,KM,:) = Q(JM,KM,:) + dt*R(JM,KM,:);
    
    %Boundary Conditions
    Q=bc8(Q,JM,KM,jmax,kmax,M_inf);
    
    %contourf(X,Y,Q(:,:,1)',50);
    %xlabel('rho');
    PlotQ=Q(:,:,1)';
    PlotQ(foil)=0;
    set(gcf,'CurrentObject',im);    
    set(im,'cdata',flipud(PlotQ),'xdata',x,'ydata',y);
    set(maxp,'string',sprintf('max_p=%2.2f',max(max(Q(:,:,1)))));
    % set(gca,'title',text('string',sprintf('b=%2.2f',max(max(max(Q(:,:,1)))),...
    %  'color','r','fontsize',14)));
    drawnow;
    pause(0.1);
    end
end