mercury6_2-dpi.for

        if (abs(cf).gt.1) cf = sign(1.d0,cf)
        f = acos(cf)
        if (rv.lt.0) f = TWOPI - f
        p = true - f
        p = mod (p + TWOPI + TWOPI, TWOPI)
      end if
c
      if (l.lt.0) l = l + TWOPI
      if (l.gt.TWOPI) l = mod (l, TWOPI)
c
c------------------------------------------------------------------------------
c
      return
      end
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c      MDT_BS1.FOR    (ErikSoft   2 March 2001)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c Author: John E. Chambers
c
c Integrates NBOD bodies (of which NBIG are Big) for one timestep H0
c using the Bulirsch-Stoer method. The accelerations are calculated using the 
c subroutine FORCE. The accuracy of the step is approximately determined 
c by the tolerance parameter TOL.
c
c N.B. Input/output must be in coordinates with respect to the central body.
c ===
c
c------------------------------------------------------------------------------
c
      subroutine mdt_bs1 (time,h0,hdid,tol,jcen,nbod,nbig,mass,x0,v0,s,
     %  rphys,rcrit,ngf,stat,dtflag,ngflag,opt,nce,ice,jce,force)
c
      implicit none
      include 'mercury.inc'
c
      real*8 SHRINK,GROW
      parameter (SHRINK=.55d0,GROW=1.3d0)
c
c Input/Output
      integer nbod, nbig, opt(8), stat(nbod), dtflag, ngflag
      integer nce, ice(nce), jce(nce)
      real*8 time,h0,hdid,tol,jcen(3),mass(nbod),x0(3,nbod),v0(3,nbod)
      real*8 s(3,nbod),ngf(4,nbod),rphys(nbod),rcrit(nbod)
      external force
c
c Local
      integer j, j1, k, n
      real*8 tmp0,tmp1,tmp2,errmax,tol2,h,hx2,h2(8)
      real*8 x(3,NMAX),v(3,NMAX),xend(3,NMAX),vend(3,NMAX)
      real*8 a(3,NMAX),a0(3,NMAX),d(6,NMAX,8),xscal(NMAX),vscal(NMAX)
c
c------------------------------------------------------------------------------
c
      tol2 = tol * tol
c
c Calculate arrays used to scale the relative error (R^2 for position and
c V^2 for velocity).
      do k = 2, nbod
        tmp1 = x0(1,k)*x0(1,k) + x0(2,k)*x0(2,k) + x0(3,k)*x0(3,k)
        tmp2 = v0(1,k)*v0(1,k) + v0(2,k)*v0(2,k) + v0(3,k)*v0(3,k)
        xscal(k) = 1.d0 / tmp1
        vscal(k) = 1.d0 / tmp2
      end do
c
c Calculate accelerations at the start of the step
      call force (time,jcen,nbod,nbig,mass,x0,v0,s,rcrit,a0,stat,ngf,
     %  ngflag,opt,nce,ice,jce)
c
 100  continue
c
c For each value of N, do a modified-midpoint integration with 2N substeps
      do n = 1, 8
        h = h0 / (2.d0 * float(n))
        h2(n) = .25d0 / (n*n)
        hx2 = h * 2.d0
c
        do k = 2, nbod
          x(1,k) = x0(1,k) + h*v0(1,k)
          x(2,k) = x0(2,k) + h*v0(2,k)
          x(3,k) = x0(3,k) + h*v0(3,k)
          v(1,k) = v0(1,k) + h*a0(1,k)
          v(2,k) = v0(2,k) + h*a0(2,k)
          v(3,k) = v0(3,k) + h*a0(3,k)
        end do
        call force (time,jcen,nbod,nbig,mass,x,v,s,rcrit,a,stat,ngf,
     %    ngflag,opt,nce,ice,jce)
        do k = 2, nbod
          xend(1,k) = x0(1,k) + hx2*v(1,k)
          xend(2,k) = x0(2,k) + hx2*v(2,k)
          xend(3,k) = x0(3,k) + hx2*v(3,k)
          vend(1,k) = v0(1,k) + hx2*a(1,k)
          vend(2,k) = v0(2,k) + hx2*a(2,k)
          vend(3,k) = v0(3,k) + hx2*a(3,k)
        end do
c
        do j = 2, n
          call force (time,jcen,nbod,nbig,mass,xend,vend,s,rcrit,a,stat,
     %      ngf,ngflag,opt,nce,ice,jce)
          do k = 2, nbod
            x(1,k) = x(1,k) + hx2*vend(1,k)
            x(2,k) = x(2,k) + hx2*vend(2,k)
            x(3,k) = x(3,k) + hx2*vend(3,k)
            v(1,k) = v(1,k) + hx2*a(1,k)
            v(2,k) = v(2,k) + hx2*a(2,k)
            v(3,k) = v(3,k) + hx2*a(3,k)
          end do
          call force (time,jcen,nbod,nbig,mass,x,v,s,rcrit,a,stat,ngf,
     %      ngflag,opt,nce,ice,jce)
          do k = 2, nbod
            xend(1,k) = xend(1,k) + hx2*v(1,k)
            xend(2,k) = xend(2,k) + hx2*v(2,k)
            xend(3,k) = xend(3,k) + hx2*v(3,k)
            vend(1,k) = vend(1,k) + hx2*a(1,k)
            vend(2,k) = vend(2,k) + hx2*a(2,k)
            vend(3,k) = vend(3,k) + hx2*a(3,k)
          end do
        end do
c
        call force (time,jcen,nbod,nbig,mass,xend,vend,s,rcrit,a,stat,
     %    ngf,ngflag,opt,nce,ice,jce)
