%
%  code01.m -- 11 jan 07 -- djm
%
%  - pendulum phase plane
%  - to run, type "code01" at the matlab prompt

%  coordinates
N = 100;
the = -2*pi:pi/N:2*pi;
phi =   -pi:pi/N:  pi;

[Gthe Gphi] = meshgrid(the,phi);

%  first integral
Fint = (1/2)*Gphi.^2 - cos(Gthe);

%  contourplot
figure(1);  clf;  hold on
contour(the,phi,Fint,[-1:1/2:6],'k')
contour(the,phi,Fint,[1 1],'r')
contour(the,phi,Fint,-1 + [1 1]*0.1,'b')
axis equal;  axis([-2*pi 2*pi -pi pi]);

title('\bf nonlinear pendulum phase plane','fontsize',12)
xlabel('\bf \theta (t)-axis','fontsize',12)
ylabel('\bf \theta'' (t)-axis','fontsize',12)