!******************************************************************************
!
!                                 DPI Library
!                       Dynamical Plug-In for Mercury 6
!                                Release 1.0
!
!                    Copyright (C) 2002-2015  Diego Turrini
!
!     This program is free software: you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation, either version 3 of the License, or
!     (at your option) any later version.
! 
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.
! 
!     You should have received a copy of the GNU General Public License
!     along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!------------------------------------------------------------------------------
!
!     DPI Library
!     Dynamical Plug-In for Mercury 6
!
!     Release: 1.0
!     Author: Diego Turrini
!     E-mail: diego.turrini_at_iaps.inaf.it
!     Last modified: January 2010
!     Modified by: Diego Turrini
!
!     Disclaimer: this library is supplied as a plug-in for the Mercury software
!                 developed by John E. Chamber: specifically, it is designed to 
!                 work with version 6 of Mercury. To facilitate its use in 
!                 conjunction with Mercury, the design of the subroutines was
!                 kept as similar as possible to that of analogous subroutines
!                 in Mercury. Mercury is developed by John E. Chambers and can
!                 be downloaded at the following URL:
!                 http://www.arm.ac.uk/~jec/home.html (see also the 
!                 Astrophysics Source Code Library, ascl id. 1201.008) 
!                 The DPI library contains some subroutines derived and adapted 
!                 from Mercury 6 and originally developed by John E. Chambers. 
!                 Information on the author and, where appropriate, the name of
!                 the original subroutines are reported in the headers.  
!
!     Acknowledgements: the author wishes to thank Patricia Eleanor Verrier
!                       for her assistance in debugging and restructuring the
!                       algorithm and the code through comparison with her
!                       MOIRAI software.
!
!     Bibliographic references:
!
!     Mercury 6 software:
!
!     Chambers J. E., 1999, Monthly Notices of the Royal Astronomical Society, 
!     vol. 304, pp. 793-799 
!
!     Symplectic mapping for S-type binary star systems:
!
!     Chambers J. E., Quintana E. V., Duncan M. J., Lissauer J. J., 2002, The
!     Astrophysical Journal, vol. 123, pp. 2884-2894
!
!     MOIRAI software:
! 
!     Verrier P. E., Evans N. W., 2007, Monthly Notices of the Royal 
!     Astronomical Society, vol. 382, pp. 1432-1446
!
!     DPI library:
!
!     Turrini D., Barbieri M., Marzari F., Thebault P., Tricarico P., 2005, 
!     Memorie della Societa' Astronomica Italiana Supplementi, vol. 6,
!     pp. 172-177
!
!     Turrini D., Barbieri M., Marzari F., Tricarico P., 2004, Memorie della
!     Societa' Astronomica Italiana Supplementi, vol. 5, pp. 127-130
!
!     Thebault P., Marzari F., Scholl H., Turrini D., Barbieri M., 2004, 
!     Astronomy & Astrophysics, vol. 427, pp. 1097-1104
!
!******************************************************************************
!
!******************************************************************************
!     Dynamical Plug-In for Mercury 6
!     Subroutine to convert from the initial, inertial coordinates to
!     S-type binary coordinates as in Chambers et al. (2002)
!     N.B.: the subroutine computes positions and pseudovelocities, not the
!           state vectors
!     N.B.: the initial coordinates are centered on the initial position of the
!           central star yet are not heliocentric
!     Author: Diego Turrini
!     Last modified: May 2009
!******************************************************************************
      subroutine dpi_h2wb(nbig,nbod,m,xh,vh,xb,pvb)
      implicit none
!     Input/Output variables
      integer nbod,nbig
      real*8 m(nbod)
      real*8 xh(3,nbod),vh(3,nbod)
      real*8 xb(3,nbod),pvb(3,nbod)
!     Local variables
      integer i,j
      real*8 mtot,mden,mvh(3),mxh(3)
!     Variables initialization
      do i=1,3
        mxh(i)=0.d0
        mvh(i)=0.d0
      end do
      do j=1,nbod
        do i=1,3
          xb(i,j)=0.d0
          pvb(i,j)=0.d0
        end do
      end do
!     Computing support variables
      mtot=m(1)
      if (nbig.gt.2) then
        do j=2,nbig-1
          mtot=mtot+m(j)
          do i=1,3
            mvh(i)=mvh(i)+m(j)*vh(i,j)
            mxh(i)=mxh(i)+m(j)*xh(i,j)
          end do
        end do
      end if
      mden=mtot
      mtot=mtot+m(nbig)