c
        do k = 2, nbod
          d(1,k,n) = .5d0*(xend(1,k) + x(1,k) + h*vend(1,k))
          d(2,k,n) = .5d0*(xend(2,k) + x(2,k) + h*vend(2,k))
          d(3,k,n) = .5d0*(xend(3,k) + x(3,k) + h*vend(3,k))
          d(4,k,n) = .5d0*(vend(1,k) + v(1,k) + h*a(1,k))
          d(5,k,n) = .5d0*(vend(2,k) + v(2,k) + h*a(2,k))
          d(6,k,n) = .5d0*(vend(3,k) + v(3,k) + h*a(3,k))
        end do
c
c Update the D array, used for polynomial extrapolation
        do j = n - 1, 1, -1
          j1 = j + 1
          tmp0 = 1.d0 / (h2(j) - h2(n))
          tmp1 = tmp0 * h2(j1)
          tmp2 = tmp0 * h2(n)
          do k = 2, nbod
            d(1,k,j) = tmp1 * d(1,k,j1)  -  tmp2 * d(1,k,j)
            d(2,k,j) = tmp1 * d(2,k,j1)  -  tmp2 * d(2,k,j)
            d(3,k,j) = tmp1 * d(3,k,j1)  -  tmp2 * d(3,k,j)
            d(4,k,j) = tmp1 * d(4,k,j1)  -  tmp2 * d(4,k,j)
            d(5,k,j) = tmp1 * d(5,k,j1)  -  tmp2 * d(5,k,j)
            d(6,k,j) = tmp1 * d(6,k,j1)  -  tmp2 * d(6,k,j)
          end do
        end do
c
c After several integrations, test the relative error on extrapolated values
        if (n.gt.3) then
          errmax = 0.d0
c
c Maximum relative position and velocity errors (last D term added)
          do k = 2, nbod
            tmp1 = max( d(1,k,1)*d(1,k,1), d(2,k,1)*d(2,k,1),
     %                  d(3,k,1)*d(3,k,1) )
            tmp2 = max( d(4,k,1)*d(4,k,1), d(5,k,1)*d(5,k,1),
     %                  d(6,k,1)*d(6,k,1) )
            errmax = max(errmax, tmp1*xscal(k), tmp2*vscal(k))
          end do
c
c If error is smaller than TOL, update position and velocity arrays, and exit
          if (errmax.le.tol2) then
            do k = 2, nbod
              x0(1,k) = d(1,k,1)
              x0(2,k) = d(2,k,1)
              x0(3,k) = d(3,k,1)
              v0(1,k) = d(4,k,1)
              v0(2,k) = d(5,k,1)
              v0(3,k) = d(6,k,1)
            end do
c
            do j = 2, n
              do k = 2, nbod
                x0(1,k) = x0(1,k) + d(1,k,j)
                x0(2,k) = x0(2,k) + d(2,k,j)
                x0(3,k) = x0(3,k) + d(3,k,j)
                v0(1,k) = v0(1,k) + d(4,k,j)
                v0(2,k) = v0(2,k) + d(5,k,j)
                v0(3,k) = v0(3,k) + d(6,k,j)
              end do
            end do
c
c Save the actual stepsize used
            hdid = h0
c
c Recommend a new stepsize for the next call to this subroutine
            if (n.eq.8) h0 = h0 * SHRINK
            if (n.lt.7) h0 = h0 * GROW
            return
          end if
        end if
c
      end do
c
c If errors were too large, redo the step with half the previous step size.
      h0 = h0 * .5d0
      goto 100
c
c------------------------------------------------------------------------------
c
      end
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c      MDT_BS2.FOR    (ErikSoft   2 March 2001)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c Author: John E. Chambers
c
c Integrates NBOD bodies (of which NBIG are Big) for one timestep H0
c using the Bulirsch-Stoer method. The accelerations are calculated using the 
c subroutine FORCE. The accuracy of the step is approximately determined 
c by the tolerance parameter TOL.
c
c N.B. This version only works for conservative systems (i.e. force is a
c ===  function of position only) !!!! Hence, non-gravitational forces
c      and post-Newtonian corrections cannot be used.
c
c N.B. Input/output must be in coordinates with respect to the central body.
c ===
c
c------------------------------------------------------------------------------
c
      subroutine mdt_bs2 (time,h0,hdid,tol,jcen,nbod,nbig,mass,x0,v0,s,
     %  rphys,rcrit,ngf,stat,dtflag,ngflag,opt,nce,ice,jce,force)
c
      implicit none
      include 'mercury.inc'
c
      real*8 SHRINK,GROW
      parameter (SHRINK=.55d0,GROW=1.3d0)
c
c Input/Output
      integer nbod, nbig, opt(8), stat(nbod), dtflag, ngflag
      real*8 time,h0,hdid,tol,jcen(3),mass(nbod),x0(3,nbod),v0(3,nbod)
      real*8 s(3,nbod),ngf(4,nbod),rphys(nbod),rcrit(nbod)
      integer nce,ice(nce),jce(nce)
      external force
c
c Local
      integer j,j1,k,n
      real*8 tmp0,tmp1,tmp2,errmax,tol2,h,h2(12),hby2,h2by2
      real*8 xend(3,NMAX),b(3,NMAX),c(3,NMAX)
      real*8 a(3,NMAX),a0(3,NMAX),d(6,NMAX,12),xscal(NMAX),vscal(NMAX)
c
c------------------------------------------------------------------------------
c
      tol2 = tol * tol
