dpi1_0.for

      end do
      do k = 1, nclo
        iclo(k) = indexs(iclo(k))
        jclo(k) = indexs(jclo(k))
      end do
      end subroutine dpi_bhkce
!******************************************************************************
!     Dynamical Plug-In for Mercury 6
!     Subroutine to integrate with Bulirsch-Stoer algorithm the Keplerian
!     propagation of the orbits of those bodies involved in close encounters
!     (derived from MDT_BS2 subroutine of Mercury 6.2)
!     Adapted by: Diego Turrini
!     Last modified: July 2009
!     N.B.: the input variable FORCE (i.e. the subroutine to compute the
!           accelerations) has been removed. Subroutine dpi_bscea has
!           been instead used explicitly in the subroutine.
!******************************************************************************
!
!     Author: John E. Chambers (subroutine MDT_BS2 in Mercury 6.2)
!
!     Integrates NBOD bodies (of which NBIG are Big) for one timestep H0
!     using the Bulirsch-Stoer method. The accelerations are calculated using
!     the subroutine FORCE. The accuracy of the step is approximately
!     determined by the tolerance parameter TOL.
!
!     N.B.: This version only works for conservative systems (i.e. force is a
!           function of position only): non-gravitational forces and
!           post-Newtonian corrections cannot be used.
!
!******************************************************************************
      subroutine dpi_bs2(time,h0,hdid,tol,jcen,nbod,nbig,mass,x0,v0,s,
     %  rphys,rcrit,ngf,stat,dtflag,ngflag,opt,nce,ice,jce)
      implicit none
      include 'mercury.inc'
!     Parameters
      real*8 SHRINK,GROW
      parameter (SHRINK=.55d0,GROW=1.3d0)
!     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 dpi_bscea
!     Local
      integer j,j1,k,n,i,l
      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)
!
      tol2 = tol * tol
      do k=1,3
        do i=1,nbod
          xend(k,i)=0.d0
          a0(k,i)=0.d0
          a(k,i)=0.d0
        end do
      end do
