%  Code to do O Delta E stability analysis for Euler Explicit
%  applied to the AA214 class Euler equation project. 
%  Includes: +- flux vector splitting 
%            1st and 2nd order one-sided differences
%
%  Modified Wave number analysis in 2 Dimensions
%  Full System form:
%  Sigma(3,3) = 
%                    +  *      -  *      +  *      -  *
%    eye(3,3) + h (-A ik_bx - A ik_fx - B ik_by - B ik_fy)
% where  
%
%    *
%  ik_bx  is the modified wave number of the 1st or 2nd order 
%         backward difference in x 
%                  *
%         e.g.  i k_bx = [(1-cos(k_x dx)) + i sin(k_x dx)]/dx
%               with k_x the wave number in x
%
%    *      *          *
%  ik_fx, ik_by, and ik_fy  are defined similarly

clear;
M_inf = 0.8;
M_inf = input('M_inf ');

%  Grid Generation
xmin = -1.0;
xmax = 3.0;
ymin = 0.0;
ymax = 2.0;
nx = 10;
nx   = input('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);

% 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;

%  Eigensystem
[XA,LA] = eig(A);
[XB,LB] = eig(B);

%  +/- Flux Jacobian Matrices
Aplus  = 0.5*XA*(LA + abs(LA))*inv(XA);
Bplus  = 0.5*XB*(LB + abs(LB))*inv(XB);
Aminus = 0.5*XA*(LA - abs(LA))*inv(XA);
Bminus = 0.5*XB*(LB - abs(LB))*inv(XB);

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

% k_x dx,   k_y dy  waves numbers
k_x = ([1:jmax])*pi/(jmax+1);
k_y = ([1:kmax])*pi/(kmax+1);

cos_k_x = cos(k_x);
sin_k_x = sin(k_x);
cos_k_y = cos(k_y);
sin_k_y = sin(k_y);
cos_k_2x = cos(2*k_x);
sin_k_2x = sin(2*k_x);
cos_k_2y = cos(2*k_y);
sin_k_2y = sin(2*k_y);

% Modified Wave numbers for each differencing operator
order = 1;
order = menu('Differencing','1st O upwind','2nd O upwind');
mod_kx_c = i*sin_k_x/dx;
mod_ky_c = i*sin_k_y/dy;
if order == 1
  mod_kx_b = ((1-cos_k_x) + i*sin_k_x)/dx;
  mod_ky_b = ((1-cos_k_y) + i*sin_k_y)/dy;
  mod_kx_f = (-(1-cos_k_x) + i*sin_k_x)/dx;
  mod_ky_f = (-(1-cos_k_y) + i*sin_k_y)/dy;
elseif order == 2
  mod_kx_b = ((3-4*cos_k_x+cos_k_2x) + i*(4*sin_k_x-sin_k_2x))/(2*dx);
  mod_ky_b = ((3-4*cos_k_y+cos_k_2y) + i*(4*sin_k_y-sin_k_2y))/(2*dy);
  mod_kx_f = (-(3-4*cos_k_x+cos_k_2x) + i*(4*sin_k_x-sin_k_2x))/(2*dx);
  mod_ky_f = (-(3-4*cos_k_y+cos_k_2y) + i*(4*sin_k_y-sin_k_2y))/(2*dy);
end

% Input a range of CFL #'s 
cfl_in_min = input('Enter min cfl ');
cfl_in_max = input('Enter max cfl ');  
num_cfl = 20; num_cfl = input('Enter # of cfl increments ');
d_cfl = (cfl_in_max-cfl_in_min)/(num_cfl-1);
if cfl_in_min == cfl_in_max num_cfl = 1; end

% Use the definition of CFL from the Projects
Eig_x = 1.0 + M_inf; Eig_xp = Eig_x; Eig_xm =  M_inf-1;
Eig_y = 1.0;         Eig_yp = Eig_y; Eig_ym = -1.0;
Eig_t = Eig_x + Eig_y;

% Sigma(3,3) is a 3x3 matrix defining sigma
% Form it and find the max eigenvalues over 
% a spectrum of k_x and k_y

for n = 1:num_cfl
  cfl_n(n) = (cfl_in_min + (n-1)*d_cfl);
  dt = min(dx,dy)*cfl_n(n)/Eig_t;   dt_n(n) = dt;

  for k = 1:kmax
    for j = 1:jmax
      Lam_h = -(Aplus*mod_kx_b(j) + Aminus*mod_kx_f(j) +  ...
	  Bplus*mod_ky_b(k) + Bminus*mod_ky_f(k) );
      Sigma = eye(3,3) + Lam_h*dt;
      % $$$       sigma(j,k,n) = max(abs(eig(Sigma)));
      sigma(j,k,n) = max(eig(Sigma));
    end
  end
  % $$$   sig_max(n) = max(max(abs(sigma(:,:,n))));
  sig_max(n) = max(max(abs(sigma(:,:,n))));
end
format short g;
disp( ['     n          cfl           dt       sigmax ']);
[ [1:num_cfl]' cfl_n' dt_n' sig_max']
figure;
i_one = 1;
for i = 1:num_cfl
  if sig_max(i) > 1
    i_one = i;
    break;
  end
end
if i_one >= 3 & i_one <= num_cfl-3 subplot(2,1,1); end
plot(cfl_n,sig_max);
hold on;
plot(cfl_n,ones(num_cfl,1),'k:');
if i_one >= 3 & i_one <= num_cfl-3
  subplot(2,1,2);
  plot(cfl_n(i_one-2:i_one+2),sig_max(i_one-2:i_one+2));
hold on;
plot(cfl_n(i_one-2:i_one+2),ones(5,1),'k:');
end