c
c Calculate arrays used to scale the relative error (R^2 for position and
c V^2 for velocity).
      do k = 2, nbod
        tmp1 = x0(1,k)*x0(1,k) + x0(2,k)*x0(2,k) + x0(3,k)*x0(3,k)
        tmp2 = v0(1,k)*v0(1,k) + v0(2,k)*v0(2,k) + v0(3,k)*v0(3,k)
        xscal(k) = 1.d0 / tmp1
        vscal(k) = 1.d0 / tmp2
      end do
c
c Calculate accelerations at the start of the step
      call force (time,jcen,nbod,nbig,mass,x0,v0,s,rcrit,a0,stat,ngf,
     %  ngflag,opt,nce,ice,jce)
c
 100  continue
c
c For each value of N, do a modified-midpoint integration with N substeps
      do n = 1, 12
        h = h0 / (dble(n))
        hby2  = .5d0 * h
        h2(n) = h * h
        h2by2 = .5d0 * h2(n)
c
        do k = 2, nbod
          b(1,k) = .5d0*a0(1,k)
          b(2,k) = .5d0*a0(2,k)
          b(3,k) = .5d0*a0(3,k)
          c(1,k) = 0.d0
          c(2,k) = 0.d0
          c(3,k) = 0.d0
          xend(1,k) = h2by2 * a0(1,k)  +  h * v0(1,k)  +  x0(1,k)
          xend(2,k) = h2by2 * a0(2,k)  +  h * v0(2,k)  +  x0(2,k)
          xend(3,k) = h2by2 * a0(3,k)  +  h * v0(3,k)  +  x0(3,k)
        end do
c
        do j = 2, n
          call force (time,jcen,nbod,nbig,mass,xend,v0,s,rcrit,a,stat,
     %      ngf,ngflag,opt,nce,ice,jce)
          tmp0 = h * dble(j)
          do k = 2, nbod
            b(1,k) = b(1,k) + a(1,k)
            b(2,k) = b(2,k) + a(2,k)
            b(3,k) = b(3,k) + a(3,k)
            c(1,k) = c(1,k) + b(1,k)
            c(2,k) = c(2,k) + b(2,k)
            c(3,k) = c(3,k) + b(3,k)
            xend(1,k) = h2(n)*c(1,k) + h2by2*a0(1,k) + tmp0*v0(1,k)
     %                + x0(1,k)
            xend(2,k) = h2(n)*c(2,k) + h2by2*a0(2,k) + tmp0*v0(2,k)
     %                + x0(2,k)
            xend(3,k) = h2(n)*c(3,k) + h2by2*a0(3,k) + tmp0*v0(3,k)
     %                + x0(3,k)
          end do
        end do
c
        call force (time,jcen,nbod,nbig,mass,xend,v0,s,rcrit,a,stat,ngf,
     %    ngflag,opt,nce,ice,jce)
c
        do k = 2, nbod
          d(1,k,n) = xend(1,k)
          d(2,k,n) = xend(2,k)
          d(3,k,n) = xend(3,k)
          d(4,k,n) = h*b(1,k) + hby2*a(1,k) + v0(1,k)
          d(5,k,n) = h*b(2,k) + hby2*a(2,k) + v0(2,k)
          d(6,k,n) = h*b(3,k) + hby2*a(3,k) + v0(3,k)
        end do
c
c Update the D array, used for polynomial extrapolation
        do j = n - 1, 1, -1
          j1 = j + 1
          tmp0 = 1.d0 / (h2(j) - h2(n))
          tmp1 = tmp0 * h2(j1)
          tmp2 = tmp0 * h2(n)
          do k = 2, nbod
            d(1,k,j) = tmp1 * d(1,k,j1)  -  tmp2 * d(1,k,j)
            d(2,k,j) = tmp1 * d(2,k,j1)  -  tmp2 * d(2,k,j)
            d(3,k,j) = tmp1 * d(3,k,j1)  -  tmp2 * d(3,k,j)
            d(4,k,j) = tmp1 * d(4,k,j1)  -  tmp2 * d(4,k,j)
            d(5,k,j) = tmp1 * d(5,k,j1)  -  tmp2 * d(5,k,j)
            d(6,k,j) = tmp1 * d(6,k,j1)  -  tmp2 * d(6,k,j)
          end do
        end do
c
c After several integrations, test the relative error on extrapolated values
        if (n.gt.3) then
          errmax = 0.d0
c
c Maximum relative position and velocity errors (last D term added)
          do k = 2, nbod
            tmp1 = max( d(1,k,1)*d(1,k,1), d(2,k,1)*d(2,k,1),
     %                  d(3,k,1)*d(3,k,1) )
            tmp2 = max( d(4,k,1)*d(4,k,1), d(5,k,1)*d(2,k,1),
     %                  d(6,k,1)*d(6,k,1) )
            errmax = max( errmax, tmp1*xscal(k), tmp2*vscal(k) )
          end do
c
c If error is smaller than TOL, update position and velocity arrays and exit
          if (errmax.le.tol2) then
            do k = 2, nbod
              x0(1,k) = d(1,k,1)
              x0(2,k) = d(2,k,1)
              x0(3,k) = d(3,k,1)
              v0(1,k) = d(4,k,1)
              v0(2,k) = d(5,k,1)
              v0(3,k) = d(6,k,1)
            end do
c
            do j = 2, n
              do k = 2, nbod
                x0(1,k) = x0(1,k) + d(1,k,j)
                x0(2,k) = x0(2,k) + d(2,k,j)
                x0(3,k) = x0(3,k) + d(3,k,j)
                v0(1,k) = v0(1,k) + d(4,k,j)
                v0(2,k) = v0(2,k) + d(5,k,j)
                v0(3,k) = v0(3,k) + d(6,k,j)
              end do
            end do
