c===========================================================
c     Implements (planar) equations of motion for restricted 
c     2-body gravitational problem.  Central mass, M, is
c     fixed at (0,0).  Mass of other object with coordinates 
c     (x_c,y_c) is gravitationally negligible. 
c     ' denotes differentiation with respect to t.
c     
c     y(1) := x_c
c     y(2) := y_c
c     y(3) := x_c'
c     y(4) := y_c'
c===========================================================
      subroutine fcn(neq,t,y,yprime)
         implicit    none

c-----------------------------------------------------------
c        Problem parameters (G, M) passed in via common 
c        block defined in 'fcn.inc'
c-----------------------------------------------------------
         include    'fcn.inc'

         integer     neq
         real*8      t,     y(neq),    yprime(neq)

         real*8      c1

         c1 = -G * M / (y(1)**2 + y(2)**2)**1.5d0

         yprime(1) = y(3)
         yprime(2) = y(4)
         yprime(3) = c1 * y(1)
         yprime(4) = c1 * y(2)

         return
      end 

c===========================================================
c     Computes mechanical energy (etot) and angular momentum 
c     about the origin (location of the gravitating mass)
c     from the dynamical variables.  "Specific" quantities
c     (i.e. normalized by the mass of the dynamical test 
c     particle) are computed.
c===========================================================
      subroutine calc_ej(y,etot,jtot)
         implicit      none
         real*8        y(4),    etot,    jtot

         include      'fcn.inc'

         etot = 0.5d0 * (y(3)**2 + y(4)**2) - 
     &          G * M / sqrt(y(1)**2 + y(2)**2)
         jtot = y(1) * y(4) - y(2) * y(3)

         return
      end

c===========================================================
c     Dummy Jacobian routine.
c===========================================================
      subroutine jac
         implicit    none

         include    'fcn.inc'

         return
      end
