c===========================================================
c     Solves EOM for orbiting dumbbell (rigid body composed 
c     of 2 point masses m1 and m2, separation d)
c
c     See class notes for equations of motion.
c
c     Canonicalization:
c
c         y(1) = xc
c         y(2) = d(xc)/dt
c         y(3) = yc 
c         y(4) = d(yx)/dt
c         y(5) = th
c         y(6) = d(th)/dt
c===========================================================
c     neq = 6
c-----------------------------------------------------------

c-----------------------------------------------------------
c     Implements differential equations.
c-----------------------------------------------------------
      subroutine fcn(neq,t,y,yprime)
         implicit    none

         include    'fcn.inc'

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

         real*8      xc,    yc,    th,    
     &               c1,    c2,
     &               r1m3,  r2m3

         xc = y(1)
         yc = y(3)
         th = y(5)

         x1 = xc + d1 * cos(th)
         y1 = yc + d1 * sin(th)
         x2 = xc - d2 * cos(th)
         y2 = yc - d2 * sin(th)

         r1m3 = 1.0d0 / (x1**2 + y1**2) ** 1.5d0
         r2m3 = 1.0d0 / (x2**2 + y2**2) ** 1.5d0

         c1 = -G * MM 
         c2 =  G * MM / d

         yprime(1) = y(2)
         yprime(2) = c1 * ((1.0d0 - mu) * x1 * r1m3 +
     &                              mu  * x2 * r2m3)
         yprime(3) = y(4)
         yprime(4) = c1 * ((1.0d0 - mu) * y1 * r1m3 +
     &                              mu  * y2 * r2m3)
         yprime(5) = y(6)
         yprime(6) = c2 * (r1m3 - r2m3) *
     &                    (sin(th) * xc - cos(th) * yc)

         xc = y(1)
         yc = y(3)
         th = y(5)

         x1 = xc + d1 * cos(th)
         y1 = yc + d1 * sin(th)
         x2 = xc - d2 * cos(th)
         y2 = yc - d2 * sin(th)

         return
      end 

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

         include    'fcn.inc'

         return
      end
