%
%  w02pplane1.m -- 12 jan 04 -- djm
%
%  - orbital s-s' phase plane

%  coordinates

s  =  0.1:0.02:2.0;
sp = -2.0:0.02:2.0;

[u up] = meshgrid(s,sp);

%  first integral

Gm = 1.0;
L  = 0.25;
fint = up.^2 - Gm./u + L./(u.^2);

%  contourplot

figure(2);  clf;  hold on
contour(s,sp,fint,[-3:0.25:0  ],'k--')
contour(s,sp,fint,[ 0:0.25:1.4],'k')
contour(s,sp,fint,[0 0],'r')
axis([0 2 -1.5 1.5])

ptitle = ['\bf sun-earth distance phase plane (L = ' num2str(L) ')'];
title(ptitle)
xlabel('\bf s-axis')
ylabel('\bf s''-axis')