%
%  rays.m -- (djm:  24 feb 2004)
%    - phase plane for sinusoidal index medium
%

%  phase plane coordinates
xp = -pi:pi/100:pi;
x = -pi:pi/100:pi;

[xg,xpg] = meshgrid(x,xp);

%  first integral
a = 0.5;
yg = (1 + a*cos(xg)) - xpg.^2;

%  plot
figure(1);  clf

%  critical contour & non-propagating traj
contour(x,xp,yg,[0 0],'r');  
hold on
contour(x,xp,yg,-[1:1:10]*2,'b:');  

%  propagating traj
contour(x,xp,yg, [1:1:3]*0.125,'b');  
contour(x,xp,yg, [4 4]*0.125,'r--');  
contour(x,xp,yg, [5:1:12]*0.125,'k');  

title('\bf optical ray phase plane')
xlabel('\bf x-axis')
ylabel('\bf x''-axis')