%
%  code01a.m -- 06 jan 04 -- djm
%
%  - gravity's phase plane
%  - to run, type "code01a" at the matlab prompt

%  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.5:4],'k')
contour(s,sp,fint,[0 0],'r')

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