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

%  parameters (damping)
global eta
eta = 0.0;

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

%  too fancy IV input (discretized for repeatability)
disp(' ');  disp(' click to choose initial point');  disp(' ');

figure(1)
in = [];  IV = [2*pi pi];
while (isempty(in))
  in = floor(ginput(1)/(pi/25))*(pi/25); 

  if(isempty(in))
    break
  else
    plot(IV(1),IV(2),'w.')
    disp(['  IV:  ' num2str(in) '  (rtn to confirm)']); 
    IV = in;  in = [];
    plot(IV(1),IV(2),'r.')
  end
end
plot(IV(1),IV(2),'w.');  plot(IV(1),IV(2),'rx');

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

%  initial first integral
First0 = (1/2)*IV(2)^2 - cos(IV(1));
disp(['  ' num2str(0) '  ' num2str(First0)]);
figure(2);  clf;  hold on
plot(0,IV(1),'b.');  plot(0,IV(2),'r.')


%  ODE solve & plot
for C = 1:cycles
  [t y] = ode45(Xode,[(C-1)*T C*T],IV,options);

  %  update IV by last point
  IV = y(end,:);  

  %  monitor first integral
  First = (1/2)*y(end,2)^2 - cos(y(end,1));
  disp(['  ' num2str(t(end)) '  ' num2str(First) '  ' num2str(First-First0)])

  figure(1);
  plot(y(end,1),y(end,2),'b.');

  figure(2);
  plot(t,y(:,1),'b')
  plot(t,y(:,2),'r')
  plot(t(end),y(end,1),'b.');
  plot(t(end),y(end,2),'r.');
  title('\bf pendulum dynamics','fontsize',12)
  xlabel('\bf time','fontsize',12)
  ylabel('\bf \theta (t) = blue ; \theta'' (t) = red','fontsize',12)

  drawnow
end