%
%  lect22a.m -- stability analysis
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

%  stab = c dt/dx

stab = 0.75

%  basic parameters

N = 128;  L = pi;   	dx = 2*L/N;
c = 1.0;  		dt = stab*(dx/c);

%  wavenumber vector

n1 = -N/2 :1: N/2;

%  amplification factor

zu = 1 - (c*dt/dx)*(1 - exp(-i*(2*pi/N)*n1));
zd = 1 + (c*dt/dx)*(1 - exp( i*(2*pi/N)*n1));
zc = 1 - (c*dt/dx/2)*(exp( i*(2*pi/N)*n1) - exp(-i*(2*pi/N)*n1));

%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
%  plot here

clf

subplot(2,1,1)
plot(n1/N,abs(zu)-1,'b')
hold on
plot(n1/N,abs(zd)-1,'b--')
plot(n1/N,abs(zc)-1,'g')

subplot(2,1,2)
plot(n1/N,unwrap(angle(zu))+stab*2*pi/N*n1 + 2*pi,'b')
hold on
plot(n1/N,unwrap(angle(zd))+stab*2*pi/N*n1,'b--')
plot(n1/N,unwrap(angle(zc))+stab*2*pi/N*n1,'g')

subplot(2,1,1)
title('\bf upwind, downwind & centered')
xlabel('\bf wavenumber (n/N)')
ylabel('\bf amplitude error (|Z|-1)')
axis([0 0.5 -1 1])

subplot(2,1,2)
title('\bf phase error')
xlabel('\bf wavenumber (n/N)')
ylabel('\bf phase error (\phi - \phi_{true})')
axis([0 0.5 -pi/2 pi])