%
%  w09guide.m -- (djm:  04 mar 2004)
%    - phasefronts for parabolic cylinder waveguide
%

%  parameters
omega = 1;  n = 12;  ll = 4.0;
asq = (2*n+1)/omega;  a = sqrt(asq);

%  coordinates
x = -1.0*a:a/100:1.0*a;
y =  0.0:pi*ll/100:2*pi*ll;

[xg,yg] = meshgrid(x,y);

%  eikonal solution
Qx = 0.5*( xg.*sqrt(a^2 - xg.^2) + (a^2)*asin(xg/a) );
phase1 =  Qx + ll*yg;
phase2 = -Qx + ll*yg;

%  plot
figure(1);  clf

%  critical contour & non-propagating traj
ph0 = 0.0;  delph = ll^2 + (a^2)/2;
contour(x,y,phase1,[0:delph*pi/8:delph*pi/2],'b');  
hold on
ph1 = (ll^2)*pi/2 - (a^2)*pi/4;
contour(x,y,phase2,ph1 + [0:delph*pi/8:delph*pi],'r');  
ph2 = (ll^2)*3*pi/2 - (a^2)*pi/4;
contour(x,y,phase1,ph2 + [0:delph*pi/8:delph*pi/2],'g');  
title('\bf phasefronts')
xlabel('\bf x-axis')
ylabel('\bf y-axis')

%  ray trajectory
sigma = 0:pi/100:2*pi;
sigma1 = -pi/2:pi/100:pi/2;

yr = ll*sigma;
yr1 = ll*(sigma1+pi/2);

plot(a*sin(sigma ),yr ,'k')
plot(a*sin(sigma1),yr1+(a^2/2/ll)*pi/2,'k:')
axis equal; axis image