%
%  w01pplane.m -- 12 jan 04 -- djm
%
%  - gravity'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;
fint = up.^2 - Gm./u;

%  contourplot

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

title('\bf phase plane (universal scaling)')
xlabel('\bf s(\tau)-axis')
ylabel('\bf s''(\tau)-axis')