!     Computing transformed positions for S type binary systems
!     Computing position of central star
      do i=1,3
        xb(i,1)=(m(1)*xh(i,1)+m(nbig)*xh(i,nbig)+mxh(i))/mtot
      end do
!     If any, computing positions of massive particles
      if (nbig.gt.2) then
        do j=2,nbig-1
          do i=1,3
            xb(i,j)=xh(i,j)-xh(i,1)
          end do
        end do
      end if
!     Computing position of binary star
      do i=1,3
        xb(i,nbig)=xh(i,nbig)-(m(1)*xh(i,1)+mxh(i))/mden
      end do
!     If any, computing positions of massless particles
      if (nbod.gt.nbig) then
        do j=nbig+1,nbod
          do i=1,3
            xb(i,j)=xh(i,j)-xh(i,1)
          end do
        end do
      end if
!     Computing transformed pseudo-velocities for S type binary systems
!     Computing pseudo-velocity of central star
      do i=1,3
        pvb(i,1)=(m(1)*vh(i,1)+m(nbig)*vh(i,nbig)+mvh(i))/mtot
      end do
!     If any, computing pseudo-velocities of massive particles
      if (nbig.gt.2) then
        do j=2,nbig-1
          do i=1,3
            pvb(i,j)=vh(i,j)-(m(1)*vh(i,1)+mvh(i))/mden
          end do
        end do
      end if
!     Computing pseudo-velocity of binary star
      do i=1,3
        pvb(i,nbig)=vh(i,nbig)-pvb(i,1)
      end do
!     If any, computing pseudo-velocities of massless particles
      if (nbod.gt.nbig) then
        do j=nbig+1,nbod
          do i=1,3
            pvb(i,j)=vh(i,j)-(m(1)*vh(i,1)+mvh(i))/mden
          end do
        end do
      end if
      end subroutine dpi_h2wb
!******************************************************************************
!     Dynamical Plug-In for Mercury 6
!     Subroutine to convert from the S-type binary coordinates to the initial,
!     inertial coordinates as in Chambers et al. (2002)
!     N.B.: the subroutine needs positions and pseudovelocities, not the
!           state vectors
!     N.B.: the initial coordinates are centered on the initial position of the
!           central star yet are not heliocentric
!     Author: Diego Turrini
!     Last modified: January 2010
!******************************************************************************
      subroutine dpi_wb2h(nbig,nbod,m,xh,vh,xb,pvb)
      implicit none
!     Input/Output variables
      integer nbod,nbig
      real*8 m(nbod)
      real*8 xh(3,nbod),vh(3,nbod)
      real*8 xb(3,nbod),pvb(3,nbod)
!     Local variables
      integer i,j
      real*8 mtot,mden,mvb(3),mxb(3)
!     Variables initialization
      do i=1,3
        mxb(i)=0.d0
        mvb(i)=0.d0
      end do
      do j=1,nbod
        do i=1,3
          xh(i,j)=0.d0
          vh(i,j)=0.d0
        end do
      end do
!     Computing support variables
      mtot=m(1)
      if (nbig.gt.2) then
        do j=2,nbig-1
          mtot=mtot+m(j)
          do i=1,3
            mxb(i)=mxb(i)+m(j)*xb(i,j)
            mvb(i)=mvb(i)+m(j)*pvb(i,j)
          end do
        end do
      end if
      mden=mtot
      mtot=mtot+m(nbig)
!     Computing heliocentric positions for S type binary systems
!     Computing position of central star
      do i=1,3
        xh(i,1)=xb(i,1)-(m(nbig)/mtot)*xb(i,nbig)-mxb(i)/mden
      end do
!     If any, computing positions of massive particles
      if (nbig.gt.2) then
        do j=2,nbig-1
          do i=1,3
            xh(i,j)=xb(i,j)+xh(i,1)
          end do
        end do
      end if
!     Computing position of binary star
      do i=1,3
        xh(i,nbig)=xb(i,1)+(mden/mtot)*xb(i,nbig)
      end do
!     Computing positions of massless particles
      if (nbod.gt.nbig) then
        do j=nbig+1,nbod
          do i=1,3
            xh(i,j)=xb(i,j)+xh(i,1)
          end do
        end do
      end if