c
c Save the actual stepsize used
            hdid = h0
c
c Recommend a new stepsize for the next call to this subroutine
            if (n.ge.8) h0 = h0 * SHRINK
            if (n.lt.7) h0 = h0 * GROW
            return
          end if
        end if
c
      end do
c
c If errors were too large, redo the step with half the previous step size.
      h0 = h0 * .5d0
      goto 100
c
c------------------------------------------------------------------------------
c
      end
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c      MDT_HY.FOR    (ErikSoft   2 March 2001)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c Author: John E. Chambers
c
c Integrates NBOD bodies (of which NBIG are Big) for one timestep H
c using a second-order hybrid-symplectic integrator algorithm
c
c DTFLAG = 0 implies first ever call to this subroutine, 
c        = 1 implies first call since number/masses of objects changed.
c        = 2 normal call
c
c N.B. Input/output must be in democratic heliocentric coordinates.
c ===
c
c------------------------------------------------------------------------------
c
      subroutine mdt_hy (time,tstart,h0,tol,rmax,en,am,jcen,rcen,nbod,
     %  nbig,m,x,v,s,rphys,rcrit,rce,stat,id,ngf,algor,opt,dtflag,
     %  ngflag,opflag,colflag,nclo,iclo,jclo,dclo,tclo,ixvclo,jxvclo,
     %  outfile,mem,lmem)
c
      implicit none
      include 'mercury.inc'
c
c Input/Output
      integer nbod,nbig,stat(nbod),algor,opt(8),dtflag,ngflag,opflag
      integer colflag,lmem(NMESS),nclo,iclo(CMAX),jclo(CMAX)
      real*8 time,tstart,h0,tol,rmax,en(3),am(3),jcen(3),rcen
      real*8 m(nbod),x(3,nbod),v(3,nbod),s(3,nbod),rphys(nbod)
      real*8 rce(nbod),rcrit(nbod),ngf(4,nbod),tclo(CMAX),dclo(CMAX)
      real*8 ixvclo(6,CMAX),jxvclo(6,CMAX)
      character*80 outfile(3),mem(NMESS)
      character*8 id(nbod)
c
c Local
      integer j,nce,ice(NMAX),jce(NMAX),ce(NMAX),iflag
      real*8 a(3,NMAX),hby2,hrec,x0(3,NMAX),v0(3,NMAX),mvsum(3),temp
      real*8 angf(3,NMAX),ausr(3,NMAX)
      external mfo_hkce
c
c------------------------------------------------------------------------------
c
      save a, hrec, angf, ausr
      hby2 = h0 * .5d0
      nclo = 0
      colflag = 0
c
c If accelerations from previous call are not valid, calculate them now
      if (dtflag.ne.2) then
        if (dtflag.eq.0) hrec = h0
        call mfo_hy (jcen,nbod,nbig,m,x,rcrit,a,stat)
        dtflag = 2
        do j = 2, nbod
          angf(1,j) = 0.d0
          angf(2,j) = 0.d0
          angf(3,j) = 0.d0
          ausr(1,j) = 0.d0
          ausr(2,j) = 0.d0
          ausr(3,j) = 0.d0
        end do
c If required, apply non-gravitational and user-defined forces
        if (opt(8).eq.1) call mfo_user (time,jcen,nbod,nbig,m,x,v,ausr)
        if (ngflag.eq.1.or.ngflag.eq.3) call mfo_ngf (nbod,x,v,angf,ngf)
      end if
c
c Advance interaction Hamiltonian for H/2
      do j = 2, nbod
        v(1,j) = v(1,j)  +  hby2 * (angf(1,j) + ausr(1,j) + a(1,j))
        v(2,j) = v(2,j)  +  hby2 * (angf(2,j) + ausr(2,j) + a(2,j))
        v(3,j) = v(3,j)  +  hby2 * (angf(3,j) + ausr(3,j) + a(3,j))
      end do
c
c Advance solar Hamiltonian for H/2
      mvsum(1) = 0.d0
      mvsum(2) = 0.d0
      mvsum(3) = 0.d0
      do j = 2, nbod
        mvsum(1) = mvsum(1)  +  m(j) * v(1,j)
        mvsum(2) = mvsum(2)  +  m(j) * v(2,j)
        mvsum(3) = mvsum(3)  +  m(j) * v(3,j)
      end do
c
      temp = hby2 / m(1)
      mvsum(1) = temp * mvsum(1)
      mvsum(2) = temp * mvsum(2)
      mvsum(3) = temp * mvsum(3)
      do j = 2, nbod
        x(1,j) = x(1,j)  +  mvsum(1)
        x(2,j) = x(2,j)  +  mvsum(2)
        x(3,j) = x(3,j)  +  mvsum(3)
      end do
c
c Save the current coordinates and velocities
      call mco_iden (time,jcen,nbod,nbig,h0,m,x,v,x0,v0,ngf,ngflag,opt)
c
c Advance H_K for H
      do j = 2, nbod
        iflag = 0
        call drift_one (m(1),x(1,j),x(2,j),x(3,j),v(1,j),v(2,j),
     %    v(3,j),h0,iflag)
      end do
c
c Check whether any object separations were < R_CRIT whilst advancing H_K
      call mce_snif (h0,2,nbod,nbig,x0,v0,x,v,rcrit,ce,nce,ice,jce)
