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

c----------------------------------------------------------------------
c     Implements differential equations.
c     Called by ODEPACK routine LSODA.
c----------------------------------------------------------------------
      subroutine fcn(neq,t,y,yprime)
         implicit    none

c----------------------------------------------------------------------
c        Problem parameters (G, M) passed in via common block defined
c        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(3)**2)**1.5d0

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

         return
      end 

c----------------------------------------------------------------------
c     Implements Jacobian  (optional) ...
c----------------------------------------------------------------------
      subroutine jac
         implicit    none

         include    'fcn.inc'

         return
      end