function swing(eps,th,u,Fmx)

% Phase portrait and Poincare map for pumped swing problem
%
% Inputs; 
%  eps    'Pumping' parameter scaling swing length oscillation from pumping
%  th, u  vectors of values defining all coordinate pairs v=(th,u)
%         at which to plot the vector F = (dth/dt,du/dt). The 
%         positive part of the plotted range of y is also used to make
%         a Poincare map
%  Fmx    sets a plotting threshold for |F| such that if |F|>Fmax,
%         we plot F*(Fmx/|F|) rather than F. This allows short
%         vectors F to show up on plot.
%  Example call:
%   eps = 0.1;
%   th = -1:0.2:1;
%   u = -1:0.2:1;
%   Fmx = 0;
%   swing(eps,th,u,Fmx)

  s = 0.4;
  
  figure(1)

  PhasePortrait(@dvdt,th,u,s,Fmx)
  axis equal
  xlim([min(th) max(th)])
  ylim([min(u) max(u)])
  xlabel('\theta (radians)')
  ylabel('u = d\theta/dt')
  title('Pumped swing')

  % Plot fixed points in black
  
  hold on
  plot(0,0,'ko')
  
  % Plot some orbits in red
  
  dt = 0.1*pi;

  thmax = max(th);
  th0 = 0.3*(1:3);
  u0 = zeros(1,3);
  tmax = 2*pi*ones(1,3);
  for j = 1:length(th0)
    tspan = 0:dt:tmax(j);
    v0 = [th0(j); u0(j)];
    [t, v] = ode15s(@dvdt,tspan,v0);
    plot(v(:,1),v(:,2),'r.')
  end
  hold off  

  figure(2)
  
  % Poincare map using events option of ode15s
  
  options = odeset('Events',@events);

  n0 = 50;
  th0 = thmax*(1:n0)/n0;
  u0 = zeros(1,n0);
  tmax = 3*pi*ones(1,n0);
  te = zeros(1,n0);
  the = zeros(1,n0);  
  for j = 1:n0
    tspan = [0 tmax(j)];
    v0 = [th0(j); u0(j)];
    [t,v,tev,ve,ie] = ode15s(@dvdt,tspan,v0,options);
    if(length(ve)>0)
      the(j) = ve(end,1);
    else
      the(j) = NaN;
    end
  end
  plot(th0,the,'.',th0,th0,'--')
  xlim([min(th0) max(th0)])
  xlabel('\theta_n')
  ylabel('\theta_{n+1}')
  title('Poincare map')
%-----------------------------------------------------------------  
function F = dvdt(t,v)

%  Calculates time derivative of v = [theta; u=theta'] based on the 
%  ODE  
%    r*theta'' + 2*r'*theta' + theta = 0
%  where r = 1 - eps*theta*theta'/(theta^2 + theta'^2).  This turns into
%  the system  
%    theta' = u
%  and  
%    r*du/dt = -2*r'*u -theta 
%  which must be rewritten using the chain rule 
%    r' = (dr/du)*u' + (dr/dtheta)*theta'
%  and theta' = u as:
%    (r + 2*(dr/du)*u)*du/dt = -2*(dr/dtheta)*u^2 - theta
  
  th = v(1,:);
  u = v(2,:);
  r = 1 - eps*th.*u./(th.^2 + u.^2);
  drdu = eps*th.*(u.^2 - th.^2)./(th.^2 + u.^2).^2;
  drdth = eps*u.*(th.^2 - u.^2)./(th.^2 + u.^2).^2;
  dudt = -(2*(drdth).*u.^2 + th)./(r + 2*(drdu).*u);
  
  F = [u; dudt];
  end
%-------------------------------------------------------------
function [value,isterminal,direction] = events(t,v)

  value = v(2);  % Detect if theta'=0
  isterminal = 1; % Stop integration 
  direction = -1; % ...but only if theta' is decreasing.
  end
  
end