c
c If objects had close encounters, advance H_K using Bulirsch-Stoer instead
      if (nce.gt.0) then
        do j = 2, nbod
          if (ce(j).ne.0) then
            x(1,j) = x0(1,j)
            x(2,j) = x0(2,j)
            x(3,j) = x0(3,j)
            v(1,j) = v0(1,j)
            v(2,j) = v0(2,j)
            v(3,j) = v0(3,j)
          end if
        end do
        call mdt_hkce (time,tstart,h0,hrec,tol,rmax,en(3),jcen,rcen,
     %    nbod,nbig,m,x,v,s,rphys,rcrit,rce,stat,id,ngf,algor,opt,
     %    ngflag,colflag,ce,nce,ice,jce,nclo,iclo,jclo,dclo,tclo,ixvclo,
     %    jxvclo,outfile,mem,lmem,mfo_hkce)
      end if
c
c Advance solar Hamiltonian for H/2
      mvsum(1) = 0.d0
      mvsum(2) = 0.d0
      mvsum(3) = 0.d0
      do j = 2, nbod
        mvsum(1) = mvsum(1)  +  m(j) * v(1,j)
        mvsum(2) = mvsum(2)  +  m(j) * v(2,j)
        mvsum(3) = mvsum(3)  +  m(j) * v(3,j)
      end do
c
      temp = hby2 / m(1)
      mvsum(1) = temp * mvsum(1)
      mvsum(2) = temp * mvsum(2)
      mvsum(3) = temp * mvsum(3)
      do j = 2, nbod
        x(1,j) = x(1,j)  +  mvsum(1)
        x(2,j) = x(2,j)  +  mvsum(2)
        x(3,j) = x(3,j)  +  mvsum(3)
      end do
c
c Advance interaction Hamiltonian for H/2
      call mfo_hy (jcen,nbod,nbig,m,x,rcrit,a,stat)
      if (opt(8).eq.1) call mfo_user (time,jcen,nbod,nbig,m,x,v,ausr)
      if (ngflag.eq.1.or.ngflag.eq.3) call mfo_ngf (nbod,x,v,angf,ngf)
c
      do j = 2, nbod
        v(1,j) = v(1,j)  +  hby2 * (angf(1,j) + ausr(1,j) + a(1,j))
        v(2,j) = v(2,j)  +  hby2 * (angf(2,j) + ausr(2,j) + a(2,j))
        v(3,j) = v(3,j)  +  hby2 * (angf(3,j) + ausr(3,j) + a(3,j))
      end do
c
c------------------------------------------------------------------------------
c
      return
      end
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c      MDT_HKCE.FOR    (ErikSoft   1 March 2001)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c Author: John E. Chambers
c
c Integrates NBOD bodies (of which NBIG are Big) for one timestep H under
c the Hamiltonian H_K, including close-encounter terms.
c
c------------------------------------------------------------------------------
c
      subroutine mdt_hkce (time,tstart,h0,hrec,tol,rmax,elost,jcen,
     %  rcen,nbod,nbig,m,x,v,s,rphy,rcrit,rce,stat,id,ngf,algor,opt,
     %  ngflag,colflag,ce,nce,ice,jce,nclo,iclo,jclo,dclo,tclo,ixvclo,
     %  jxvclo,outfile,mem,lmem,force)
c
      implicit none
      include 'mercury.inc'
c
c Input/Output
      integer nbod,nbig,nce,ice(nce),jce(nce),stat(nbod),ngflag,ce(nbod)
      integer algor,opt(8),colflag,lmem(NMESS),nclo,iclo(CMAX)
      integer jclo(CMAX)
      real*8 time,tstart,h0,hrec,tol,rmax,elost,jcen(3),rcen
      real*8 m(nbod),x(3,nbod),v(3,nbod),s(3,nbod)
      real*8 rce(nbod),rphy(nbod),rcrit(nbod),ngf(4,nbod)
      real*8 tclo(CMAX),dclo(CMAX),ixvclo(6,CMAX),jxvclo(6,CMAX)
      character*80 outfile(3),mem(NMESS)
      character*8 id(nbod)
      external force
c
c Local
      integer iback(NMAX),index(NMAX),ibs(NMAX),jbs(NMAX),nclo_old
      integer i,j,k,nbs,nbsbig,statbs(NMAX)
      integer nhit,ihit(CMAX),jhit(CMAX),chit(CMAX),nowflag,dtflag
      real*8 tlocal,hlocal,hdid,tmp0
      real*8 mbs(NMAX),xbs(3,NMAX),vbs(3,NMAX),sbs(3,NMAX)
      real*8 rcritbs(NMAX),rcebs(NMAX),rphybs(NMAX)
      real*8 ngfbs(4,NMAX),x0(3,NMAX),v0(3,NMAX)
      real*8 thit(CMAX),dhit(CMAX),thit1,temp
      character*8 idbs(NMAX)
c
c------------------------------------------------------------------------------
c
c N.B. Don't set nclo to zero!!
      nbs = 1
      nbsbig = 0
      mbs(1) = m(1)
      if (algor.eq.11) mbs(1) = m(1) + m(2)
      sbs(1,1) = s(1,1)
      sbs(2,1) = s(2,1)
      sbs(3,1) = s(3,1)
c
c Put data for close-encounter bodies into local arrays for use with BS routine
      do j = 2, nbod
        if (ce(j).ne.0) then
          nbs = nbs + 1
          if (j.le.nbig) nbsbig = nbs
          mbs(nbs)   = m(j)
          xbs(1,nbs) = x(1,j)
          xbs(2,nbs) = x(2,j)
          xbs(3,nbs) = x(3,j)
          vbs(1,nbs) = v(1,j)
          vbs(2,nbs) = v(2,j)
          vbs(3,nbs) = v(3,j)
          sbs(1,nbs) = s(1,j)
          sbs(2,nbs) = s(2,j)
          sbs(3,nbs) = s(3,j)
          rcebs(nbs) = rce(j)
          rphybs(nbs) = rphy(j)
          statbs(nbs) = stat(j)
          rcritbs(nbs) = rcrit(j)
          idbs(nbs) = id(j)
          index(nbs) = j
          iback(j) = nbs
        end if
      end do
