%
%  pperiod.m -- djm -- 11 jan 2007
%
%  - period
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - -

%  parameters
global eta
eta = 0.0;

%  ODE, timing & accuracy (events on)
Xode = 'Xpend';  T = 8*pi;  tol = 1e-3;
options = odeset('reltol',tol*1e-3,'abstol',tol*1e-6, ...
		'maxstep',T/5,'stats','off','events','on');

%  loop of ODE solves    
figure(3);  clf;  hold on
period = [];  dIV = pi/25;

for IV = dIV:dIV:pi-dIV
    [t y] = ode45(Xode,[0 T],[IV 0],options);
    
    plot(IV,0,'rx');  plot(y(end,1),y(end,2),'ro')
    plot(y(:,1),y(:,2))
    axis equal;  axis([-pi pi -2 2])
    period = [period ; IV t(end)];
    drawnow
    
    disp(['amp = ' num2str(IV) ' ; period = ' num2str(t(end))])
end
title('\bf periodic orbits','fontsize',12)
xlabel('\bf \theta (t)','fontsize',12)
ylabel('\bf \theta'' (t)','fontsize',12)

%  final plot
figure(4);  clf;  hold on
plot(period(:,1),period(:,2),'.')
plot(0,2*pi,'r.')
plot([pi pi],[0 max(period(:,2))],'k--')
title('\bf period vs initial amplitude','fontsize',12)
xlabel('\bf initial amplitude','fontsize',12)
ylabel('\bf period','fontsize',12)