%This code will attempt to find the stability bounds on the Euler Explicit
%and RK4 methods

%Our difference Operator has the form
%
%dQ
%--=-A+dbx(Q)-A+dfx(Q)-B+dby(Q)-B-dfy(Q)
%dt
%
%Fourier Analsyis
%dQ
%--=(-(ik*A+)-(ik*A-)-(ik*B+)-(ik*B-))
%dt
clear;
tmeth=input('Time March\n 0 for Euler Explicit; 1 for RK4:');
spmeth=input('Space Method\n 0 for Second order upwind; 1 for compact:');

M_inf=.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,:)=[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);

%Figure Out the maximum sigma

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

kxmax=length(k_x);   kymax=length(k_y);

maxsigmas=zeros(kxmax,kymax);
maxabssigmas=zeros(kxmax,kymax);
cfl=(0.1:.2:3);
maxsignums=zeros(length(cfl),3);

Eig_t=2.0+M_inf;
% 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;

for cflind=1:length(cfl)
    
    dt=min(dx,dy)*cfl(cflind)/Eig_t;

    for kxind=1:kxmax
        for kyind=1:kymax
        
            if spmeth==0
              %Get values for the ik*'s
              ikstararr=ikstar(dx,dy,k_x(kxind),k_y(kyind),spmeth);
        
              %Form the Spatial Operator using thes ik*'s
              Lambda_matrix=-(ikstararr(1)*Am+ikstararr(2)*Ap+...
                          ikstararr(3)*Bm+ikstararr(4)*Bp);
            elseif spmeth==1
              %Get values for the ik*'s
              ikstararr=ikstar(dx,dy,k_x(kxind),k_y(kyind),spmeth);
              Lambda_matrix=(ikstararr(1)*A+ikstararr(2)*B);
           end
            
           if tmeth==0
             %Form the sigma matrix          
             sigma=eye(3)+dt*Lambda_matrix;
           elseif tmeth==1
             sigma=eye(3)+dt*Lambda_matrix+.5*dt^2*Lambda_matrix^2+...
                 (1/6)*dt^3*Lambda_matrix^3+(1/24)*dt^4*Lambda_matrix^4;
        
           end
           [evec,eval]=eig(sigma);
           maxsigmas(kxind+1,kyind+1)=eval(find(max(max(abs(eval)))));
           maxabssigmas(kxind+1,kyind+1)=max(max(abs(eval)));
           
        end
    end
    maxsignums(cflind,3)=max(max(maxabssigmas));    
    maxsignums(cflind,2)=cfl(cflind);
    maxsignums(cflind,1)=dt;
    
end

if spmeth==0
    spstring='Second Order Backward/Forward';
elseif spmeth==1
    spstring='Compact Scheme';
end

if tmeth==0
    tstring='Euler Explicit';
elseif tmeth==1
    tstring='RK4'
end

format short g;
disp( ['      cfl            dt          sigmax ']);
[ maxsignums(:,2) maxsignums(:,1) maxsignums(:,3)]
plot(cfl,maxsignums(:,3));

title(sprintf('Stability for %s with %s with nx=10',tstring,spstring));
xlabel('CFL number');
ylabel('Maximum value of sigma');

%get the max and min and their corresponding CFL number
maxcfl=find(maxsignums(:,3)==max(maxsignums(:,3)));
mincfl=find(maxsignums(:,3)==min(maxsignums(:,3)))
maxcfl=maxsignums(maxcfl,2);
mincfl=maxsignums(mincfl,2);

minsig=text(.1,.85,'Min Sig=','units','normalized');
maxsig=text(.1,.9,'Max Sig=','units','normalized');
set(minsig,'string',sprintf('Min Sig=%2.12f at CFL=%2f',min(maxsignums(:,3)),maxcfl));
set(maxsig,'string',sprintf('Max Sig=%2.12f at CFL=%2f',max(maxsignums(:,3)),mincfl));


%plot(real(maxsigmas),imag(maxsigmas),'*');