c
      do k = 1, nce
        ibs(k) = iback(ice(k))
        jbs(k) = iback(jce(k))
      end do
c
      tlocal = 0.d0
      hlocal = sign(hrec,h0)
c
c Begin the Bulirsch-Stoer integration
  50  continue
        tmp0 = abs(h0) - abs(tlocal)
        hrec = hlocal
        if (abs(hlocal).gt.tmp0) hlocal = sign (tmp0, h0)
c
c Save old coordinates and integrate
        call mco_iden (time,jcen,nbs,0,h0,mbs,xbs,vbs,x0,v0,ngf,ngflag,
     %    opt)
        call mdt_bs2 (time,hlocal,hdid,tol,jcen,nbs,nbsbig,mbs,xbs,vbs,
     %    sbs,rphybs,rcritbs,ngfbs,statbs,dtflag,ngflag,opt,nce,
     %    ibs,jbs,force)
        tlocal = tlocal + hdid
c
c Check for close-encounter minima
        nclo_old = nclo
        temp = time + tlocal
        call mce_stat (temp,hdid,rcen,nbs,nbsbig,mbs,x0,v0,xbs,vbs,
     %    rcebs,rphybs,nclo,iclo,jclo,dclo,tclo,ixvclo,jxvclo,nhit,ihit,
     %    jhit,chit,dhit,thit,thit1,nowflag,statbs,outfile(3),mem,lmem)
c
c If collisions occurred, resolve the collision and return a flag
        if (nhit.gt.0.and.opt(2).ne.0) then
          do k = 1, nhit
            if (chit(k).eq.1) then
              i = ihit(k)
              j = jhit(k)
              call mce_coll (thit(k),tstart,elost,jcen,i,j,nbs,nbsbig,
     %          mbs,xbs,vbs,sbs,rphybs,statbs,idbs,opt,mem,lmem,
     %          outfile(3))
              colflag = colflag + 1
            end if
          end do
        end if
c
c If necessary, continue integrating objects undergoing close encounters
      if ((tlocal - h0)*h0.lt.0) goto 50
c
c Return data for the close-encounter objects to global arrays
      do k = 2, nbs
        j = index(k)
        m(j)   = mbs(k)
        x(1,j) = xbs(1,k)
        x(2,j) = xbs(2,k)
        x(3,j) = xbs(3,k)
        v(1,j) = vbs(1,k)
        v(2,j) = vbs(2,k)
        v(3,j) = vbs(3,k)
        s(1,j) = sbs(1,k)
        s(2,j) = sbs(2,k)
        s(3,j) = sbs(3,k)
        stat(j) = statbs(k)
      end do
      do k = 1, nclo
        iclo(k) = index(iclo(k))
        jclo(k) = index(jclo(k))
      end do
c
c------------------------------------------------------------------------------
c
      return
      end
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c      MDT_MVS.FOR    (ErikSoft   28 March 2001)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c Author: John E. Chambers
c
c Integrates NBOD bodies (of which NBIG are Big) for one timestep H
c using a second-order mixed-variable symplectic integrator.
c
c DTFLAG = 0 implies first ever call to this subroutine, 
c        = 1 implies first call since number/masses of objects changed.
c        = 2 normal call
c
c N.B. Input/output must be in coordinates with respect to the central body.
c ===
c
c------------------------------------------------------------------------------
c
      subroutine mdt_mvs (time,tstart,h0,tol,rmax,en,am,jcen,rcen,nbod,
     %  nbig,m,x,v,s,rphys,rcrit,rce,stat,id,ngf,algor,opt,dtflag,
     %  ngflag,opflag,colflag,nclo,iclo,jclo,dclo,tclo,ixvclo,jxvclo,
     %  outfile,mem,lmem)
c
      implicit none
      include 'mercury.inc'
c
c Input/Output
      integer nbod,nbig,stat(nbod),algor,opt(8),dtflag,ngflag,opflag
      integer colflag,lmem(NMESS),nclo,iclo(CMAX),jclo(CMAX)
      real*8 time,tstart,h0,tol,rmax,en(3),am(3),jcen(3),rcen
      real*8 m(nbod),x(3,nbod),v(3,nbod),s(3,nbod),rphys(nbod)
      real*8 rce(nbod),rcrit(nbod),ngf(4,nbod),tclo(CMAX),dclo(CMAX)
      real*8 ixvclo(6,CMAX),jxvclo(6,CMAX)
      character*80 outfile(3),mem(NMESS)
      character*8 id(nbod)
c
c Local
      integer j,iflag,nhit,ihit(CMAX),jhit(CMAX),chit(CMAX),nowflag
      real*8 xj(3,NMAX),vj(3,NMAX),a(3,NMAX),gm(NMAX),hby2,thit1,temp
      real*8 msofar,minside,x0(3,NMAX),v0(3,NMAX),dhit(CMAX),thit(CMAX)
      real*8 angf(3,NMAX),ausr(3,NMAX)
c
c------------------------------------------------------------------------------
c
      save a, xj, gm, angf, ausr
      hby2 = .5d0 * h0
      nclo = 0
c
c If accelerations from previous call are not valid, calculate them now,
c and also the Jacobi coordinates XJ, and effective central masses GM.
      if (dtflag.ne.2) then
        dtflag = 2
        call mco_h2j (time,jcen,nbig,nbig,h0,m,x,v,xj,vj,ngf,ngflag,opt)
        call mfo_mvs (jcen,nbod,nbig,m,x,xj,a,stat)
