%
%  code08.m -- steepest descent paths (djm, 02 oct 01)
%
%  overlay for "code07.m"

zz = 9;

%  parameters
Ns = 100;  Ls = 1.5;  ds = Ls/Ns;  s = -Ls:ds:Ls;

%  steepest descent quantities
 
k_sd = sqrt(-zz);  
iphi_sd = i*k_sd * (2*zz/3);  C = imag(iphi_sd);

v1 = sqrt(  zz)*cosh(s);
u1 = sqrt(3*zz)*sinh(s);

v2 =    sqrt(zz)*sinh(s);
u2 = -2*sqrt(zz)*cosh(s).*cos( (sign(s).*acos(sech(s).^3)+2*pi)/3);

v3 =    sqrt(zz)*sinh(-s);
u3 =  2*sqrt(zz)*cosh(-s).*cos( (sign(-s).*acos(sech(-s).^3)+2*pi)/3);

figure(1);  plot(u1,v1,'r--');
figure(2);  plot(u2,v2,'r--');  plot(u3,v3,'g--');


figure(3);  clf

subplot(2,1,1)
k1 = u1 + i*v1;  phi = i*(zz*k1 + (k1.^3)/3);
integrand = exp(real(phi)-max(real(phi)));
phase = imag(phi) - 0;

plot(s,phase); hold on;  plot(s,integrand,'r');
axis([-Ls Ls -0.1 1.1])

title(['\bf Airy steepest descent (normalized) integrand & (constant) phase:  z = ' num2str(zz)])
xlabel('\bf s-axis');  ylabel('\bf integrand (red) & phase (blue)');  

subplot(2,1,2)
k2 = u2 + i*v2;  phi = i*(-zz*k2 + (k2.^3)/3);
integrand = exp(real(phi)-max(real(phi)));
phase = imag(phi) + (2/3)*(zz)^(3/2);

plot(s+Ls,phase);  hold on;  plot(s+Ls,integrand,'r');

k3 = u3 + i*v3;  phi = i*(-zz*k3 + (k3.^3)/3);
integrand = exp(real(phi)-max(real(phi)));
phase = imag(phi) - (2/3)*(zz)^(3/2);

plot(s-Ls,phase);  hold on;  plot(s-Ls,integrand,'g');
axis([-2*Ls 2*Ls -0.1 1.1])

title(['\bf Airy steepest descent (normalized) integrand & (constant) phase:  z = -' num2str(zz)])
xlabel('\bf s-axis');  ylabel('\bf integrand (green/red) & phase (blue)');