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===========================================================
      subroutine fcn(neq,t,y,yprime)
         implicit    none

         include    'fcn.inc'

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

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

c-----------------------------------------------------------
c        Define some auxiliary quantities to make 
c        computation of RHSs more transparent.
c-----------------------------------------------------------
         xc = y(1)
         yc = y(3)
         th = y(5)
         om = y(6)

c-----------------------------------------------------------
c        Compute positions of two components of the 
c        dumbbell.
c-----------------------------------------------------------
         x1 = xc + d1 * cos(th)
         y1 = yc + d1 * sin(th)
         x2 = xc - d2 * cos(th)
         y2 = yc - d2 * sin(th)

c-----------------------------------------------------------
c        Compute some additional quantities that appear in
c        the equations of motion.
c-----------------------------------------------------------
         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

c-----------------------------------------------------------
c        Compute the RHS's of the ODEs
c-----------------------------------------------------------
         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)

c-----------------------------------------------------------
c        Compute the various energy quantities.
c-----------------------------------------------------------
         ketrans = 0.5d0 * (m1 + m2) * 
     &                    (y(2)**2 + y(4)**2)
         kerot   = 0.5d0 * (m1 * m2) / (m1 + m2) * 
     &                    (d * y(6))**2
         pegrav  = - G * MM * 
     &                    (m1 / sqrt(x1**2 + y1**2) +
     &                     m2 / sqrt(x2**2 + y2**2))
         etot    = ketrans + kerot + pegrav

         return
      end 

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

         include    'fcn.inc'

         return
      end