c
        minside = m(1)
        do j = 2, nbig
          msofar = minside + m(j)
          gm(j) = m(1) * msofar / minside
          minside = msofar
          angf(1,j) = 0.d0
          angf(2,j) = 0.d0
          angf(3,j) = 0.d0
          ausr(1,j) = 0.d0
          ausr(2,j) = 0.d0
          ausr(3,j) = 0.d0
        end do
c If required, apply non-gravitational and user-defined forces
        if (opt(8).eq.1) call mfo_user (time,jcen,nbod,nbig,m,x,v,ausr)
        if (ngflag.eq.1.or.ngflag.eq.3) call mfo_ngf (nbod,x,v,angf,ngf)
      end if
c
c Advance interaction Hamiltonian for H/2
      do j = 2, nbod
        v(1,j) = v(1,j)  +  hby2 * (angf(1,j) + ausr(1,j) + a(1,j))
        v(2,j) = v(2,j)  +  hby2 * (angf(2,j) + ausr(2,j) + a(2,j))
        v(3,j) = v(3,j)  +  hby2 * (angf(3,j) + ausr(3,j) + a(3,j))
      end do
c
c Save current coordinates and velocities
      call mco_iden (time,jcen,nbod,nbig,h0,m,x,v,x0,v0,ngf,ngflag,opt)
c
c Advance Keplerian Hamiltonian (Jacobi/helio coords for Big/Small bodies)
      call mco_h2j (time,jcen,nbig,nbig,h0,m,x,v,xj,vj,ngf,ngflag,opt)
      do j = 2, nbig
        iflag = 0
        call drift_one (gm(j),xj(1,j),xj(2,j),xj(3,j),vj(1,j),
     %    vj(2,j),vj(3,j),h0,iflag)
      end do
      do j = nbig + 1, nbod
        iflag = 0
        call drift_one (m(1),x(1,j),x(2,j),x(3,j),v(1,j),v(2,j),
     %    v(3,j),h0,iflag)
      end do
      call mco_j2h (time,jcen,nbig,nbig,h0,m,xj,vj,x,v,ngf,ngflag,opt)
c
c Check for close-encounter minima during drift step
      temp = time + h0
      call mce_stat (temp,h0,rcen,nbod,nbig,m,x0,v0,x,v,rce,rphys,nclo,
     %  iclo,jclo,dclo,tclo,ixvclo,jxvclo,nhit,ihit,jhit,chit,dhit,thit,
     %  thit1,nowflag,stat,outfile(3),mem,lmem)
c
c Advance interaction Hamiltonian for H/2
      call mfo_mvs (jcen,nbod,nbig,m,x,xj,a,stat)
      if (opt(8).eq.1) call mfo_user (time,jcen,nbod,nbig,m,x,v,ausr)
      if (ngflag.eq.1.or.ngflag.eq.3) call mfo_ngf (nbod,x,v,angf,ngf)
c
      do j = 2, nbod
        v(1,j) = v(1,j)  +  hby2 * (angf(1,j) + ausr(1,j) + a(1,j))
        v(2,j) = v(2,j)  +  hby2 * (angf(2,j) + ausr(2,j) + a(2,j))
        v(3,j) = v(3,j)  +  hby2 * (angf(3,j) + ausr(3,j) + a(3,j))
      end do
c
c------------------------------------------------------------------------------
c
      return
      end
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c      MDT_RA15.FOR    (ErikSoft   2 March 2001)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c Author: John E. Chambers
c
c Integrates NBOD bodies (of which NBIG are Big) for one timestep H0 using
c Everhart's RA15 integrator algorithm. The accelerations are calculated
c using the subroutine FORCE. The accuracy of the step is approximately 
c determined by the tolerance parameter TOL.
c
c Based on RADAU by E. Everhart, Physics Department, University of Denver.
c Comments giving equation numbers refer to Everhart (1985) ``An Efficient
c Integrator that Uses Gauss-Radau Spacings'', in The Dynamics of Comets:
c Their Origin and Evolution, p185-202, eds. A. Carusi & G. B. Valsecchi,
c pub Reidel. (A listing of the original subroutine is also given in this 
c paper.)
c
c DTFLAG = 0 implies first ever call to this subroutine, 
c        = 1 implies first call since number/masses of objects changed.
c        = 2 normal call
c
c N.B. Input/output must be in coordinates with respect to the central body.
c ===
c
c------------------------------------------------------------------------------
c
      subroutine mdt_ra15 (time,t,tdid,tol,jcen,nbod,nbig,mass,x1,v1,
     %  spin,rphys,rcrit,ngf,stat,dtflag,ngflag,opt,nce,ice,jce,force)
c
      implicit none
      include 'mercury.inc'
c
c Input/Output
      integer nbod,nbig,dtflag,ngflag,opt(8),stat(nbod)
      integer nce,ice(nce),jce(nce)
      real*8 time,t,tdid,tol,jcen(3),mass(nbod)
      real*8 x1(3*nbod),v1(3*nbod),spin(3*nbod)
      real*8 ngf(4,nbod),rphys(nbod),rcrit(nbod)
      external force
c
c Local
      integer nv,niter,j,k,n
      real*8 x(3*NMAX),v(3*NMAX),a(3*NMAX),a1(3*NMAX)
      real*8 g(7,3*NMAX),b(7,3*NMAX),e(7,3*NMAX)
      real*8 h(8),xc(8),vc(7),c(21),d(21),r(28),s(9)
      real*8 q,q2,q3,q4,q5,q6,q7,temp,gk
