subroutine rkint(dt,u)

*	Runge-Kutta integrator.
*       -----------------------

*       Author:  Rosemary Mardling (3/98).

      parameter (n=6)
      implicit real*8 (A-H,O-Z)
      real*8 u0(n),ut(n),du(n),u(n)
      real*8 a(4),b(4)

      iq = 0
      a(1)=dt/2.0d0
      a(2)=a(1)
      a(3)=dt
      a(4)=0.0d0

      b(1)=dt/6.0d0
      b(2)=dt/3.0d0
      b(3)=b(2)
      b(4)=b(1)

      do i=1,n
         u0(i)=u(i)
         ut(i)=0.0d0
      enddo

      do j=1,4
         call deriv(u,du,iq)

         do i=1,n
            u(i)=u0(i)+a(j)*du(i)
            ut(i)=ut(i)+b(j)*du(i)
         enddo
      enddo

      if(iq.gt.0)then
         do i=1,n
            ut(i)=0.0d0
         enddo
      endif

      do i=1,n
         u(i)=u0(i)+ut(i)
      enddo

      end


      subroutine deriv(u,du,iq)

*	Differential equations for hierarchical binary.
*       ----------------------------------------------

*       Author:  Rosemary Mardling (3/98).

      implicit real*8 (A-H,M,O-Z)
      common/rksave/  coeff,HOhat(3),e,a,hh,mb
      common/tidal/  cq(2),ct(2),cgr,dedt
      real*8  ehat(3),hhat(3),qhat(3),edot(3),hdot(3),u(6),du(6)
      SAVE ITIME,IGR
      DATA ITIME,IGR /0,0/
	

*       Save initial basic elements.
      e0 = e
      a0 = a
      h0 = hh
 
*       Update eccentricity and angular momentum.
      e=sqrt(u(1)**2+u(2)**2+u(3)**2)
      if(e.ge.1.0)then
         iq = 1
         goto 10
      endif
      hh=sqrt(u(4)**2+u(5)**2+u(6)**2)
      a=hh**2/mb/(1.0-e**2)

*       Define unit vectors.
      do i=1,3
         ehat(i)=u(i)/e
         hhat(i)=u(i+3)/hh
      enddo

*	Calculate unit vector orthogonal to ehat and hhat (Peter's q).

      call cross(hhat,ehat,qhat)

*	Calculate components of Peter's Sij tensor (third body average).

      S11=-coeff*(3*(dot(HOhat,ehat))**2 - 1)
      S12=-3*coeff*dot(HOhat,ehat)*dot(HOhat,qhat)
      S13=-3*coeff*dot(HOhat,ehat)*dot(HOhat,hhat)
      S22=-coeff*(3*(dot(HOhat,qhat))**2 - 1)
      S23=-3*coeff*dot(HOhat,qhat)*dot(HOhat,hhat)

*	Calculate rate of change of evec and hvec.

      do i=1,3
         edot(i)=(e*a*hh/2/mb)*
     &           (-5*S12*ehat(i)+(4*S11-S22)*qhat(i)-S23*hhat(i))

         hdot(i)=(a**2/2)*((1-e**2)*S23*ehat(i)
     &          -(1+4*e**2)*S13*qhat(i)+5*e**2*S12*hhat(i))
      enddo

      EKD = edot(1)*ehat(1)+edot(2)*ehat(2)+edot(3)*ehat(3)
      TAV = e/sqrt(edot(1)**2+edot(2)**2+edot(3)**2)
*
*       Adopt eccentricity functions from Hut theory (A & A 99, 126, 1981).
      ge=30*e+45*e**3+3.75*e**5
      ge=0.5*ge
      he=(1.0+3.75*e**2+1.875*e**4+5.0/64.0*e**6)/(1.0-e**2)**6.5

*       Specify scaling factors for quadrupole and tidal interaction.
      slr=a*(1.0-e**2)
      zq=ge/(slr**5*a*sqrt(a))
      zz=he/a**8
      zt=9.0*e*zz*(ct(1) + ct(2))
      zq=zq*(cq(1) + cq(2))
*       Correct for wrong eccentricity dependence.
      zq=zq*(1.0 - e**2)
*       Note that angular momentum derivative is of same form as de/dt.
*     zh=hh*zz*(ct(1) + ct(2))
*
*       Form scaled expression for GR precession.
      zg = cgr/(slr*a*sqrt(a))

*     zq = 0.0
*     zg = 0.0
*     zt = 0.0
*     zh = 0.0
*     IF (e.GT.0.995) THEN
*     WRITE (6,9)  (edot(i),i=1,3),(zt*ehat(i),i=1,3)
*   9 FORMAT (' edot zt*ehat   ',1P,3E10.2,2X,3E10.2)
*     END IF
*
*       Include tidal and quadrupole terms for each star and also GR.
      do i=1,3
         edot(i) = edot(i) + zt*ehat(i) + (zq + zg)*qhat(i)
*        hdot(i) = hdot(i) + zh*hhat(i)
      enddo
*       Save scalar value of eccentricity derivative for decision-making.
      dedt = edot(1)*ehat(1)+edot(2)*ehat(2)+edot(3)*ehat(3)
*
*     IF (e.GT.0.995) THEN
*     EDF = edot(1)*ehat(1)+edot(2)*ehat(2)+edot(3)*ehat(3)
*     WRITE (6,3)  e, EKD, (EKD - EDF)/EKD, e/zt, e/zq, e/zg
*   3 FORMAT (' E EKD DED/ED TT TQ TGR  ',F9.5,1P,5E9.1)
*     END IF
*     IF (IGR.LT.10.AND.zg.GT.zq*(cq(1) + cq(2)).AND.e.GT.0.99) THEN
*     IF (e.GT.0.995) THEN
*     ZD = ABS(EKD)
*     ZD = MAX(ZD,zq*(cq(1)+cq(2)))
*     if (e.GT.0.05.AND.zg.GT.ZD) THEN
*         WRITE (6,4)  e, a, a*(1.0-e), zg, ZD, 6.28/zg
*   4     FORMAT ('  GR    E A QP ZGR QUAD P_gr ',
*    &                                F8.4,1P,5E10.2)
*         IGR = IGR + 1
*     END IF
*
*     TAU = e/sqrt(edot(1)**2+edot(2)**2+edot(3)**2)
*     WRITE (6,5)  e, TAV, TAU, edot, zt, zq
*   5 FORMAT (' DERIV    e TAV TAU ed zt zt  ',F9.5,2F9.3,1P,5E10.2)
      ITIME = ITIME + 1
      IF (e.GT.0.980.AND.ITIME.EQ.-4) THEN
      EDT = zt
      EDQ = zq
      TF = e/EDT
      TQ = e/EDQ
      EDTOT = edot(1)*ehat(1)+edot(2)*ehat(2)+edot(3)*ehat(3)
      IF (ITIME.EQ.4) WRITE (6,6)  e0,TF,TQ,EDT,EDQ,EDTOT,EKD
    6 FORMAT (' DERIV    e TF TQ EDT EDQ ED EK ',F9.5,1P,6E10.2)
      END IF
      IF (ITIME.GE.4) ITIME = 0

*       Copy derivatives for Runge-Kutta integrator.
      do i=1,3
         du(i)=edot(i)
         du(i+3)=hdot(i)
      enddo
 10   return
      end