!     Computing heliocentric velocities for S type binary systems
!     Computing velocity of central star
      do i=1,3
        vh(i,1)=pvb(i,1)-(m(nbig)/mden)*pvb(i,nbig)-mvb(i)/m(1)
      end do
!     If any, computing velocities of massive particles
      if (nbig.gt.2) then
        do j=2,nbig-1
          do i=1,3
            vh(i,j)=pvb(i,j)+pvb(i,1)-(m(nbig)/mden)*pvb(i,nbig)
          end do
        end do
      end if
!     Computing velocity of binary star
      do i=1,3
        vh(i,nbig)=pvb(i,1)+pvb(i,nbig)
      end do
!     If any, computing velocities of massless particles
      if (nbod.gt.nbig) then
        do j=nbig+1,nbod
          do i=1,3
            vh(i,j)=pvb(i,j)+pvb(i,1)-(m(nbig)/mden)*pvb(i,nbig)
          end do
        end do
      end if
      end subroutine dpi_wb2h
!******************************************************************************
!     Dynamical Plug-In for Mercury 6
!     Subroutine to convert S-type velocities for Keplerian propagation back
!     to pseudovelocities (based on the subroutines developed by Patricia
!     Author: Diego Turrini
!     Last modified: May 2009
!******************************************************************************
      subroutine dpi_vb2psv_k(nbig,nbod,m,xb,vb,psv)
      implicit none
!     Input/Output variables
      integer nbod,nbig
      real*8 m(nbod)
      real*8 xb(3,nbod),vb(3,nbod),psv(3,nbod)
!     Local variables
      integer i,j
      real*8 mtot,mden
!     Variables initialization
      mtot=m(1)
      if (nbig.gt.2) then
        do j=2,nbig-1
          mtot=mtot+m(j)
        end do
      end if
      mden=mtot
      mtot=mtot+m(nbig)
!     Computing pseudo-velocity of central star
      do i=1,3
        psv(i,1)=vb(i,1)
      end do
!     If any, computing pseudo-velocities of massive particles
      if (nbig.gt.2) then
        do j=2,nbig-1
          do i=1,3
            psv(i,j)=vb(i,j)
          end do
        end do
      end if
!     Computing pseudo-velocity of binary star
      do i=1,3
        psv(i,nbig)=vb(i,nbig)/(mtot/mden)
      end do
!     If any, computing pseudo-velocities of massless particles
      if (nbod.gt.nbig) then
        do j=nbig+1,nbod
          do i=1,3
            psv(i,j)=vb(i,j)
          end do
        end do
      end if
      end subroutine dpi_vb2psv_k
!******************************************************************************
!     Dynamical Plug-In for Mercury 6
!     Subroutine to convert pseudovelocities to S-type velocities for Keplerian
!     propagation
!     Author: Diego Turrini
!     Last modified: May 2009
!******************************************************************************
      subroutine dpi_psv2vb_k(nbig,nbod,m,xb,vb,psv)
      implicit none
!     Input/Output variables
      integer nbod,nbig
      real*8 m(nbod)
      real*8 xb(3,nbod),vb(3,nbod),psv(3,nbod)
!     Local variables
      integer i,j
      real*8 mtot,mden
!     Variables initialization
      mtot=m(1)
      if (nbig.gt.2) then
        do j=2,nbig-1
          mtot=mtot+m(j)
        end do
      end if
      mden=mtot
      mtot=mtot+m(nbig)
!     Computing transformed velocity of central star
      do i=1,3
        vb(i,1)=psv(i,1)
      end do
!     If any, computing transformed velocities of massive particles
      if (nbig.gt.2) then
        do j=2,nbig-1
          do i=1,3
            vb(i,j)=psv(i,j)
          end do
        end do
      end if
!     Computing transformed velocity of binary star
      do i=1,3
        vb(i,nbig)=psv(i,nbig)*mtot/mden
      end do
!     If any, computing transformed velocities of massless particles
      if (nbod.gt.nbig) then
        do j=nbig+1,nbod
          do i=1,3
            vb(i,j)=psv(i,j)
          end do
        end do
      end if
      end subroutine dpi_psv2vb_k
!******************************************************************************
!     Dynamical Plug-In for Mercury 6
!     Subroutine to compute the H_jump part of the transformed Hamiltonian for
!     S-type binary systems
!     Author: Diego Turrini
!     Last modified: May 2009
!******************************************************************************
      subroutine dpi_wbjump(nbig,nbod,dt,m,xb,pvb)
      implicit none
!     Input/Output variables
      integer nbod,nbig
      real*8 m(nbod),dt
      real*8 xb(3,nbod),pvb(3,nbod)
!     Local variables
      integer i,j
      real*8 pbsum(3)
!     Variables initialization
      do i=1,3
        pbsum(i)=0.d0
      end do
      if (nbig.gt.2) then
        do j=2,nbig-1
          do i=1,3
            pbsum(i)=pbsum(i)+pvb(i,j)*m(j)
          end do
        end do
        do i=1,3
          pbsum(i)=dt*pbsum(i)/m(1)
        end do
!       Advancing massive particles
        do j=2,nbig-1
          do i=1,3
            xb(i,j)=xb(i,j)+pbsum(i)
          end do
        end do
      end if
!     Advancing massless particles
      if (nbod.gt.nbig) then
        do j=nbig+1,nbod
          do i=1,3
            xb(i,j)=xb(i,j)+pbsum(i)
          end do
        end do
      end if
      end subroutine dpi_wbjump
!******************************************************************************
!     Dynamical Plug-In for Mercury 6
!     Subroutine to compute the accelerations due to the H_int part of the
!     transformed Hamiltonian for S-type binary systems
!     Author: Diego Turrini
!     Last modified: January 2010
!******************************************************************************
      subroutine dpi_wbfor(nbig,nbod,xb,a,m,rcrit)
      implicit none
!     Input/Output variables
      integer nbod,nbig
      real*8 xb(3,nbod),m(nbod),a(3,nbod),rcrit(nbod)
!     Local variables
      integer i,j,l
      real*8 mden,mtot,s(3),den1,den2,den3(nbod)
      real*8 factor1(3),factor2(3,nbod),rb(nbod)
      real*8 drb(3),cenorm,r,r2,r3,rc,rc2,q,q3,q4,q5
      real*8 cube,vmod,vsqr
      external cube,vmod,vsqr
      intrinsic max,sqrt
!     Variables initialization
      do i=1,3
        factor1(i)=0.d0
        s(i)=0.d0
      end do
      do j=1,nbod
        rb(j)=0.d0
        do i=1,3
          a(i,j)=0.d0
          factor2(i,j)=0.d0
        end do
      end do
!     Computation of support variables
      mtot=m(1)
      if (nbig.gt.2) then
        do j=2,nbig-1
          mtot=mtot+m(j)
          do i=1,3
            s(i)=s(i)+m(j)*xb(i,j)
          end do
        end do
        ! N.B.: here mtot = mtot-m(nbod) = mden
        do i=1,3
          s(i)=s(i)/mtot
        end do
      end if
      mden=mtot
      mtot=mtot+m(nbig)
      do i=1,3
        factor1(i)=xb(i,nbig)+s(i)
      end do
      rb(nbig)=vmod(xb(1,nbig),xb(2,nbig),xb(3,nbig))
      den1=cube(rb(nbig))
      den2=cube(vmod(factor1(1),factor1(2),factor1(3)))
      if (nbig.gt.2) then
        do j=2,nbig-1
          do i=1,3
            factor2(i,j)=xb(i,nbig)-xb(i,j)+s(i)
          end do
          rb(j)=vmod(xb(1,j),xb(2,j),xb(3,j))
          den3(j)=cube(vmod(factor2(1,j),factor2(2,j),factor2(3,j)))
        end do
      end if
      if (nbod.gt.nbig) then
        do j=nbig+1,nbod
          do i=1,3
            factor2(i,j)=xb(i,nbig)-xb(i,j)+s(i)
          end do
          rb(j)=vmod(xb(1,j),xb(2,j),xb(3,j))
          den3(j)=cube(vmod(factor2(1,j),factor2(2,j),factor2(3,j)))
        end do
      end if
!     Computation of the acceleration of the binary companion
      do i=1,3
        a(i,nbig)=a(i,nbig)+m(1)*(xb(i,nbig)/den1)
        a(i,nbig)=a(i,nbig)-m(1)*factor1(i)/den2
      end do
      if (nbig.gt.2) then
        do j=2,nbig-1
          do i=1,3
            a(i,nbig)=a(i,nbig)+m(j)*(xb(i,nbig)/den1)
            a(i,nbig)=a(i,nbig)-m(j)*factor2(i,j)/den3(j)
          end do
        end do
      end if
!     Computation of the accelerations of the massive particles
      if (nbig.gt.2) then
        do j=2,nbig-1
          do i=1,3
            a(i,j)=a(i,j)-(m(1)*m(nbig)/mden)*(factor1(i)/den2)
            a(i,j)=a(i,j)+m(nbig)*factor2(i,j)/den3(j)
          end do
        end do
        do j=2,nbig-1
          do l=2,nbig-1
            do i=1,3
              a(i,j)=a(i,j)-(m(nbig)*m(l)/mden)*factor2(i,l)/den3(l)
            end do
          end do
        end do
      end if
      if (nbig.gt.3) then
        do j=2,nbig-1
          do l=j+1,nbig-1
            drb(1)=xb(1,j)-xb(1,l)
            drb(2)=xb(2,j)-xb(2,l)
            drb(3)=xb(3,j)-xb(3,l)
            r2=vsqr(drb(1),drb(2),drb(3))
            r=sqrt(r2)
            r3=r*r*r
            rc=max(rcrit(j),rcrit(l))
            rc2=rc*rc
            do i=1,3
              if (r2.ge.rc2) then
                cenorm=1.d0
              else if (r2.le.0.01d0*rc2) then
                cenorm=0.d0
              else
                q=(r-0.1d0*rc)/(0.9d0*rc)
                q3=q*q*q
                q4=q3*q
                q5=q4*q
                cenorm=(10.d0*q3-15.d0*q4+6.d0*q5)
              end if
              a(i,j)=a(i,j)-cenorm*m(l)*drb(i)/r3
              a(i,l)=a(i,l)+cenorm*m(j)*drb(i)/r3
            end do
          end do
        end do
      end if
!     Computation of the accelerations of the massless particles
      if (nbod.gt.nbig) then
        do j=nbig+1,nbod
          do i=1,3
            a(i,j)=a(i,j)-(m(1)*m(nbig)/mden)*(factor1(i)/den2)
            a(i,j)=a(i,j)+m(nbig)*factor2(i,j)/den3(j)
          end do
        end do
      end if
      if (nbig.gt.2.and.nbod.gt.nbig) then
        do j=nbig+1,nbod
          do l=2,nbig-1
            do i=1,3
              a(i,j)=a(i,j)-(m(nbig)*m(l)/mden)*factor2(i,l)/den3(l)
            end do
          end do
        end do
        do l=nbig+1,nbod
          do j=2,nbig-1
            drb(1)=xb(1,l)-xb(1,j)
            drb(2)=xb(2,l)-xb(2,j)
            drb(3)=xb(3,l)-xb(3,j)
            r2=vsqr(drb(1),drb(2),drb(3))
            r=sqrt(r2)
            r3=r*r*r
            rc=max(rcrit(j),rcrit(l))
            rc2=rc*rc
            do i=1,3
              if (r2.ge.rc2) then
                cenorm=1.d0
              else if (r2.le.0.01d0*rc2) then
                cenorm=0.d0
              else
                q=(r-0.1d0*rc)/(0.9d0*rc)
                q3=q*q*q
                q4=q3*q
                q5=q4*q
                cenorm=(10.d0*q3-15.d0*q4+6.d0*q5)
              end if
              a(i,l)=a(i,l)-cenorm*m(j)*drb(i)/r3
            end do
          end do
        end do
      end if
      end subroutine dpi_wbfor
!******************************************************************************
!     Dynamical Plug-In for Mercury 6
!     Subroutine to advance of one time step the N-Body system for S-type
!     binary systems
!     Author: Diego Turrini
!     Last modified: January 2010
!******************************************************************************
      subroutine dpi_wbstep(time,tstart,h0,tol,rmax,en,am,jcen,rcen,
     & nbod,nbig,m,xh,vh,s,rphys,rcrit,rce,stat,id,ngf,algor,opt,dtflag
     & ,ngflag,opflag,colflag,nclo,iclo,jclo,dclo,tclo,ixvclo,jxvclo,
     & outfile,mem,lmem)
      implicit none
      include 'mercury.inc'
!     Input/Output variables
      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),xh(3,nbod),vh(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)
