clear all;

% Set some parameters of problem

% In wall:

%dt = 6.7e-10;
%cfl   = dt*0.0e0;     % This is u*dt/dx
%sigma = dt*2.9e4;     % This is k*dt/dx^2
%AR    = 160.0;        % This is dx/dy

% Near the wall:

%dt = 1.77e-7;
%cfl   = dt*5.0e3;  % This is u*dt/dx
%sigma = dt*1.1e2;  % This is k*dt/dx^2
%AR    = 160.0;     % This is dx/dy

% Farfield:

dt = 2e-5;
cfl   = dt*3.3e3;  % This is u*dt/dx
sigma = dt*1.1e2;  % This is k*dt/dx^2
AR    = 1.5;       % This is dx/dy

i = sqrt(-1);

% Plot stability boundary for Forward Euler time integration
theta = linspace(0,2*pi,101);
g = exp(i*theta);
z = g - 1;
plot(z); 
hold on;
xlabel('Real (lambda*dt)');
ylabel('Imag (lambda*dt)');
stitle = sprintf('sigma = %4.2f, cfl = %4.2f, AR = %4.2f',sigma, ...
		 cfl, AR);
title(stitle);
grid on;
axis('equal');

% Plot eigenvalues of 2-d condif with central difference
% approximations

Nxy = 21;
betax = linspace(-pi,pi,Nxy);
betay = linspace(-pi,pi,Nxy);

% Calculate some useful quantities


for ii = 1:Nxy,
  for jj = 1:Nxy,

    bx = betax(ii);
    by = betay(jj);
    lamdt(ii,jj) = 2*sigma*(cos(bx)-1 + AR^2*(cos(by)-1)) ...
	         - cfl*(1.5 - 2*exp(i*bx) + 0.5*exp(2*i*bx));
    
  end
end

plot(real(lamdt),imag(lamdt),'*');