!     Calculate arrays used to scale the relative error (R^2 for position and
!     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
!     Calculate accelerations at the start of the step
      call dpi_bscea(time,jcen,nbod,nbig,mass,x0,v0,s,rcrit,a0,stat,
     %        ngf,ngflag,opt,nce,ice,jce)
 100  continue
!     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)
        do k = 2, nbod
          do l=1,3
            b(l,k) = .5d0*a0(l,k)
            c(l,k) = 0.d0
            xend(l,k) = h2by2 * a0(l,k)  +  h * v0(l,k)  +  x0(l,k)
          end do
        end do
        do j = 2, n
          call dpi_bscea(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
            do l=1,3
              b(l,k) = b(l,k) + a(l,k)
              c(l,k) = c(l,k) + b(l,k)
              xend(l,k) = h2(n)*c(l,k) + h2by2*a0(l,k) + tmp0*v0(l,k)
     %             + x0(l,k)
            end do
          end do
        end do
        call dpi_bscea(time,jcen,nbod,nbig,mass,xend,v0,s,rcrit,a,
     %       stat,ngf,ngflag,opt,nce,ice,jce)
        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
!       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
!       After several integrations, test the relative error on extrapolated values
        if (n.gt.3) then
          errmax = 0.d0
!         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
!         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
            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
!           Save the actual stepsize used
            hdid = h0
!           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
      end do
!     If errors were too large, redo the step with half the previous step size.
      h0 = h0 * .5d0
      goto 100
      end subroutine dpi_bs2
!******************************************************************************
!     Dynamical Plug-In for Mercury 6
!     Subroutine to compute planetary accelerations for those bodies involved
!     in close encounters (derived from MFO_HKCE in Mercury 6.2).
!     Author: Diego Turrini
!     Last modified: July 2009
!******************************************************************************
!
!     Author: John E. Chambers (subroutine MFO_HKCE in Mercury 6.2)
!
!     Calculates accelerations due to the Keplerian part of the Hamiltonian
!     of a hybrid symplectic integrator, when close encounters are taking
!     place, for a set of NBOD bodies (NBIG of which are Big). Note that Small
!     Bodies do not interact with one another.
!
!******************************************************************************
      subroutine dpi_bscea(time,jcen,nbod,nbig,m,x,v,s,rcrit,a,stat,
     %  ngf,ngflag,opt,nce,ice,jce)
      implicit none
      include 'mercury.inc'
!     Input/Output
      integer nbod,nbig,stat(nbod),ngflag,opt(8),nce,ice(nce),jce(nce)
      real*8 time,jcen(3),rcrit(nbod),ngf(4,nbod),m(nbod)
      real*8 x(3,nbod),v(3,nbod),a(3,nbod),s(3,nbod)
!     Local Variables
      integer i,j,k
      real*8 r,q,q2,q3,q4,q5,tmp2
      real*8 rc,rc2,dx,dy,dz,faci,facj
      real*8 s_1,s_3,s2
      real*8 vsqr
      external vsqr
!     N.B.: here nbod=nbs from dpi_bhkce (not equal to integration nbod)
!     Variables initialization
      do j=1,nbod
        do i=1,3
          a(i,j)=0.d0
        end do
      end do
!     Computing planetary acceleration terms
      do k = 1, nce
        i = ice(k)
        j = jce(k)
        dx = x(1,j) - x(1,i)
        dy = x(2,j) - x(2,i)
        dz = x(3,j) - x(3,i)
        s2 = vsqr(dx,dy,dz)
        rc = max (rcrit(i), rcrit(j))
        rc2 = rc * rc
        if (s2.lt.rc2) then
          s_1 = 1.d0 / sqrt(s2)
          s_3 = s_1 * s_1 * s_1
          if (s2.le.0.01*rc2) then
            tmp2 = s_3
          else
            r = 1.d0 / s_1
            q = (r - 0.1d0*rc) / (0.9d0 * rc)
            q2 = q * q
            q3 = q * q2
            q4 = q2 * q2
            q5 = q2 * q3
            tmp2 = (1.d0 - 10.d0*q3 + 15.d0*q4 - 6.d0*q5) * s_3
          end if
          faci = tmp2 * m(i)
          facj = tmp2 * m(j)
          a(1,j) = a(1,j)  -  faci * dx
          a(2,j) = a(2,j)  -  faci * dy
          a(3,j) = a(3,j)  -  faci * dz
          a(1,i) = a(1,i)  +  facj * dx
          a(2,i) = a(2,i)  +  facj * dy
          a(3,i) = a(3,i)  +  facj * dz
        end if
      end do
c Solar terms
      do j = 2, nbod
        s2 = vsqr(x(1,j),x(2,j),x(3,j))
        s_1 = 1.d0 / sqrt(s2)
        tmp2 = m(1) * s_1 * s_1 * s_1
        do i= 1, 3
          a(i,j) = a(i,j)  -  tmp2 * x(i,j)
        end do
      end do
      end subroutine dpi_bscea
!******************************************************************************
!     Dynamical Plug-In for Mercury 6
!     Subroutine to identify those objects involved in close encounters with
!     other massive bodies during a time step (derived from subroutine MCE_SNIF
!     in Mercury 6.2 and adapted to S-type binary systems)
!     Author: Diego Turrini
!     Last modified: July 2009
!******************************************************************************
!
!     Author: John E. Chambers (subroutine MCE_SNIF in Mercury 6.2)
!
!     Given initial and final coordinates and velocities X and V and a timestep
!     H, the routine estimates which objects were involved in a close encounter
!     during the timestep.
!     The routine examines all objects with indices I >= I0
!
!     Returns an array CE, which for each object is:
!                           0 if it will undergo no encounters
!                           2 if it will pass within RCRIT of a Big body
!
!     Also returns arrays ICE and JCE, containing the indices of each pair of
!     objects estimated to have undergone an encounter.
!
!     N.B. All coordinates must be with respect to the central body
!
!******************************************************************************
      subroutine dpi_snif(h,i0,nbod,nbig,x0,v0,x1,v1,rcrit,ce,nce,ice,
     %  jce)
      implicit none
      include 'mercury.inc'
!     Input/Output
      integer i0,nbod,nbig,ce(nbod),nce,ice(NMAX),jce(NMAX)
      real*8 x0(3,nbod),v0(3,nbod),x1(3,nbod),v1(3,nbod),h,rcrit(nbod)
!     Local
      integer i,j
      real*8 d0,d1,d0t,d1t,d2min,temp,tmin,rc,rc2
      real*8 dx0,dy0,dz0,du0,dv0,dw0,dx1,dy1,dz1,du1,dv1,dw1
      real*8 xmin(NMAX),xmax(NMAX),ymin(NMAX),ymax(NMAX)
!
      if (i0.le.0) i0 = 2
      nce = 0
      do j = 2, nbod
        ce(j) = 0
      end do
!     Calculate maximum and minimum values of x and y coordinates
      call mce_box (nbod,h,x0,v0,x1,v1,xmin,xmax,ymin,ymax)
!     Adjust values for the Big bodies by symplectic close-encounter distance
!     N.B.:
      do j = i0, nbig-1
        xmin(j) = xmin(j) - rcrit(j)
        xmax(j) = xmax(j) + rcrit(j)
        ymin(j) = ymin(j) - rcrit(j)
        ymax(j) = ymax(j) + rcrit(j)
      end do
!     Identify pairs whose X-Y boxes overlap, and calculate minimum separation
!     First loop (over massive bodies)
      do i = i0, nbig-1
        do j = i + 1, nbig-1
          if (xmax(i).ge.xmin(j).and.xmax(j).ge.xmin(i)
     %      .and.ymax(i).ge.ymin(j).and.ymax(j).ge.ymin(i)) then
!           Determine the maximum separation that would qualify as an encounter
            rc = max(rcrit(i), rcrit(j))
            rc2 = rc * rc
!           Calculate initial and final separations
            dx0 = x0(1,i) - x0(1,j)
            dy0 = x0(2,i) - x0(2,j)
            dz0 = x0(3,i) - x0(3,j)
            dx1 = x1(1,i) - x1(1,j)
            dy1 = x1(2,i) - x1(2,j)
            dz1 = x1(3,i) - x1(3,j)
            d0 = dx0*dx0 + dy0*dy0 + dz0*dz0
            d1 = dx1*dx1 + dy1*dy1 + dz1*dz1
!           Check for a possible minimum in between
            du0 = v0(1,i) - v0(1,j)
            dv0 = v0(2,i) - v0(2,j)
            dw0 = v0(3,i) - v0(3,j)
            du1 = v1(1,i) - v1(1,j)
            dv1 = v1(2,i) - v1(2,j)
            dw1 = v1(3,i) - v1(3,j)
            d0t = (dx0*du0 + dy0*dv0 + dz0*dw0) * 2.d0
            d1t = (dx1*du1 + dy1*dv1 + dz1*dw1) * 2.d0
!           If separation derivative changes sign, find the minimum separation
            d2min = HUGE
            if (d0t*h.le.0.and.d1t*h.ge.0) call mce_min (d0,d1,d0t,d1t,
     %        h,d2min,tmin)
!           If minimum separation is small enough, flag this as a possible encounter
            temp = min (d0,d1,d2min)
            if (temp.le.rc2) then
              ce(i) = 2
              ce(j) = 2
              nce = nce + 1
              ice(nce) = i
              jce(nce) = j
            end if
          end if
        end do
      end do
!     Second loop (over massless particles)
      do i = i0, nbig-1
        do j = nbig+1, nbod
          if (xmax(i).ge.xmin(j).and.xmax(j).ge.xmin(i)
     %      .and.ymax(i).ge.ymin(j).and.ymax(j).ge.ymin(i)) then
!           Determine the maximum separation that would qualify as an encounter
            rc = max(rcrit(i), rcrit(j))
            rc2 = rc * rc
!           Calculate initial and final separations
            dx0 = x0(1,i) - x0(1,j)
            dy0 = x0(2,i) - x0(2,j)
            dz0 = x0(3,i) - x0(3,j)
            dx1 = x1(1,i) - x1(1,j)
            dy1 = x1(2,i) - x1(2,j)
            dz1 = x1(3,i) - x1(3,j)
            d0 = dx0*dx0 + dy0*dy0 + dz0*dz0
            d1 = dx1*dx1 + dy1*dy1 + dz1*dz1
!           Check for a possible minimum in between
            du0 = v0(1,i) - v0(1,j)
            dv0 = v0(2,i) - v0(2,j)
            dw0 = v0(3,i) - v0(3,j)
            du1 = v1(1,i) - v1(1,j)
            dv1 = v1(2,i) - v1(2,j)
            dw1 = v1(3,i) - v1(3,j)
            d0t = (dx0*du0 + dy0*dv0 + dz0*dw0) * 2.d0
            d1t = (dx1*du1 + dy1*dv1 + dz1*dw1) * 2.d0
!           If separation derivative changes sign, find the minimum separation
            d2min = HUGE
            if (d0t*h.le.0.and.d1t*h.ge.0) call mce_min (d0,d1,d0t,d1t,
     %        h,d2min,tmin)
!           If minimum separation is small enough, flag this as a possible encounter
            temp = min (d0,d1,d2min)
            if (temp.le.rc2) then
              ce(i) = 2
              ce(j) = 2
              nce = nce + 1
              ice(nce) = i
              jce(nce) = j
            end if
          end if
        end do
      end do
      end subroutine dpi_snif
!******************************************************************************
!     Dynamical Plug-In for Mercury 6
!     Subroutine to resolve collisions between bodies (not the stars) in S-type
!     binary systems (derived from subroutine MCE_COLL in Mercury 6.2)
!     N.B.: only perfectly inelastic collisions are allowed in present version
!     Author: Diego Turrini
!     Last modified: August 2009
!******************************************************************************
!
!     Author: John E. Chambers
!
!     Resolves a collision between two objects, using the collision model
!     chose by the user. Also writes a message to the information file, and
!     updates the value of ELOST, the change in energy due to collisions and
!     ejections.
!
!     N.B. All coordinates and velocities must be with respect to central body.
!
!******************************************************************************
      subroutine dpi_coll (time,tstart,elost,jcen,i,j,nbod,nbig,m,xh,
     %  vh,s,rphys,stat,id,opt,mem,lmem,outfile)
      implicit none
      include 'mercury.inc'
!     Input/Output
      integer i,j,nbod,nbig,stat(nbod),opt(8),lmem(NMESS)
      real*8 time,tstart,elost,jcen(3)
      real*8 m(nbod),xh(3,nbod),vh(3,nbod),s(3,nbod),rphys(nbod)
      character*80 outfile,mem(NMESS)
      character*8 id(nbod)
!     Local
      integer year,month,itmp
      real*8 t1
      character*38 flost,fcol
      character*6 tstring
      external dpi_merg
!     If two bodies collided, check that the less massive one is removed
!     (unless the more massive one is a Small body)
      if (i.ne.0) then
        if (m(j).gt.m(i).and.j.le.nbig) then
          itmp = i
          i = j
          j = itmp
        end if
      end if
!     Write message to info file (I=0 implies collision with the central body)
  10  open (23, file=outfile, status='old', access='append', err=10)
      if (opt(3).eq.1) then
        call mio_jd2y (time,year,month,t1)
        if (i.eq.0) then
          flost = '(1x,a8,a,i10,1x,i2,1x,f8.5)'
          write (23,flost) id(j),mem(67)(1:lmem(67)),year,month,t1
        else
          fcol  = '(1x,a8,a,a8,a,i10,1x,i2,1x,f4.1)'
          write (23,fcol) id(i),mem(69)(1:lmem(69)),id(j),
     %      mem(71)(1:lmem(71)),year,month,t1
        end if
      else
        if (opt(3).eq.3) then
          t1 = (time - tstart) / 365.25d0
          tstring = mem(2)
          flost = '(1x,a8,a,f18.7,a)'
          fcol  = '(1x,a8,a,a8,a,1x,f14.3,a)'
        else
          if (opt(3).eq.0) t1 = time
          if (opt(3).eq.2) t1 = time - tstart
          tstring = mem(1)
          flost = '(1x,a8,a,f18.5,a)'
          fcol  = '(1x,a8,a,a8,a,1x,f14.1,a)'
        end if
        if (i.eq.0.or.i.eq.1) then
          write (23,flost) id(j),mem(67)(1:lmem(67)),t1,tstring
        else
          write (23,fcol) id(i),mem(69)(1:lmem(69)),id(j),
     %      mem(71)(1:lmem(71)),t1,tstring
        end if
      end if
      close (23)
!     Do the collision (inelastic merger)
      call dpi_merg (jcen,i,j,nbod,nbig,m,xh,vh,s,stat,elost)
      end subroutine dpi_coll
!******************************************************************************
!     Dynamical Plug-In for Mercury 6
!     Subroutine to inelastically merge bodies involved in collisions (not the
!     stars) in S-type binary systems (derived from subroutine MCE_MERG in
!     Mercury 6.2)
!     N.B.: since this subroutine internally uses the initial inertial
!           coordinates the primary star is treated as any other body
!     Author: Diego Turrini
!     Last modified: August 2009
!******************************************************************************
!
!     Author: John E. Chambers
!
!     Merges objects I and J inelastically to produce a single new body by
!     conserving mass and linear momentum.
!     If J <= NBIG, then J is a Big body
!     If J >  NBIG, then J is a Small body
!     If I = 0, then I is the central body
!
!     N.B. All coordinates and velocities must be with respect to central body.
!
!******************************************************************************
      subroutine dpi_merg (jcen,i,j,nbod,nbig,m,x,v,s,stat,elost)
      implicit none
      include 'mercury.inc'
!     Input/Output
      integer i, j, nbod, nbig, stat(nbod)
      real*8 jcen(3),m(nbod),x(3,nbod),v(3,nbod),s(3,nbod),elost
!     Local
      real*8 tmp1, tmp2, dx, dy, dz, du, dv, dw, msum, mredu, msum_1
      real*8 xh(3,nbod), vh(3,nbod)
      external dpi_h2wb,dpi_wb2h
!     Checking collisions with central body
      if (i.le.1) i=1
!     Change back to inertial reference system
      call dpi_wb2h(nbig,nbod,m,xh,vh,x,v)
!     Resolve collisions
      msum   = m(i) + m(j)
      msum_1 = 1.d0 / msum
      mredu  = m(i) * m(j) * msum_1
      dx = xh(1,i) - xh(1,j)
      dy = xh(2,i) - xh(2,j)
      dz = xh(3,i) - xh(3,j)
      du = vh(1,i) - vh(1,j)
      dv = vh(2,i) - vh(2,j)
      dw = vh(3,i) - vh(3,j)
!     Calculate energy loss due to the collision
      elost = elost  +  .5d0 * mredu * (du*du + dv*dv + dw*dw)
     %        -  m(i) * m(j) / sqrt(dx*dx + dy*dy + dz*dz)
!     Calculate spin angular momentum of the new body
      s(1,i) = s(1,i)  +  s(1,j)  +  mredu * (dy * dw  -  dz * dv)
      s(2,i) = s(2,i)  +  s(2,j)  +  mredu * (dz * du  -  dx * dw)
      s(3,i) = s(3,i)  +  s(3,j)  +  mredu * (dx * dv  -  dy * du)
!     Calculate new coords and velocities by conserving centre of mass & momentum
      tmp1 = m(i) * msum_1
      tmp2 = m(j) * msum_1
      xh(1,i) = xh(1,i) * tmp1  +  xh(1,j) * tmp2
      xh(2,i) = xh(2,i) * tmp1  +  xh(2,j) * tmp2
      xh(3,i) = xh(3,i) * tmp1  +  xh(3,j) * tmp2
      vh(1,i) = vh(1,i) * tmp1  +  vh(1,j) * tmp2
      vh(2,i) = vh(2,i) * tmp1  +  vh(2,j) * tmp2
      vh(3,i) = vh(3,i) * tmp1  +  vh(3,j) * tmp2
      m(i) = msum
!     Flag the lost body for removal, and move it away from the new body
      stat(j) = -2
      xh(1,j) = -xh(1,j)
      xh(2,j) = -xh(2,j)
      xh(3,j) = -xh(3,j)
      vh(1,j) = -vh(1,j)
      vh(2,j) = -vh(2,j)
      vh(3,j) = -vh(3,j)
      m(j)   = 0.d0
      s(1,j) = 0.d0
      s(2,j) = 0.d0
      s(3,j) = 0.d0
!     Update coordinates in S-type reference system
      call dpi_h2wb(nbig,nbod,m,xh,vh,x,v)
      end subroutine dpi_merg
!******************************************************************************
!     Dynamical Plug-In for Mercury 6
!     Subroutine for the numerical integration of S-type binary stars systems
!     (derived from MAL-HCON subroutine of Mercury 6.2)
!     Adapted by: Diego Turrini
!     Last modified: August 2009
!******************************************************************************
!
!     Author: J. E. Chambers (subroutine MAL_HCON in Mercury 6.2)
!
!     Does an integration using an integrator with a constant stepsize H.
!     Input and output to this routine use coordinates XH, and velocities VH,
!     with respect to the central body, but the integration algorithm uses
!     its own internal coordinates X, and velocities V.
!
!******************************************************************************
      subroutine dpi_devbhc(time,tstart,tstop,dtout,algor,h0,tol,jcen,
     %  rcen,rmax,en,am,cefac,ndump,nfun,nbod,nbig,m,xh,vh,s,rho,rceh,
     %  stat,id,ngf,opt,opflag,ngflag,outfile,dumpfile,mem,lmem,nbin)
      implicit none
      include 'mercury.inc'
!     Input/Output variables
      integer algor,nbod,nbig,stat(nbod),opt(8),opflag,ngflag
      integer lmem(NMESS),ndump,nfun,nbin
      real*8 time,tstart,tstop,dtout,h0,tol,jcen(3),rcen,rmax
      real*8 en(3),am(3),cefac,m(nbod),xh(3,nbod),vh(3,nbod)
      real*8 s(3,nbod),rho(nbod),rceh(nbod),ngf(4,nbod)
      character*8 id(nbod)
      character*80 outfile(3),dumpfile(4),mem(NMESS)
!     Local variables
      integer i,j,k,n,itmp,nclo,nhit,jhit(CMAX),iclo(CMAX),jclo(CMAX)
      integer dtflag,ejflag,stopflag,colflag,nstored
      real*8 x(3,NMAX),v(3,NMAX),x0(3,NMAX),v0(3,NMAX)
      real*8 xh0(3,NMAX),vh0(3,NMAX),jac(nbod)
      real*8 rce(NMAX),rphys(NMAX),rcrit(NMAX),epoch(NMAX)
      real*8 hby2,tout,tmp0,tdump,tfun,tlog,dtdump,dtfun
      real*8 dclo(CMAX),tclo(CMAX),dhit(CMAX),thit(CMAX)
      real*8 ixvclo(6,CMAX),jxvclo(6,CMAX),a(NMAX)
      external dpi_wbstep,dpi_en,dpi_ejec,dpi_coll
!     Initialize variables. DTFLAG = 0/2: first call ever/normal
      dtout  = abs(dtout)
      dtdump = abs(h0) * ndump
      dtfun  = abs(h0) * nfun
      dtflag = 0
      nstored = 0
      hby2 = 0.5d0 * abs(h0)
      stopflag=0
!     Initialize state vectors to be used in the numerical integration
!     N.B.: after the first time step, the state vectors x and v do not express
!           anymore heliocentric coordinates and velocities. The variables x
!           and v are measured in a inertial coordinate system centered on the
!           origin O of the cartesian axes but the primary star (i.e. the Sun)
!           moves respect to the origin O (i.e. it orbits the center of mass of
!           the system). Heliocentric xh and vh are updated after each time step.
      do j=1,nbod
        do i=1,3
          x(i,j)=xh(i,j)
          v(i,j)=vh(i,j)
          jac(j)=0.d0
        end do
      end do
!     Calculate close-encounter limits and physical radii
!     N.B.: heliocentric state vectors should be used with this subroutine
      call mce_init (tstart,algor,h0,jcen,rcen,rmax,cefac,nbod,nbig,
     %  m,xh,vh,s,rho,rceh,rphys,rce,rcrit,id,opt,outfile(2),1)
!     Set up time of next output, times of previous dump, log and periodic effect
      if (opflag.eq.-1) then
        tout = tstart
      else
        n = int (abs (time-tstart) / dtout) + 1
        tout = tstart  +  dtout * sign (dble(n), tstop - tstart)
        if ((tstop-tstart)*(tout-tstop).gt.0) tout = tstop
      end if
      tdump = time
      tfun  = time
      tlog  = time
!
!     MAIN  LOOP  STARTS  HERE
!
 100  continue
!     Is it time for output ?
      if (abs(tout-time).le.hby2.and.opflag.ge.-1) then
!     Beware: the integration may change direction at this point!!!!
        if (opflag.eq.-1.and.dtflag.ne.0) dtflag = 1
!       Output data for all bodies
!       N.B.: heliocentric state vectors should be used with this subroutine
        call mio_out (time,jcen,rcen,rmax,nbod,nbig,m,xh,vh,s,rho,
     %    stat,id,opt,opflag,algor,outfile(1))
        tmp0 = tstop - tout
        tout = tout + sign( min( abs(tmp0), abs(dtout) ), tmp0 )
!       Update the data dump files
        do j = 2, nbod
          epoch(j) = time
        end do
!       N.B.: heliocentric state vectors should be used with this subroutine
        call mio_dump (time,tstart,tstop,dtout,algor,h0,tol,jcen,rcen,
     %    rmax,en,am,cefac,ndump,nfun,nbod,nbig,m,xh,vh,s,rho,rceh,stat,
     %    id,ngf,epoch,opt,opflag,dumpfile,mem,lmem)
        tdump = time
      end if
!     If integration has finished, return
      if (abs(tstop-time).le.hby2.and.opflag.ge.0) then
!     Copy inertial state vectors in output variables for the last
!     computation of energy and angular momentum done by the program
        do j=1,nbod
          do i=1,3
            xh(i,j)=x(i,j)
            vh(i,j)=v(i,j)
          end do
        end do
        return
      end if
!     Make sure the integration is heading in the right direction
 150  continue
      tmp0 = tstop - time
      if (opflag.eq.-1) tmp0 = tstart - time
      h0 = sign (h0, tmp0)
!     Save the current inertial coordinates and velocities
!     N.B.: inertial state vectors should be used with this subroutine
      call mco_iden (time,jcen,nbod,nbig,h0,m,x,v,x0,v0,ngf,
     %ngflag,opt)
!     Advance N-Body system of one time step
!     N.B.: inertial state vectors should be used with this subroutine
      if (algor.eq.12) then
        call dpi_wbstep(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)
      else
!       Check that only the implemented algorithm is used
        stop
      end if
!     Update heliocentric state vectors from the new inertial ones
      do j=1,nbod
        do i=1,3
          xh(i,j)=x(i,j)-x(i,1)
          vh(i,j)=v(i,j)-v(i,1)
        end do
      end do
      time = time + h0
!
!     CLOSE ENCOUNTERS AND COLLISIONS
!
!     CLOSE  ENCOUNTERS
!     If encounter minima occurred, output details and decide whether to stop
      if (nclo.gt.0.and.opflag.ge.-1) then
        itmp = 1
        if (colflag.ne.0) itmp = 0
        if (stopflag.eq.1) return
      end if
!     COLLISIONS
!     If collisions occurred, output details and remove lost objects
      if (colflag.ne.0) then
!       Reindex the surviving objects
!       N.B.: inertial state vectors should be used with this subroutine
        call mxx_elim (nbod,nbig,m,x,v,s,rho,rceh,rcrit,ngf,stat,
     %    id,mem,lmem,outfile(3),itmp)
!       Update heliocentric state vectors from the new inertial ones
        do j=1,nbod
          do i=1,3
            xh(i,j)=x(i,j)-x(i,1)
            vh(i,j)=v(i,j)-v(i,1)
          end do
        end do
!       Reset flags, and calculate new Hill radii and physical radii
!       N.B.: heliocentric state vectors should be used with this subroutine
        dtflag = 1
        if (opflag.ge.0) opflag = 1
        call mce_init (tstart,algor,h0,jcen,rcen,rmax,cefac,nbod,nbig,
     %    m,xh,vh,s,rho,rceh,rphys,rce,rcrit,id,opt,outfile(2),1)
      end if
!     COLLISIONS  WITH  CENTRAL  BODY
!     Check for collisions with the central body
      itmp = 2
      if (algor.eq.11) itmp = 3
!     Compute old heliocentric state vectors from the saved inertial ones
      do j=1,nbod
        do i=1,3
          xh0(i,j)=x0(i,j)-x0(i,1)
          vh0(i,j)=v0(i,j)-v0(i,1)
        end do
      end do
!     N.B.: heliocentric state vectors should be used with this subroutine
      call mce_cent (time,h0,rcen,jcen,itmp,nbod,nbig,m,xh0,vh0,xh,
     %     vh,nhit,jhit,thit,dhit,algor,ngf,ngflag)
!     If something hit the central body, restore the inertial coords prior to this step
      if (nhit.gt.0) then
         call mco_iden (time,jcen,nbod,nbig,h0,m,x0,v0,x,v,ngf,
     %        ngflag,opt)
         time = time - h0
!        Merge the object(s) with the central body
         do k = 1, nhit
            i = 1
            j = jhit(k)
!           N.B.: inertial state vectors should be used with this subroutine
            call dpi_coll (thit(k),tstart,en(3),jcen,i,j,nbod,nbig,m,x,
     %           v,s,rphys,stat,id,opt,mem,lmem,outfile(3))
        end do
!       Remove lost objects, reset flags and recompute Hill and physical radii
!       N.B.: inertial state vectors should be used with this subroutine
        call mxx_elim (nbod,nbig,m,x,v,s,rho,rceh,rcrit,ngf,stat,
     %    id,mem,lmem,outfile(3),itmp)
!       Update heliocentric state vectors from the new inertial ones
        do j=1,nbod
          do i=1,3
            xh(i,j)=x(i,j)-x(i,1)
            vh(i,j)=v(i,j)-v(i,1)
          end do
        end do
        if (opflag.ge.0) opflag = 1
        dtflag = 1
!       N.B.: heliocentric state vectors should be used with this subroutine
        call mce_init (tstart,algor,h0,jcen,rcen,rmax,cefac,nbod,nbig,
     %    m,xh,vh,s,rho,rceh,rphys,rce,rcrit,id,opt,outfile(2),0)
!       N.B.: inertial state vectors should be used with this subroutine
        call mco_iden (time,jcen,nbod,nbig,h0,m,x,v,x0,v0,ngf,
     %   ngflag,opt)
!       Redo that integration time step
        goto 150
      end if
!
!     I/O OPERATIONS
!
!     DATA  DUMP  AND  PROGRESS  REPORT
!     Do the data dump
      if (abs(time-tdump).ge.abs(dtdump).and.opflag.ge.-1) then
        do j = 2, nbod
          epoch(j) = time
        end do
!       N.B.: heliocentric state vectors should be used with this subroutine
        call mio_dump (time,tstart,tstop,dtout,algor,h0,tol,jcen,rcen,
     %    rmax,en,am,cefac,ndump,nfun,nbod,nbig,m,xh,vh,s,rho,rceh,stat,
     %    id,ngf,epoch,opt,opflag,dumpfile,mem,lmem)
        tdump = time
      end if
!     Update energy and angular momentum values and write a progress
!     report to the log file
      if (abs(time-tlog).ge.abs(dtdump).and.opflag.ge.0) then
!       N.B.: inertial state vectors should be used with this subroutine
        call dpi_en (nbod,nbig,m,x,v,en(2),am(2))
!         call mxx_jac(jcen,nbod,nbig,m,xh,vh,jac)
!         write(166,*) time/365.25d0,(jac(i),i=nbig+1,nbod)
!       N.B.: energy and angular momentum should be computed in the
!             S-type binary reference system
        call mio_log (time,tstart,en,am,opt,mem,lmem)
        tlog = time
      end if
!
!  EJECTIONS AND PERIODIC EFFECTS
!
      if (abs(time-tfun).ge.abs(dtfun).and.opflag.ge.-1) then
!       Recompute close encounter limits, to allow for changes in Hill radii
!       N.B.: heliocentric state vectors should be used with this subroutine
        call mce_hill (nbod,m,xh,vh,rce,a)
        do j = 2, nbod
          rce(j) = rce(j) * rceh(j)
        end do
!       Check for ejections
        ejflag=0
        itmp = 2
        if (algor.eq.11) itmp = 3
!       N.B.: inertial state vectors should be used with this subroutine
        call dpi_ejec (time,tstart,rmax,en,am,jcen,itmp,nbod,nbig,m,x,
     %    v,s,stat,id,opt,ejflag,outfile(3),mem,lmem)
!       Remove ejected objects, reset flags, calculate new Hill and physical radii
        if (ejflag.ne.0) then
!         N.B.: inertial state vectors should be used with this subroutine
          call mxx_elim (nbod,nbig,m,x,v,s,rho,rceh,rcrit,ngf,stat,
     %      id,mem,lmem,outfile(3),itmp)
!         Update heliocentric state vectors from the new inertial ones
          do j=1,nbod
            do i=1,3
              xh(i,j)=x(i,j)-x(i,1)
              vh(i,j)=v(i,j)-v(i,1)
            end do
          end do
!         N.B.: inertial state vectors should be used with this subroutine
          call dpi_en(nbod,nbig,m,x,v,en(2),am(2))
          if (opflag.ge.0) opflag = 1
          dtflag = 1
!         N.B.: heliocentric state vectors should be used with this subroutine
          call mce_init (tstart,algor,h0,jcen,rcen,rmax,cefac,nbod,nbig,
     %      m,xh,vh,s,rho,rceh,rphys,rce,rcrit,id,opt,outfile(2),0)
       end if
       tfun = time
      end if
!     Go on to the next time step
      goto 100
      end subroutine dpi_devbhc
!******************************************************************************
!     Dynamical Plug-In for Mercury 6
!     Subroutine for check massive bodies and test particles for ejections and
!     update the values of energy and angular momentum (derived from MXX_EJEC
!     subroutine of Mercury 6.2)
!     Adapted by: Diego Turrini
!     Last modified: August 2009
!******************************************************************************
!
!     Author: John E. Chambers
!
!     Calculates the distance from the central body of each object with index
!     I >= I0. If this distance exceeds RMAX, the object is flagged for
!     ejection (STAT set to -3). If any object is to be ejected, EJFLAG = 1 on
!     exit, otherwise EJFLAG = 0.
!
!     Also updates the values of EN(3) and AM(3)---the change in energy and
!     angular momentum due to collisions and ejections.
!
!     N.B. All coordinates must be with respect to the central body!!
!
!******************************************************************************
      subroutine dpi_ejec (time,tstart,rmax,en,am,jcen,i0,nbod,nbig,m,x,
     %  v,s,stat,id,opt,ejflag,outfile,mem,lmem)
      implicit none
      include 'mercury.inc'
!     Input/Output
      integer i0, nbod, nbig, stat(nbod), opt(8), ejflag, lmem(NMESS)
      real*8 time, tstart, rmax, en(3), am(3), jcen(3)
      real*8 m(nbod), x(3,nbod), v(3,nbod), s(3,nbod)
      character*80 outfile, mem(NMESS)
      character*8 id(nbod)
!     Local
      integer j, year, month
      real*8 r2,rmax2,t1,e,l
      character*38 flost
      character*6 tstring
      real*8 vsqr
      external dpi_en,vsqr
      if (i0.le.0) i0 = 2
      ejflag = 0
      rmax2 = rmax * rmax
!     Accessing output unit
  20  open  (23,file=outfile,status='old',access='append',err=20)
!     Calculate initial energy and angular momentum
      call dpi_en(nbod,nbig,m,x,v,e,l)
!     Flag each object which is ejected, and set its mass to zero
      if (i0.le.nbig) then
!       Processing both massive bodies and massless particles
!       Check massive bodies
        do j = i0, nbig-1
          r2 = vsqr(x(1,j)-x(1,1),x(2,j)-x(2,1),x(3,j)-x(3,1))
          if (r2.gt.rmax2) then
            ejflag = 1
            stat(j) = -3
            m(j) = 0.d0
            s(1,j) = 0.d0
            s(2,j) = 0.d0
            s(3,j) = 0.d0
!           Write message to information file
!   20        open  (23,file=outfile,status='old',access='append',err=20)
            if (opt(3).eq.1) then
              call mio_jd2y (time,year,month,t1)
              flost = '(1x,a8,a,i10,1x,i2,1x,f8.5)'
              write (23,flost) id(j),mem(68)(1:lmem(68)),year,month,t1
            else
              if (opt(3).eq.3) then
                t1 = (time - tstart) / 365.25d0
                tstring = mem(2)
                flost = '(1x,a8,a,f18.7,a)'
              else
                if (opt(3).eq.0) t1 = time
                if (opt(3).eq.2) t1 = time - tstart
                tstring = mem(1)
                flost = '(1x,a8,a,f18.5,a)'
              end if
              write (23,flost) id(j),mem(68)(1:lmem(68)),t1,tstring
            end if
!             close (23)
          end if
        end do
!       Check massless particles
        do j = nbig+1, nbod
          r2 = vsqr(x(1,j)-x(1,1),x(2,j)-x(2,1),x(3,j)-x(3,1))
          if (r2.gt.rmax2) then
            ejflag = 1
            stat(j) = -3
            m(j) = 0.d0
            s(1,j) = 0.d0
            s(2,j) = 0.d0
            s(3,j) = 0.d0
!           Write message to information file
!   20        open  (23,file=outfile,status='old',access='append',err=20)
            if (opt(3).eq.1) then
              call mio_jd2y (time,year,month,t1)
              flost = '(1x,a8,a,i10,1x,i2,1x,f8.5)'
              write (23,flost) id(j),mem(68)(1:lmem(68)),year,month,t1
            else
              if (opt(3).eq.3) then
                t1 = (time - tstart) / 365.25d0
                tstring = mem(2)
                flost = '(1x,a8,a,f18.7,a)'
              else
                if (opt(3).eq.0) t1 = time
                if (opt(3).eq.2) t1 = time - tstart
                tstring = mem(1)
                flost = '(1x,a8,a,f18.5,a)'
              end if
              write (23,flost) id(j),mem(68)(1:lmem(68)),t1,tstring
            end if
!             close (23)
          end if
        end do
      else
!       Processing only massless particles
        do j = i0, nbod
          r2 = vsqr(x(1,j)-x(1,1),x(2,j)-x(2,1),x(3,j)-x(3,1))
          if (r2.gt.rmax2) then
            ejflag = 1
            stat(j) = -3
            m(j) = 0.d0
            s(1,j) = 0.d0
            s(2,j) = 0.d0
            s(3,j) = 0.d0
!           Write message to information file
!   20        open  (23,file=outfile,status='old',access='append',err=20)
            if (opt(3).eq.1) then
              call mio_jd2y (time,year,month,t1)
              flost = '(1x,a8,a,i10,1x,i2,1x,f8.5)'
              write (23,flost) id(j),mem(68)(1:lmem(68)),year,month,t1
            else
              if (opt(3).eq.3) then
                t1 = (time - tstart) / 365.25d0
                tstring = mem(2)
                flost = '(1x,a8,a,f18.7,a)'
              else
                if (opt(3).eq.0) t1 = time
                if (opt(3).eq.2) t1 = time - tstart
                tstring = mem(1)
                flost = '(1x,a8,a,f18.5,a)'
              end if
              write (23,flost) id(j),mem(68)(1:lmem(68)),t1,tstring
            end if
!             close (23)
          end if
        end do
      end if
!     Exiting output unit
      close(23)
!     If ejections occurred, update ELOST and LLOST
      if (ejflag.ne.0) then
        call dpi_en(nbod,nbig,m,x,v,en(2),am(2))
        en(3) = en(3) + (e - en(2))
        am(3) = am(3) + (l - am(2))
      end if
      end subroutine dpi_ejec