!     Local variables
      integer i,j,iflag,nce,ice(nbod),jce(nbod),ce(nbod)
      real*8 a(3,NMAX),x0(3,nbod),v0(3,nbod),hrec,hby2
      real*8 xb(3,NMAX),pvb(3,NMAX),vb(3,NMAX)
      real*8 mden,mtot,gm(NMAX),mb(NMAX)
      external dpi_h2wb,dpi_wb2h,dpi_psv2vb_k,dpi_vb2psv_k
      external dpi_wbfor,dpi_wbjump,drift_one,dpi_bhkce
      save hrec,gm,mb,xb,pvb,a,mtot,mden
!     N.B.: this subroutine internally uses the type S (wide binary)
!           coordinate system described in Chambers et al. (2002)
      hby2=h0*0.5d0
      nclo = 0
      colflag = 0
      if (dtflag.ne.2) then
        if (dtflag.eq.0) hrec=h0
        call dpi_h2wb(nbig,nbod,m,xh,vh,xb,pvb)
!       Computing accelerations
        do j=1,nbod
          do i=1,3
            a(i,j)=0.d0
          end do
        end do
        if (nbod.gt.2) call dpi_wbfor(nbig,nbod,xb,a,m,rcrit)
        gm(1)=0.d0
        mb(1)=0.d0
        if (nbod.lt.3) then
          mden=m(1)
          mtot=m(1)+m(2)
          gm(2)=mtot
          mb(2)=m(nbod)*m(1)/mtot
        else
          mtot=m(1)
          do j=2,nbig-1
            mtot=mtot+m(j)
            gm(j)=m(1)
            mb(j)=m(j)
          end do
          mden=mtot
          mtot=mtot+m(nbig)
          gm(nbig)=mtot
          mb(nbig)=m(nbig)*mden/mtot
          do j=nbig+1,nbod
            gm(j)=m(1)
            mb(j)=m(j)
          end do
        end if
        dtflag = 2
      end if
      if (nbod.gt.2) then
