%
%  w09herm.m -- (djm:  03 mar 2004)
%    - hermite-gaussian solution
%

clear

%  coordinates
omega = 1;  n = 24;
asq = (2*n+1)/omega;  a = sqrt(asq);
L = 1.5*a;

eta = -L:0.01:L;  Npts = size(eta,2);
x = eta/sqrt(omega);

%  gaussian factor & hermite
gauss = exp(-0.5*eta.^2);

hvects = zeros(n+1,Npts);
for j = 1:Npts
    hvects(:,j) = hermite(n,eta(j))';
end

%  plot eigenfunctions
figure(1);  clf;
subplot(2,1,1);  hold on

for k = n+1:n+1
    y = hvects(k,:).*gauss;  y = y/y((Npts-1)/2);
    plot(x,y,'b')
end
axis([-L L -2 2])
grid

title(['\bf hermite gaussian'])
ylabel(['\bf n = ' num2str(k-1) ' ef'])

subplot(2,1,2);  hold on

for k = 5:5
    y = hvects(k,:).*gauss;  y = y/y((Npts-1)/2);
    plot(x,y,'b')
end
axis([-L L -2 2])
grid

ylabel(['\bf n = ' num2str(k-1) ' ef'])
xlabel('\bf x-axis')

% - - - - - - 
%
%  wkb solution
x0 = -a:a/1001:a;  x0 = x0(2:end-1);
Qx =  0.5*x0.*sqrt(asq - x0.^2) + 0.5*asq*asin(x0/a);
amp = (asq./(asq-x0.^2)).^(1/4);

wkb1 = cos(omega*Qx);
wkb2 = amp.*wkb1;

subplot(2,1,1)
plot(x0,wkb2,'r--')
plot([ a  a],[-2 2],'k--')
plot([-a -a],[-2 2],'k--')

%  turning point connection
ep = (2*a*omega^2)^(-1/3);

y = -10:0.1:10;

%  plot turning points
x1 = a + ep*y;

turn1 = real(airy(0,y));
turn2 = sqrt(pi)*((a/2/ep)^(1/4)) * turn1;

plot(x1,turn2,'m--')
axis([-L L -max(turn2)*1.2 max(turn2)*1.2])