% Euler_template V1.0
%  Final Project 
%
% THP 04/2K+1
% 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    |

%         M_inf :  Free Stream Mach number
%         nx    :  Number of points on airfoil (x=0,1)
%
%  Inputs:
%         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
%
clear;
M_inf = 0.8;  % remove comment % $$$ 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(' X of Leading edge point = %g  X of trailing edge point = %g \n',x(jle),x(jte));

% Surface Definition
tau = 0.1; % tau = input('Enter 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
cfl  = input('Enter cfl = ');
nmax = input('Enter nmax = ');

dt = min(dx,dy)*cfl/(2.0+M_inf);
fprintf(' AT CFL = %g  dt = %g \n',cfl,dt);

for nstep = 1:nmax

  %  rhs Central differencing
  Qx_c(JM,KM,:) = 0.5*(Q(JM+1,KM,:) - Q(JM-1,KM,:))/dx;
  Qy_c(JM,KM,:) = 0.5*(Q(JM,KM+1,:) - Q(JM,KM-1,:))/dy;
  for m = 1:3
    R(JM,KM,m) = 0.0;
    for n = 1:3
      R(JM,KM,m) = R(JM,KM,m) + ...
	  A(m,n)*Qx_c(JM,KM,n) + ...
	  B(m,n)*Qy_c(JM,KM,n) ;
    end
  end

  % $$$   You can add these next 4 lines and run the code with cfl = 0.1
  % $$$   for 1000 steps to see what a reasonable solution would be.
  % $$$   I added this in to give yuou at a sample working code.  
  % $$$   Remove % line#
  % line1   epse = 0.1;
  % line2   Qx_d(JM,KM,:) = (Q(JM+1,KM,:) -2*Q(JM,KM,:) + Q(JM-1,KM,:))/dx;
  % line3   Qy_d(JM,KM,:) = (Q(JM,KM+1,:) -2*Q(JM,KM,:) + Q(JM,KM-1,:))/dy;
  % line4   R(JM,KM,:) = R(JM,KM,:)-epse*(Qx_d(JM,KM,:) + Qy_d(JM,KM,:));
  
  resid(nstep) = ...
      (norm(R(:,:,1)) + norm(R(:,:,2)) + norm(R(:,:,3)))/(3*jmax*kmax);
  if nstep == 1 resid_norm = resid(1);  end
  resid(nstep) = resid(nstep)/resid_norm;
  fprintf(' nstep = %g  Resid = %g  \n',nstep,resid(nstep));

  Q(JM,KM,:) = Q(JM,KM,:) - dt*R(JM,KM,:); % Euler Explicit Step

  % BC
  % At j = 1 if M-inf < 1 Extrap rho
  if M_inf < 1.0 
    Q(1,   KM,1) = Q(2,      KM,1); 
  end
  % At jmax DQ/Dx = 0 using backward diff. 
  Q(jmax,KM,1) = Q(jmax-1, KM,1);
  Q(jmax,KM,2) = Q(jmax-1, KM,2);
  Q(jmax,KM,3) = Q(jmax-1, KM,3);

  % At k=1 D rho/Dy = 0, Du/Dy = 0, fix v
  RHO = Q(JM,1,1);
  Q(JM,1,1) = Q(JM,2,1);
  Q(JM,1,2) = Q(JM,2,2)./Q(JM,2,1).*Q(JM,1,1);
  % fixing v, not rho v so this little rescaling is necessary
  Q(JM,1,3) = Q (JM,1,3)./RHO.*Q(JM,1,1);  

  % At kmax fix all do nothing 

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);
U = Q(:,:,2)./Q(:,:,1);
V = Q(:,:,3)./Q(:,:,1);
surf(X,Y,U');
xlabel('x');
ylabel('y');
zlabel('u');
view(-25,50);


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

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

title(['Euler Explicit','  CFL = ',num2str(cfl),...
      '  M_{\infty} = ',num2str(M_inf)]);