!     Advancing the Interaction Hamiltonian
        do j=2,nbod
          do i=1,3
            pvb(i,j)=pvb(i,j)+hby2*a(i,j)
          end do
        end do
!       Advancing the Jump Hamiltonian
        call dpi_wbjump(nbig,nbod,hby2,m,xb,pvb)
      end if
!     Computing transformed velocities
      call dpi_psv2vb_k(nbig,nbod,m,xb,vb,pvb)
!     Save the current coordinates and velocities
      call mco_iden(time,jcen,nbod,nbig,h0,m,xb,vb,x0,v0,ngf,ngflag,opt)
!     Advancing the Keplerian Hamiltonian
      do j=2,nbod
        iflag=0
        call drift_one(gm(j),xb(1,j),xb(2,j),xb(3,j),vb(1,j),vb(2,j),
     &       vb(3,j),h0,iflag)
      end do
!     Check whether any object separations were < R_CRIT whilst advancing H_K
      call dpi_snif(h0,2,nbod,nbig,x0,v0,xb,vb,rcrit,ce,nce,ice,jce)
!     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
            do i=1,3
              xb(i,j) = x0(i,j)
              vb(i,j) = v0(i,j)
            end do
          end if
        end do
        call dpi_bhkce(time,tstart,h0,hrec,tol,rmax,en(3),jcen,rcen,
     &               nbod,nbig,m,xb,vb,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)
      end if
