% Euler V1.0
% THP 11/2K
% Linearized Euler Code for A Thin Circular Arc Airfoil
%
% Flow is Assumed to be rho = 1, u = M_inf , v = 0 at Inflow
% Linearization is about rho = 1, u = M_inf, v = 0
% Thin Airfoil BC 
%
% Euler equations are written as 
%    ( Pressure = rho is used i.e. Constant Temperature)
%
%   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;

M_inf = 0.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);
[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   = 0.02; 
tau = 0.1;
%tau = input('tau ');
ys    = zeros(jmax,1);
s     = [jle:jte];
ys(s) = tau*0.5*x(s).*(1.0-x(s));

% Initialize Q's
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

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

A(1,2) = 1.0;
A(2,1) = -M_inf^2 + 1;
A(2,2) = 2.0*M_inf;
A(3,3) = M_inf;

B(1,3) = 1;
B(2,3) = M_inf;
B(3,1) = 1;

%  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 = 100;
nmax = input('Enter nmax ');
nmax_tot = nmax;


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

meth = 0;
meth = input('Enter 0 for Euler Implicit AF, 1 for Euler Explicit:  ');

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;
fprintf(' AT CFL = %g  dt = %g \n',cfl,dt);

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

% Initial Residual Computation
resid_0 = (norm(R(:,:,1)) + norm(R(:,:,2)) + norm(R(:,:,3)))/(3*jmax*kmax);

figure;
ncount = 0;
cpu0 = cputime;

%  only need to compute these once, they are used over and over again
%  and don't change with time
if meth == 0 
  [a_x,b_x,c_x,d_x,e_x] = fill_x(JM,KM,A,dx,dt,jmax,kmax,epse);
  [a_y,b_y,c_y,d_y,e_y] = fill_y(JM,KM,B,dy,dt,jmax,kmax,epse);
end

for nstep = 1:nmax_tot

  ncount = ncount + 1;

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

  resid(nstep) = ...
      (norm(R(:,:,1)) + norm(R(:,:,2)) + norm(R(:,:,3)))/(3*jmax*kmax);
  resid(nstep) = resid(nstep)/resid_0;
  fprintf(' nstep = %g  Resid = %g  \n',nstep,resid(nstep));
  cpu(nstep) = cputime-cpu0;

  if meth == 0
      %  Delta Form Implicit Approximate Factorization
  
      % [I + dt A Del_cen_x ]^(-1) R
      % First Step of Implicit Factored Scheme : x block tridiagonal step
      R = dt*R;
      R = btrix_3x3(2,jmax-1,KM,b_x,c_x,d_x,R);
      
      % [I + dt B Del_cen_y ]^(-1) R
      % Second Step of Implicit Factored Scheme : y block tridiagonal step
      DelQ = btriy_3x3(2,kmax-1,JM,b_y,c_y,d_y,R);
      % Q^(n+1) = Q^n + Delta Q^n
      Q(JM,KM,:) = Q(JM,KM,:) + DelQ(JM,KM,:);
  end
  if meth == 1
      Q(JM,KM,:) = Q(JM,KM,:) + dt*R(JM,KM,:);
  end
  % BC
  Q = bc(Q,JM,KM,jmax,kmax,M_inf);

 
  P = Q(:,2,1);
  cp = (P - 1.0)/(0.5*M_inf^2);
  subplot(3,1,1);
  plot(x,-cp);
  xlabel('x');
  ylabel('-cp');
  subplot(3,1,2);
  contour(X,Y,Q(:,:,1)');
  xlabel('rho');
  subplot(3,1,3);
  semilogy([1:nstep],resid);
  xlabel('n');
  ylabel('resid');
    
  pause(0.1);
 
end

% Fix up corners for plot
Q(1,1,:) = Q(2,1,:); Q(jmax,1,:) = Q(jmax-1,1,:); 
Q(1,kmax,:) = Q(2,kmax,:); Q(jmax,kmax,:) = Q(jmax-1,kmax,:); 
figure;

subplot(2,2,1);
surf(X,Y,Q(:,:,1)');
xlabel('x');
ylabel('y');
zlabel('rho');
view(-25,50);

subplot(2,2,2);
semilogy([1:nmax],resid);
xlabel('n');
ylabel('resid');
subplot(2,2,4);
semilogy(cpu,resid);
xlabel('cpu');
ylabel('resid');

P = Q(:,2,1);
cp = (P - 1.0)/(0.5*M_inf^2);
subplot(2,2,3);
plot(x,-cp);
xlabel('x');
ylabel('-cp');