c
c------------------------------------------------------------------------------
c
      save h,xc,vc,c,d,r,b,e
c
c Gauss-Radau spacings for substeps within a sequence, for the 15th order 
c integrator. The sum of the H values should be 3.733333333333333
c
      data h/          0.d0,.0562625605369221d0,.1802406917368924d0,
     %  .3526247171131696d0,.5471536263305554d0,.7342101772154105d0,
     %  .8853209468390958d0,.9775206135612875d0/
c
c Constant coefficients used in series expansions for X and V
c  XC: 1/2,  1/6,  1/12, 1/20, 1/30, 1/42, 1/56, 1/72
c  VC: 1/2,  1/3,  1/4,  1/5,  1/6,  1/7,  1/8
      data xc/.5d0,.1666666666666667d0,.08333333333333333d0,.05d0,
     %  .03333333333333333d0,.02380952380952381d0,.01785714285714286d0,
     %  .01388888888888889d0/
      data vc/.5d0,.3333333333333333d0,.25d0,.2d0,
     %  .1666666666666667d0,.1428571428571429d0,.125d0/
c
c If this is first call to the subroutine, set values of the constant arrays
c (R = R21, R31, R32, R41, R42, R43 in Everhart's paper.)
      if (dtflag.eq.0) then
        n = 0
        do j = 2, 8
          do k = 1, j - 1
            n = n + 1
            r(n) = 1.d0 / (h(j) - h(k))
          end do
        end do
c
c Constants to convert between B and G arrays (C = C21, C31, C32, C41, C42...)
        c(1) = - h(2)
        d(1) =   h(2)
        n = 1
        do j = 3, 7
          n = n + 1
          c(n) = -h(j) * c(n-j+2)
          d(n) =  h(2) * d(n-j+2)
          do k = 3, j - 1
            n = n + 1
            c(n) = c(n-j+1)  -  h(j) * c(n-j+2)
            d(n) = d(n-j+1)  +  h(k) * d(n-j+2)
          end do
          n = n + 1
          c(n) = c(n-j+1) - h(j)
          d(n) = d(n-j+1) + h(j)
        end do
c
        dtflag = 1
      end if
c
      nv = 3 * nbod
  100 continue
c
c If this is first call to subroutine since number/masses of objects changed
c do 6 iterations and initialize B, E arrays, otherwise do 2 iterations.
      if (dtflag.eq.1) then
        niter = 6
        do j = 4, nv
          do k = 1, 7
            b (k,j) = 0.d0
            e (k,j) = 0.d0
          end do
        end do
      else
        niter = 2
      end if
c
c Calculate forces at the start of the sequence
      call force (time,jcen,nbod,nbig,mass,x1,v1,spin,rcrit,a1,stat,ngf,
     %  ngflag,opt,nce,ice,jce)
c
c Find G values from B values predicted at the last call (Eqs. 7 of Everhart)
      do k = 4, nv
        g(1,k) = b(7,k)*d(16) + b(6,k)*d(11) + b(5,k)*d(7)
     %         + b(4,k)*d(4)  + b(3,k)*d(2)  + b(2,k)*d(1)  + b(1,k)
        g(2,k) = b(7,k)*d(17) + b(6,k)*d(12) + b(5,k)*d(8)
     %         + b(4,k)*d(5)  + b(3,k)*d(3)  + b(2,k)
        g(3,k) = b(7,k)*d(18) + b(6,k)*d(13) + b(5,k)*d(9)
     %         + b(4,k)*d(6)  + b(3,k)
        g(4,k) = b(7,k)*d(19) + b(6,k)*d(14) + b(5,k)*d(10) + b(4,k)
        g(5,k) = b(7,k)*d(20) + b(6,k)*d(15) + b(5,k)
        g(6,k) = b(7,k)*d(21) + b(6,k)
        g(7,k) = b(7,k)
      end do
c
c------------------------------------------------------------------------------
c
c  MAIN  LOOP  STARTS  HERE
c
c For each iteration (six for first call to subroutine, two otherwise)...
      do n = 1, niter
c
c For each substep within a sequence...
        do j = 2, 8
c
c Calculate position predictors using Eqn. 9 of Everhart
          s(1) = t * h(j)
          s(2) = s(1) * s(1) * .5d0
          s(3) = s(2) * h(j) * .3333333333333333d0
          s(4) = s(3) * h(j) * .5d0
          s(5) = s(4) * h(j) * .6d0
          s(6) = s(5) * h(j) * .6666666666666667d0
          s(7) = s(6) * h(j) * .7142857142857143d0
          s(8) = s(7) * h(j) * .75d0
          s(9) = s(8) * h(j) * .7777777777777778d0
c
          do k = 4, nv
            x(k) = s(9)*b(7,k) + s(8)*b(6,k) + s(7)*b(5,k)
     %           + s(6)*b(4,k) + s(5)*b(3,k) + s(4)*b(2,k)
     %           + s(3)*b(1,k) + s(2)*a1(k)  + s(1)*v1(k) + x1(k)
          end do
c
c If necessary, calculate velocity predictors too, from Eqn. 10 of Everhart
          if (ngflag.ne.0) then
            s(1) = t * h(j)
            s(2) = s(1) * h(j) * .5d0
            s(3) = s(2) * h(j) * .6666666666666667d0
            s(4) = s(3) * h(j) * .75d0
            s(5) = s(4) * h(j) * .8d0
            s(6) = s(5) * h(j) * .8333333333333333d0
            s(7) = s(6) * h(j) * .8571428571428571d0
            s(8) = s(7) * h(j) * .875d0
c
            do k = 4, nv