!     Computing transformed pseudovelocities
      call dpi_vb2psv_k(nbig,nbod,m,xb,vb,pvb)
      if (nbod.gt.2) then
!       Advancing the Jump Hamiltonian
        call dpi_wbjump(nbig,nbod,hby2,m,xb,pvb)
!       Computing accelerations
        call dpi_wbfor(nbig,nbod,xb,a,m,rcrit)
!       Advancing the Interaction Hamiltonian
        do j=2,nbod
          do i=1,3
            pvb(i,j)=pvb(i,j)+hby2*a(i,j)
          end do
        end do
      end if
!     Updating inertial (input) coordinates
      call dpi_wb2h(nbig,nbod,m,xh,vh,xb,pvb)
      end subroutine dpi_wbstep
!******************************************************************************
!     Dynamical Plug-In for Mercury 6
!     Subroutine to computate of energy and angular momentum in the transformed
!     coordinate system for S-type binary systems
!     Author: Diego Turrini
!     Last modified: June 2009
!******************************************************************************
      subroutine dpi_en(nbod,nbig,m,xh,vh,energy,angmom)
      implicit none
!     Input/Output variables
      integer nbod,nbig
      real*8 m(nbod)
      real*8 xh(3,nbod),vh(3,nbod)
      real*8 energy,angmom
