function out = pcrtbp(t,y)
% 
%   Usage:
%   y' = pcrtbp(t,y)  or  H = pcrtbp(y)  or  y0 = pcrtbp()      
%
%   planar circular restricted three-body problem. y = [q1 q2 p1 p2].
%
%   See pages 140 and 142-144  of
%   L.Brugnano, F.Iavernaro.
%   Line Integral Methods for Conservative Problems.
%   Chapman and Hall/CRC - 2015
%   Series: Monographs and Research Notes in Mathematics.
%
mt = 0.3;
ms = 1-mt;
if nargin==0      %% initial condition
 out = [6.955415851757456e-01  ...  
        6.456133089565784e-03  ...
        -7.663758431829502e+00 ...
        4.716929963135788e+00];
elseif nargin==2   %% vector field 
   q1     = y(:,1);
   q2     = y(:,2);
   p1     = y(:,3);
   p2     = y(:,4);
   ro1    = sqrt( ( q1 +mt ).^2 +q2.^2 );
   ro2    = sqrt( ( q1 -ms ).^2 +q2.^2 );
   dro1q1 = -( q1 +mt )./ro1.^3;
   dro1q2 = -q2./ro1.^3;
   dro2q1 = -( q1 -ms )./ro2.^3;
   dro2q2 = -q2./ro2.^3;
   dq1    = p1 +q2;
   dq2    = p2 -q1; 
   dp1    =  p2 +ms.*dro1q1 +mt.*dro2q1;
   dp2    = -p1 +ms.*dro1q2 +mt.*dro2q2;
   out    = [dq1 dq2 dp1 dp2];  
elseif nargin==1  %% Hamiltonian
   q1    = t(:,1);
   q2    = t(:,2);
   p1    = t(:,3);
   p2    = t(:,4);
   ro1   = sqrt( ( q1 +mt ).^2 +q2.^2 );
   ro2   = sqrt( ( q1 -ms ).^2 +q2.^2 );
   omega = .5*( q1.^2 +q2.^2 ) +ms./ro1 +mt./ro2;
   out   = [( 0.5.*( p1.^2 +p2.^2 ) +p1.*q2 -p2.*q1 -ms./ro1 -mt./ro2 )];
end
return