function poincare_pendulum(kappa_in,epsilon_in,drive_periods,pmax,markersize);

% poincare_pendulum(kappa,epsilon,plot_points,markersize);
%
% Calculates and plots a Poincare section for the vertically driven pendulum:
%
%     d^2 phi / dt^2  = - kappa sin(phi) (1 - epsilon cos t)
%
% Required inputs:
% kappa = (omega_0/omega_D)^2 : omega_0 is SHO frequency of pendulum, 
%                               omega_D is driving frequency
% epsilon                     : Strength of driving potential 
%
% Optional inputs:
% drive_periods               : number of points to plot (default 1000)
% marker_size                 : size of points in Poincare section (default 6)
%
% You choose initial conditions by clicking on the figure with the
% left mouse button.  The computer then solves Hamilton's equations
% for the number of driving periods.
%
% The right mouse button will then continuing solving for the same 
% period of time if you want more points on your figure.
%
% However, you can now select another initial condition with the left
% mouse button, and slowly build up your Poincare section. 
%
% To quit, click outside the axes of the figure.
%
% Matthew Davis
% 21st October 2006.

global kappa epsilon
kappa = kappa_in;
epsilon = epsilon_in;

% Number of points to plot if not set by user
if (nargin < 3)
  drive_periods = 1000;
end;
if (nargin < 4)
  pmax = 3;
end;
if (nargin < 5)
  marker_size = 6;
end;

% Parameter range for Poincare section
qmax = pi;
pmax = abs(pmax);


% Set up figure
fig=figure(123);clf
set(fig,'units','normalized','position',[0.1 0.1 .8 .8]);
plot(0,0,'w');
axis([-qmax qmax -pmax pmax]);
xlabel('Angle','fontsize',30)
ylabel('Angular momentum','fontsize',30)
title(['Poincare section for \kappa = ',num2str(kappa),...
       ', \epsilon = ',num2str(epsilon)],'fontsize',30)
hold on

% Parameters for integration
N = 400;
dt = 2*pi/N;
tn = zeros(drive_periods+1,1);
qn = tn;
pn = tn;

% Variables used in control loop
colours = 'bgrcmk';
count = 0;
stop = 0;

while ~stop
  % Get q and p from user input
  [q,p,button] = ginput(1);
  t = 0;
  % If not a left click, use the final point from last time
  if ((button ~= 1) & (count > 0))
    q = qn(end);
    p = pn(end);
    t = tn(end);
  else
    % Otherwise change the colour for the plot points
    count = count + 1;
    if (count > length(colours))
      count = 1;
    end;
  end;

  % If the user has clicked outside the plot region then quit
  if ((button == 1 ) & ((abs(q) > qmax) | (abs(p)>pmax)))
    stop = 1;
  else
  % Else integrate Hamilton's equations
    tn(1) = t;
    qn(1) = q;
    pn(1) = p;
    for jj = 2:(drive_periods+1);
      [tn(jj),qn(jj),pn(jj)] = dyrting_step(dt,tn(jj-1),qn(jj-1),pn(jj-1),N);
    end;
    % Shift the q variable so that it is in the range (-pi,pi)
    qn = mod(qn,2*pi);
    index = find(qn>pi);
    qn(index) = qn(index) - 2*pi;
    % Plot the points
    plot(qn,pn,['.',colours(count)],'Markersize',marker_size);
  end;
end;

% No output is needed
y=[];



function [t,q,p] = dyrting_step(dt,t,q,p,N);
global kappa epsilon

% This is a symplectic integrator for Hamilton's equations for the vertically
% driven pendulum.  
%
% See Sigurd Dyrting PhD thesis, University of Queensland, 1995.

b = 2^(1/3);
d1 = dt/(2-b);
d2 = -b*d1;
s1 = 0.5*d1;
s2 = (1-b)*s1;

for jj = 1:N;
q = q + p*s1;t = t + s1;
p = p - d1*kappa*sin(q)*(1-epsilon*cos(t));
q = q + p*s2;t = t + s2;
p = p - d2*kappa*sin(q)*(1-epsilon*cos(t));
q = q + p*s2;t = t + s2;
p = p - d1*kappa*sin(q)*(1-epsilon*cos(t));
q = q + p*s1;t = t + s1;
end;