!     Local variables
      integer i,j,l
      real*8 xb(3,nbod),pvb(3,nbod)
      real*8 mden,mtot,rp,drp(3),amom(3)
      real*8 r(nbod),p2(nbod),pbsum(3),sp(3)
      real*8 Hkep,Hint,Hjump
      real*8 vmod,vsqr
      external vmod,vsqr
!     Computing S type variables
      call dpi_h2wb(nbig,nbod,m,xh,vh,xb,pvb)
!     Variables initialization
      Hkep=0.d0
      Hint=0.d0
      Hjump=0.d0
      mtot=m(1)
      if (nbig.gt.2) then
        do j=2,nbig-1
          mtot=mtot+m(j)
        end do
      end if
      mden=mtot
      mtot=mtot+m(nbig)
!     Computing Keplerian energy
      r(1)=0.d0
      p2(1)=0.d0!vsqr(pvb(1,1),pvb(2,1),pvb(3,1))
      do j=2,nbig
        r(j)=vmod(xb(1,j),xb(2,j),xb(3,j))
        p2(j)=vsqr(pvb(1,j),pvb(2,j),pvb(3,j))
      end do
      if (nbig.gt.2) then
        do j=2,nbig-1
          Hkep=Hkep+0.5d0*m(j)*p2(j)
          Hkep=Hkep-m(j)*m(1)/r(j)
        end do
      end if
      Hkep=Hkep+m(nbig)*mtot*p2(nbig)/(2.0d0*mden)
      Hkep=Hkep-m(nbig)*mden/r(nbig)
      if (nbig.gt.2) then
!       Computation of support variables
        do i=1,3
          sp(i)=0.d0
        end do
        do j=2,nbig-1
          do i=1,3
            sp(i)=sp(i)+m(j)*xb(i,j)
          end do
        end do
        do i=1,3
          sp(i)=sp(i)/mden
        end do
!       Computation of the Interaction part of the Hamiltonian
        Hint=Hint+m(nbig)*m(1)/r(nbig)
        do j=2,nbig-1
          Hint=Hint+m(nbig)*m(j)/r(nbig)
        end do
        Hint=Hint-m(nbig)*m(1)/vmod((xb(1,nbig)+sp(1)),
     &  (xb(2,nbig)+sp(2)),(xb(3,nbig)+sp(3)))
        do j=2,nbig-1
          Hint=Hint-m(nbig)*m(j)/vmod((xb(1,nbig)-xb(1,j)+sp(1)),
     &    (xb(2,nbig)-xb(2,j)+sp(2)),(xb(3,nbig)-xb(3,j)+sp(3)))
        end do
        if (nbig.gt.3) then
          do j=2,nbig-1
            do l=j+1,nbig-1
              drp(1)=xb(1,j)-xb(1,l)
              drp(2)=xb(2,j)-xb(2,l)
              drp(3)=xb(3,j)-xb(3,l)
              rp=vmod(drp(1),drp(2),drp(3))
              Hint=Hint-(m(j)*m(l)/rp)
            end do
          end do
        end if
!       Computation of the Jump part of the Hamiltonian
        do i=1,3
          pbsum(i)=0.d0
        end do
        do j=2,nbig-1
          do i=1,3
            pbsum(i)=pbsum(i)+pvb(i,j)*m(j)
          end do
        end do
        Hjump=vsqr(pbsum(1),pbsum(2),pbsum(3))/(2.d0*m(1))
      end if
!     Computing total energy (including contribution from central star)
      energy=Hkep+Hint+Hjump!+(0.5d0*mtot*p2(1))
!     Computing angular momentum
      amom(1)=mtot*(xb(2,1)*pvb(3,1)-xb(3,1)*pvb(2,1))
      amom(2)=mtot*(xb(3,1)*pvb(1,1)-xb(1,1)*pvb(3,1))
      amom(3)=mtot*(xb(1,1)*pvb(2,1)-xb(2,1)*pvb(1,1))
      do j=2,nbig
        amom(1)=amom(1)+m(j)*(xb(2,j)*pvb(3,j)-xb(3,j)*pvb(2,j))
        amom(2)=amom(2)+m(j)*(xb(3,j)*pvb(1,j)-xb(1,j)*pvb(3,j))
        amom(3)=amom(3)+m(j)*(xb(1,j)*pvb(2,j)-xb(2,j)*pvb(1,j))
      end do
      angmom=vmod(amom(1),amom(2),amom(3))
      end subroutine dpi_en
