%
%  pend.m -- 11 jan 07 -- djm
%
%  - ODE solver for pendulum
%
%- - - - - - - - - - - - - - - - - - - - - - - - 

clear;  code01

%  parameters (damping)
global eta No
eta = 0.0;  No = 20;

%  ODE, timing & accuracy
Xode = 'XpendN';  T = 2*pi;  cycles = 1;  tol = 1e-3;

%  circle of IVs
ang = (0:No-1)*(2*pi/No);  rad = 0.2;
the0 = pi/2 + rad*cos(ang);
phi0 =    0 + rad*sin(ang);

figure(1)
plot(the0,phi0,'r.')

IV = [the0 phi0];

%  ODE solver options
options = odeset('reltol',tol*1e-3,'abstol',tol*1e-6, ...
		'maxstep',T/5,'stats','off','events','off');

%  parallel ODE solve & plot
[t y] = ode45(Xode,[0 T],IV,options);

%  evolution of circle points
theF = y(end,1:No);  phiF = y(end,No+1:end);
plot(theF,phiF,'b.');

%  area calculation
area0 = pi*rad^2

the1  = [theF theF(1)];  phi1 = [phiF phiF(1)];
areaF = sum(the1.*[phi1(2:end) phi1(1)] - phi1.*[the1(2:end) the1(1)])/2