dpi1_0.for

!******************************************************************************
!
!                                 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)