!******************************************************************************
!     Dynamical Plug-In for Mercury 6
!     Function to computate the squared modulus of a vector (i.e. its dot
!     product against itself)
!     Author: Diego Turrini
!     Last modified: May 2009
!******************************************************************************
      function vsqr(x,y,z)
      implicit none
!     Input/Output variables
      real*8 vsqr,x,y,z
!     Computing scalar product
      vsqr=(x*x+y*y+z*z)
      end function vsqr
!******************************************************************************
!     Dynamical Plug-In for Mercury 6
!     Function to computate of the modulus of a vector (i.e. its root square
!     dot product)
!     Author: Diego Turrini
!     Last modified: May 2009
!******************************************************************************
      function vmod(x,y,z)
      implicit none
!     Input/Output variables
      real*8 vmod,x,y,z
      intrinsic sqrt
!     Computing vectorial module
      vmod=sqrt(x*x+y*y+z*z)
      end function vmod
!******************************************************************************
!     Dynamical Plug-In for Mercury 6
!     Function to compute the cubic power of a scalar variable
!     Author: Diego Turrini
!     Last modified: May 2009
!******************************************************************************
      function cube(x)
      implicit none
!     Input/Output variables
      real*8 cube,x
!     Computing cube of scalar x
      cube=x*x*x
      end function cube
!******************************************************************************
!     Dynamical Plug-In for Mercury 6
!     Subroutine for the numerical integration of the Keplerian part of the
!     Hamiltonian for S-type binary stars systems during close encounters
!     (derived from MDT_HKCE subroutine of Mercury 6.2)
!     Adapted by: Diego Turrini
!     Last modified: July 2009
!******************************************************************************
!
!     Author: John E. Chambers (subroutine MDT_HKCE in Mercury 6.2)
!
!     Integrates NBOD bodies (of which NBIG are Big) for one timestep H under
!     the Hamiltonian H_K, including close-encounter terms.
!
!******************************************************************************
      subroutine dpi_bhkce(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)
      implicit none
      include 'mercury.inc'
!     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)
!     Local
      integer iback(NMAX),indexs(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)
!     N.B.: Don't set nclo to zero!!
      nbs = 1
      nbsbig = 0
      mbs(1) = m(1)
      do k=1,3
        sbs(k,1) = s(k,1)
      end do
!     P-type binary algorithm still to be implemented, stop integration
      if (algor.eq.11) then
!         if (algor.eq.11) mbs(1) = m(1) + m(2)
        write(6,*) "P-type binary algorithm still to be implemented"
        write(6,*) "Stopping integration"
        stop
      end if
!     S-type binary algorithm does not allow collisions with binary companion
      if ((algor.eq.12).and.(ce(nbig).ne.0)) stop
!     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)
          do k=1,3
            xbs(k,nbs) = x(k,j)
            vbs(k,nbs) = v(k,j)
            sbs(k,nbs) = s(k,j)
          end do
          rcebs(nbs) = rce(j)
          rphybs(nbs) = rphy(j)
          statbs(nbs) = stat(j)
          rcritbs(nbs) = rcrit(j)
          idbs(nbs) = id(j)
          indexs(nbs) = j
          iback(j) = nbs
        end if
      end do
      do k = 1, nce
        ibs(k) = iback(ice(k))
        jbs(k) = iback(jce(k))
      end do
      tlocal = 0.d0
      hlocal = sign(hrec,h0)
!     Begin the Bulirsch-Stoer integration
  50  continue
        tmp0 = abs(h0) - abs(tlocal)
        hrec = hlocal
        if (abs(hlocal).gt.tmp0) hlocal = sign (tmp0, h0)
!       Save old coordinates and integrate
        call mco_iden (time,jcen,nbs,0,h0,mbs,xbs,vbs,x0,v0,ngf,ngflag,
     %    opt)
        call dpi_bs2 (time,hlocal,hdid,tol,jcen,nbs,nbsbig,mbs,xbs,vbs,
     %    sbs,rphybs,rcritbs,ngfbs,statbs,dtflag,ngflag,opt,nce,
     %    ibs,jbs)
        tlocal = tlocal + hdid
!       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)
!       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 dpi_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
!     If necessary, continue integrating objects undergoing close encounters
      if ((tlocal - h0)*h0.lt.0) goto 50
!     Return data for the close-encounter objects to global arrays
      do k = 2, nbs
        j = indexs(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) = 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