mercury6_2-dpi.for

c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c      MERCURY6_1.FOR    (ErikSoft   3 May 2002)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c Author: John E. Chambers
c
c Mercury is a general-purpose N-body integration package for problems in
c celestial mechanics.
c
c------------------------------------------------------------------------------
c This package contains some subroutines taken from the Swift integration 
c package by H.F.Levison and M.J.Duncan (1994) Icarus, vol 108, pp18.
c Routines taken from Swift have names beginning with `drift' or `orbel'.
c
c The standard symplectic (MVS) algorithm is described in J.Widsom and
c M.Holman (1991) Astronomical Journal, vol 102, pp1528.
c
c The hybrid symplectic algorithm is described in J.E.Chambers (1999)
c Monthly Notices of the RAS, vol 304, pp793.
c
c RADAU is described in E.Everhart (1985) in ``The Dynamics of Comets:
c Their Origin and Evolution'' p185-202, eds. A.Carusi & G.B.Valsecchi,
c pub. Reidel.
c
c The Bulirsch-Stoer algorithms are described in W.H.Press et al. (1992)
c ``Numerical Recipes in Fortran'', pub. Cambridge.
c------------------------------------------------------------------------------
c
c Variables:
c ---------
c  M      = mass (in solar masses)
c  XH     = coordinates (x,y,z) with respect to the central body (AU)
c  VH     = velocities (vx,vy,vz) with respect to the central body (AU/day)
c  S      = spin angular momentum (solar masses AU^2/day)
c  RHO    = physical density (g/cm^3)
c  RCEH   = close-encounter limit (Hill radii)
c  STAT   = status (0 => alive, <>0 => to be removed)
c  ID     = name of the object (8 characters)
c  CE     = close encounter status
c  NGF    = (1-3) cometary non-gravitational (jet) force parameters
c   "     =  (4)  beta parameter for radiation pressure and P-R drag
c  EPOCH  = epoch of orbit (days)
c  NBOD  = current number of bodies (INCLUDING the central object)
c  NBIG  =    "       "    " big bodies (ones that perturb everything else)
c  TIME  = current epoch (days)
c  TOUT  = time of next output evaluation
c  TDUMP = time of next data dump
c  TFUN  = time of next periodic effect (e.g. next check for ejections)
c  H     = current integration timestep (days)
c  EN(1) = initial energy of the system
c  " (2) = current    "    "  "    "
c  " (3) = energy change due to collisions, ejections etc.
c  AM(1,2,3) = as above but for angular momentum
c
c Integration Parameters :
c ----------------------
c  ALGOR = 1  ->  Mixed-variable symplectic
c          2  ->  Bulirsch-Stoer integrator
c          3  ->         "           "      (conservative systems only)
c          4  ->  RA15 `radau' integrator
c          10 ->  Hybrid MVS/BS (democratic-heliocentric coords)
c          11 ->  Close-binary hybrid (close-binary coords)
c          12 ->  Wide-binary hybrid (wide-binary coords)
c
c TSTART = epoch of first required output (days)
c TSTOP  =   "      final required output ( "  )
c DTOUT  = data output interval           ( "  )
c DTDUMP = data-dump interval             ( "  )
c DTFUN  = interval for other periodic effects (e.g. check for ejections)
c  H0    = initial integration timestep (days)
c  TOL   = Integrator tolerance parameter (approx. error per timestep)
c  RMAX  = heliocentric distance at which objects are considered ejected (AU)
c  RCEN  = radius of central body (AU)
c  JCEN(1,2,3) = J2,J4,J6 for central body (units of RCEN^i for Ji)
c
c Options:
c  OPT(1) = close-encounter option (0=stop after an encounter, 1=continue)
c  OPT(2) = collision option (0=no collisions, 1=merge, 2=merge+fragment)
c  OPT(3) = time style (0=days 1=Greg.date 2/3=days/years w/respect to start)
c  OPT(4) = o/p precision (1,2,3 = 4,9,15 significant figures)
c  OPT(5) = < Not used at present >
c  OPT(6) = < Not used at present >
c  OPT(7) = apply post-Newtonian correction? (0=no, 1=yes)
c  OPT(8) = apply user-defined force routine mfo_user? (0=no, 1=yes)
c
c File variables :
c --------------
c  OUTFILE  (1) = osculating coordinates/velocities and masses
c     "     (2) = close encounter details
c     "     (3) = information file
c  DUMPFILE (1) = Big-body data
c     "     (2) = Small-body data
c     "     (3) = integration parameters
c     "     (4) = restart file
c
c Flags :
c -----
c  NGFLAG = do any bodies experience non-grav. forces?
c                            ( 0 = no non-grav forces)
c                              1 = cometary jets only
c                              2 = radiation pressure/P-R drag only
c                              3 = both
c  OPFLAG = integration mode (-2 = synchronising epochs)
c                             -1 = integrating towards start epoch
c                              0 = main integration, normal output
c                              1 = main integration, full output
c
c------------------------------------------------------------------------------
c
      implicit none
      include 'mercury.inc'
c
      integer j,algor,nbod,nbig,opt(8),stat(NMAX),lmem(NMESS)
      integer opflag,ngflag,ndump,nfun
      real*8 m(NMAX),xh(3,NMAX),vh(3,NMAX),s(3,NMAX),rho(NMAX)
      real*8 rceh(NMAX),epoch(NMAX),ngf(4,NMAX),rmax,rcen,jcen(3)
      real*8 cefac,time,tstart,tstop,dtout,h0,tol,en(3),am(3)
      character*8 id(NMAX)
      character*80 outfile(3), dumpfile(4), mem(NMESS)
      external mdt_mvs, mdt_bs1, mdt_bs2, mdt_ra15, mdt_hy
      external mco_dh2h,mco_h2dh
      external mco_b2h,mco_h2b,mco_h2mvs,mco_mvs2h,mco_iden
! External subroutines and functions supplied by the DPI library
      external dpi_devbhc, dpi_en
c
      data opt/0,1,1,2,0,1,0,0/
c
c------------------------------------------------------------------------------
c
c Get initial conditions and integration parameters
      call mio_in (time,tstart,tstop,dtout,algor,h0,tol,rmax,rcen,jcen,
     %  en,am,cefac,ndump,nfun,nbod,nbig,m,xh,vh,s,rho,rceh,stat,id,
     %  epoch,ngf,opt,opflag,ngflag,outfile,dumpfile,lmem,mem)
!
! Algorithms implemented by the DPI library
!
      if (algor.eq.11) then
        write(6,*) "P-type binary algorithm not yet implemented."
        write(6,*) "Integration stopped."
        stop
      end if
      if (algor.eq.12) then
!       Compute (and overwrite) initial energy and angular momentum for
!       S-type binary systems using the transformed Hamiltonian
        call dpi_en(nbod,nbig,m,xh,vh,en(1),am(1))
      end if
!
c
c If this is a new integration, integrate all the objects to a common epoch.
      if (opflag.eq.-2) then
  20    open (23,file=outfile(3),status='old',access='append',err=20)
        write (23,'(/,a)') mem(55)(1:lmem(55))
        write (*,'(a)') mem(55)(1:lmem(55))
        call mxx_sync (time,tstart,h0,tol,jcen,nbod,nbig,m,xh,vh,s,rho,
     %    rceh,stat,id,epoch,ngf,opt,ngflag)
        write (23,'(/,a,/)') mem(56)(1:lmem(56))
        write (*,'(a)') mem(56)(1:lmem(56))
        opflag = -1
        close (23)
      end if
c
c Main integration
      if (algor.eq.1) call mal_hcon (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,mdt_mvs,mco_h2mvs,mco_mvs2h)
c
      if (algor.eq.9) call mal_hcon (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,mdt_mvs,mco_iden,mco_iden)
c
      if (algor.eq.2) call mal_hvar (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,mdt_bs1)
c
      if (algor.eq.3) call mal_hvar (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,mdt_bs2)
c
      if (algor.eq.4) call mal_hvar (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,mdt_ra15)
c
      if (algor.eq.10) call mal_hcon (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,mdt_hy,mco_h2dh,mco_dh2h)
!
! DPI algorithm for S-type binary systems
!
      if (algor.eq.12) call 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)
!                               
c
c Do a final data dump
      do j = 2, nbod
        epoch(j) = time
      end do
      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)
c
c Calculate and record the overall change in energy and ang. momentum
  50  open  (23, file=outfile(3), status='old', access='append',
     %  err=50)
      write (23,'(/,a)') mem(57)(1:lmem(57))
      call mxx_en (jcen,nbod,nbig,m,xh,vh,s,en(2),am(2))
c
      write (23,231) mem(58)(1:lmem(58)), 
     %  abs((en(2) + en(3) - en(1)) / en(1))
      write (23,232) mem(59)(1:lmem(59)), 
     %  abs((am(2) + am(3) - am(1)) / am(1))
      write (23,231) mem(60)(1:lmem(60)), abs(en(3) / en(1))
      write (23,232) mem(61)(1:lmem(61)), abs(am(3) / am(1))
      close (23)
      write (*,'(a)') mem(57)(1:lmem(57))
c
c------------------------------------------------------------------------------
c
 231  format (/,a,1p1e12.5)
 232  format (a,1p1e12.5)
      stop
      end
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c      MFO_USER.FOR    (ErikSoft   2 March 2001)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c Author: John E. Chambers
c
c Applies an arbitrary force, defined by the user.
c
c If using with the symplectic algorithm MAL_MVS, the force should be
c small compared with the force from the central object.
c If using with the conservative Bulirsch-Stoer algorithm MAL_BS2, the
c force should not be a function of the velocities.
c
c N.B. All coordinates and velocities must be with respect to central body
c ===
c------------------------------------------------------------------------------
c
      subroutine mfo_user (time,jcen,nbod,nbig,m,x,v,a)
c
      implicit none
      include 'mercury.inc'
c
c Input/Output
      integer nbod, nbig
      real*8 time,jcen(3),m(nbod),x(3,nbod),v(3,nbod),a(3,nbod)
c
c Local
      integer j
c
c------------------------------------------------------------------------------
c
      do j = 1, nbod
        a(1,j) = 0.d0
        a(2,j) = 0.d0
        a(3,j) = 0.d0
      end do
c
c------------------------------------------------------------------------------
c
      return
      end
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c      MAL_HVAR.FOR    (ErikSoft   4 March 2001)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c Author: John E. Chambers
c
c Does an integration using a variable-timestep integration algorithm. The
c particular integrator routine is ONESTEP and the algorithm must use
c coordinates with respect to the central body.
c
c N.B. This routine is also called by the synchronisation routine mxx_sync,
c ===  in which case OPFLAG = -2. Beware when making changes involving OPFLAG.
c
c------------------------------------------------------------------------------
c
      subroutine mal_hvar (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,onestep)
c
      implicit none
      include 'mercury.inc'
c
c Input/Output
      integer algor,nbod,nbig,stat(nbod),opt(8),opflag,ngflag,ndump,nfun
      integer lmem(NMESS)
      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)
c
c Local
      integer i,j,k,n,itmp,nhit,ihit(CMAX),jhit(CMAX),chit(CMAX)
      integer dtflag,ejflag,nowflag,stopflag,nstored,ce(NMAX)
      integer nclo,iclo(CMAX),jclo(CMAX),nce,ice(NMAX),jce(NMAX)
      real*8 tmp0,h,hdid,tout,tdump,tfun,tlog,tsmall,dtdump,dtfun
      real*8 thit(CMAX),dhit(CMAX),thit1,x0(3,NMAX),v0(3,NMAX)
      real*8 rce(NMAX),rphys(NMAX),rcrit(NMAX),a(NMAX)
      real*8 dclo(CMAX),tclo(CMAX),epoch(NMAX)
      real*8 ixvclo(6,CMAX),jxvclo(6,CMAX)
      external mfo_all,onestep
c
c------------------------------------------------------------------------------
c
c Initialize variables. DTFLAG = 0 implies first ever call to ONESTEP
      dtout  = abs(dtout)
      dtdump = abs(h0) * ndump
      dtfun  = abs(h0) * nfun
      dtflag = 0
      nstored = 0
      tsmall = h0 * 1.d-8
      h = h0
      do j = 2, nbod
        ce(j) = 0.d0
      end do
c
c Calculate close-encounter limits and physical radii for massive bodies
      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)
c
c 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
c
c------------------------------------------------------------------------------
c
c  MAIN  LOOP  STARTS  HERE
c
 100  continue
c
c Is it time for output ?
      if (abs(tout-time).lt.abs(tsmall).and.opflag.ge.-1) then
c
c Beware: the integration may change direction at this point!!!!
        if (opflag.eq.-1) dtflag = 0
c
c Output data for all bodies
        call mio_out (time,jcen,rcen,rmax,nbod,nbig,m,xh,vh,s,rho,
     %    stat,id,opt,opflag,algor,outfile(1))
        call mio_ce (time,tstart,rcen,rmax,nbod,nbig,m,stat,id,
     %    0,iclo,jclo,opt,stopflag,tclo,dclo,ixvclo,jxvclo,mem,lmem,
     %    outfile,nstored,0)
        tmp0 = tstop - tout
        tout = tout + sign( min( abs(tmp0), abs(dtout) ), tmp0 )
c
c Update the data dump files
        do j = 2, nbod
          epoch(j) = time
        end do
        call mio_dump (time,tstart,tstop,dtout,algor,h,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
c
c If integration has finished return to the main part of programme
      if (abs(tstop-time).le.abs(tsmall).and.opflag.ne.-1) return
c
c Set the timestep
      if (opflag.eq.-1) tmp0 = tstart - time
      if (opflag.eq.-2) tmp0 = tstop  - time
      if (opflag.ge.0)  tmp0 = tout   - time
      h = sign ( max( min( abs(tmp0), abs(h) ), tsmall), tmp0 )
c
c Save the current coordinates and velocities
      call mco_iden (time,jcen,nbod,nbig,h,m,xh,vh,x0,v0,ngf,ngflag,opt)
c
c Advance one timestep
      call onestep (time,h,hdid,tol,jcen,nbod,nbig,m,xh,vh,s,rphys,
     %  rcrit,ngf,stat,dtflag,ngflag,opt,nce,ice,jce,mfo_all)
      time = time + hdid
c
c Check if close encounters or collisions occurred
      nclo = 0
      call mce_stat (time,h,rcen,nbod,nbig,m,x0,v0,xh,vh,rce,rphys,
     %  nclo,iclo,jclo,dclo,tclo,ixvclo,jxvclo,nhit,ihit,jhit,
     %  chit,dhit,thit,thit1,nowflag,stat,outfile(3),mem,lmem)
c
c------------------------------------------------------------------------------
c
c  CLOSE  ENCOUNTERS
c
c If encounter minima occurred, output details and decide whether to stop
      if (nclo.gt.0.and.opflag.ge.-1) then
        itmp = 1
        if (nhit.ne.0) itmp = 0
        call mio_ce (time,tstart,rcen,rmax,nbod,nbig,m,stat,id,nclo,
     %    iclo,jclo,opt,stopflag,tclo,dclo,ixvclo,jxvclo,mem,lmem,
     %    outfile,nstored,itmp)
        if (stopflag.eq.1) return
      end if
c
c------------------------------------------------------------------------------
c
c  COLLISIONS
c
c If a collision occurred, output details and resolve the collision
      if (nhit.gt.0.and.opt(2).ne.0) then
        do k = 1, nhit
          if (chit(k).eq.1) then
            i = ihit(k)
            j = jhit(k)
            call mce_coll (thit(k),tstart,en(3),jcen,i,j,nbod,nbig,m,xh,
     %        vh,s,rphys,stat,id,opt,mem,lmem,outfile(3))
          end if
        end do
c
c Remove lost objects, reset flags and recompute Hill and physical radii
        call mxx_elim (nbod,nbig,m,xh,vh,s,rho,rceh,rcrit,ngf,stat,
     %    id,mem,lmem,outfile(3),itmp)
        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
c
c------------------------------------------------------------------------------
c
c  COLLISIONS  WITH  CENTRAL  BODY
c
c Check for collisions
      call mce_cent (time,hdid,rcen,jcen,2,nbod,nbig,m,x0,v0,xh,vh,nhit,
     %  jhit,thit,dhit,algor,ngf,ngflag)
c
c Resolve the collisions
      if (nhit.gt.0) then
        do k = 1, nhit
          i = 1
          j = jhit(k)
          call mce_coll (thit(k),tstart,en(3),jcen,i,j,nbod,nbig,m,xh,
     %      vh,s,rphys,stat,id,opt,mem,lmem,outfile(3))
        end do
c
c Remove lost objects, reset flags and recompute Hill and physical radii
        call mxx_elim (nbod,nbig,m,xh,vh,s,rho,rceh,rcrit,ngf,stat,
     %    id,mem,lmem,outfile(3),itmp)
        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),0)
      end if
c
c------------------------------------------------------------------------------
c
c  DATA  DUMP  AND  PROGRESS  REPORT
c
c 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
        call mio_ce (time,tstart,rcen,rmax,nbod,nbig,m,stat,id,
     %    0,iclo,jclo,opt,stopflag,tclo,dclo,ixvclo,jxvclo,mem,lmem,
     %    outfile,nstored,0)
        call mio_dump (time,tstart,tstop,dtout,algor,h,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
c
c Write a progress report to the log file
      if (abs(time-tlog).ge.abs(dtdump).and.opflag.ge.0) then
        call mxx_en (jcen,nbod,nbig,m,xh,vh,s,en(2),am(2))
        call mio_log (time,tstart,en,am,opt,mem,lmem)
        tlog = time
      end if
c
c------------------------------------------------------------------------------
c
c  CHECK  FOR  EJECTIONS  AND  DO  OTHER  PERIODIC  EFFECTS
c
      if (abs(time-tfun).ge.abs(dtfun).and.opflag.ge.-1) then
c
c Recompute close encounter limits, to allow for changes in Hill radii
        call mce_hill (nbod,m,xh,vh,rce,a)
        do j = 2, nbod
          rce(j) = rce(j) * rceh(j)
        end do
c
c Check for ejections
        call mxx_ejec (time,tstart,rmax,en,am,jcen,2,nbod,nbig,m,xh,vh,
     %    s,stat,id,opt,ejflag,outfile(3),mem,lmem)
c
c Remove lost objects, reset flags and recompute Hill and physical radii
        if (ejflag.ne.0) then
          call mxx_elim (nbod,nbig,m,xh,vh,s,rho,rceh,rcrit,ngf,stat,
     %      id,mem,lmem,outfile(3),itmp)
          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),0)
        end if
        tfun = time
      end if
c
c Go on to the next time step
      goto 100
c
c------------------------------------------------------------------------------
c
      end
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c      MAL_HCON.FOR    (ErikSoft   28 March 2001)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c Author: John E. Chambers
c
c Does an integration using an integrator with a constant stepsize H.
c Input and output to this routine use coordinates XH, and velocities VH,
c with respect to the central body, but the integration algorithm uses
c its own internal coordinates X, and velocities V.
c
c The programme uses the transformation routines COORD and BCOORD to change
c to and from the internal coordinates, respectively.
c
c------------------------------------------------------------------------------
c
      subroutine mal_hcon (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,onestep,
     %  coord,bcoord)
c
      implicit none
      include 'mercury.inc'
c
c Input/Output
      integer algor,nbod,nbig,stat(nbod),opt(8),opflag,ngflag
      integer lmem(NMESS),ndump,nfun
      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)
c
c Local
      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),xh0(3,NMAX),vh0(3,NMAX)
      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 onestep,coord,bcoord
c
c------------------------------------------------------------------------------
c
c 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.500001d0 * abs(h0)
c
c Calculate close-encounter limits and physical radii
      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)
c
c 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
c
c Convert to internal coordinates and velocities
      call coord (time,jcen,nbod,nbig,h0,m,xh,vh,x,v,ngf,ngflag,opt)
c
c------------------------------------------------------------------------------
c
c  MAIN  LOOP  STARTS  HERE
c
 100  continue
c
c Is it time for output ?
      if (abs(tout-time).le.hby2.and.opflag.ge.-1) then
c
c Beware: the integration may change direction at this point!!!!
        if (opflag.eq.-1.and.dtflag.ne.0) dtflag = 1
c
c Convert to heliocentric coordinates and output data for all bodies
        call bcoord (time,jcen,nbod,nbig,h0,m,x,v,xh,vh,ngf,ngflag,opt)
        call mio_out (time,jcen,rcen,rmax,nbod,nbig,m,xh,vh,s,rho,
     %    stat,id,opt,opflag,algor,outfile(1))
        call mio_ce (time,tstart,rcen,rmax,nbod,nbig,m,stat,id,
     %    0,iclo,jclo,opt,stopflag,tclo,dclo,ixvclo,jxvclo,mem,lmem,
     %    outfile,nstored,0)
        tmp0 = tstop - tout
        tout = tout + sign( min( abs(tmp0), abs(dtout) ), tmp0 )
c
c Update the data dump files
        do j = 2, nbod
          epoch(j) = time
        end do
        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
c
c If integration has finished, convert to heliocentric coords and return
      if (abs(tstop-time).le.hby2.and.opflag.ge.0) then
        call bcoord (time,jcen,nbod,nbig,h0,m,x,v,xh,vh,ngf,ngflag,opt)
        return
      end if
c
c 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)
c
c Save the current heliocentric coordinates and velocities
      if (algor.eq.1) then
        call mco_iden (time,jcen,nbod,nbig,h0,m,x,v,xh0,vh0,ngf,ngflag,
     %    opt)
      else
        call bcoord(time,jcen,nbod,nbig,h0,m,x,v,xh0,vh0,ngf,ngflag,opt)
      end if
      call onestep (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)
      time = time + h0
c
c------------------------------------------------------------------------------
c
c  CLOSE  ENCOUNTERS
c
c 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
        call mio_ce (time,tstart,rcen,rmax,nbod,nbig,m,stat,id,nclo,
     %    iclo,jclo,opt,stopflag,tclo,dclo,ixvclo,jxvclo,mem,lmem,
     %    outfile,nstored,itmp)
        if (stopflag.eq.1) return
      end if
c
c------------------------------------------------------------------------------
c
c  COLLISIONS
c
c If collisions occurred, output details and remove lost objects
      if (colflag.ne.0) then
c
c Reindex the surviving objects
        call bcoord (time,jcen,nbod,nbig,h0,m,x,v,xh,vh,ngf,ngflag,opt)
        call mxx_elim (nbod,nbig,m,xh,vh,s,rho,rceh,rcrit,ngf,stat,
     %    id,mem,lmem,outfile(3),itmp)
c
c Reset flags, and calculate new Hill radii and physical radii
        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)
        call coord (time,jcen,nbod,nbig,h0,m,xh,vh,x,v,ngf,ngflag,opt)
      end if
c
c------------------------------------------------------------------------------
c
c  COLLISIONS  WITH  CENTRAL  BODY
c
c Check for collisions with the central body
      if (algor.eq.1) then
        call mco_iden(time,jcen,nbod,nbig,h0,m,x,v,xh,vh,ngf,ngflag,opt)
      else
        call bcoord (time,jcen,nbod,nbig,h0,m,x,v,xh,vh,ngf,ngflag,opt)
      end if
      itmp = 2
      if (algor.eq.11.or.algor.eq.12) itmp = 3
      call mce_cent (time,h0,rcen,jcen,itmp,nbod,nbig,m,xh0,vh0,xh,vh,
     %  nhit,jhit,thit,dhit,algor,ngf,ngflag)
c
c If something hit the central body, restore the coords prior to this step
      if (nhit.gt.0) then
        call mco_iden (time,jcen,nbod,nbig,h0,m,xh0,vh0,xh,vh,ngf,
     %    ngflag,opt)
        time = time - h0
c
c Merge the object(s) with the central body
        do k = 1, nhit
          i = 1
          j = jhit(k)
          call mce_coll (thit(k),tstart,en(3),jcen,i,j,nbod,nbig,m,xh,
     %      vh,s,rphys,stat,id,opt,mem,lmem,outfile(3))
        end do
c
c Remove lost objects, reset flags and recompute Hill and physical radii
        call mxx_elim (nbod,nbig,m,xh,vh,s,rho,rceh,rcrit,ngf,stat,
     %    id,mem,lmem,outfile(3),itmp)
        if (opflag.ge.0) opflag = 1
        dtflag = 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),0)
        if (algor.eq.1) then
          call mco_iden (time,jcen,nbod,nbig,h0,m,xh,vh,x,v,ngf,ngflag,
     %      opt)
        else
          call coord (time,jcen,nbod,nbig,h0,m,xh,vh,x,v,ngf,ngflag,opt)
        end if
c
c Redo that integration time step
        goto 150
      end if
c
c------------------------------------------------------------------------------
c
c  DATA  DUMP  AND  PROGRESS  REPORT
c
c Convert to heliocentric coords and do the data dump
      if (abs(time-tdump).ge.abs(dtdump).and.opflag.ge.-1) then
        call bcoord (time,jcen,nbod,nbig,h0,m,x,v,xh,vh,ngf,ngflag,opt)
        do j = 2, nbod
          epoch(j) = time
        end do
        call mio_ce (time,tstart,rcen,rmax,nbod,nbig,m,stat,id,
     %    0,iclo,jclo,opt,stopflag,tclo,dclo,ixvclo,jxvclo,mem,lmem,
     %    outfile,nstored,0)
        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
c
c Convert to heliocentric coords and write a progress report to the log file
      if (abs(time-tlog).ge.abs(dtdump).and.opflag.ge.0) then
        call bcoord (time,jcen,nbod,nbig,h0,m,x,v,xh,vh,ngf,ngflag,opt)
        call mxx_en (jcen,nbod,nbig,m,xh,vh,s,en(2),am(2))
        call mio_log (time,tstart,en,am,opt,mem,lmem)
        tlog = time
      end if
c
c------------------------------------------------------------------------------
c
c  CHECK  FOR  EJECTIONS  AND  DO  OTHER  PERIODIC  EFFECTS
c
      if (abs(time-tfun).ge.abs(dtfun).and.opflag.ge.-1) then
        if (algor.eq.1) then
          call mco_iden (time,jcen,nbod,nbig,h0,m,x,v,xh,vh,ngf,ngflag,
     %      opt)
        else
          call bcoord(time,jcen,nbod,nbig,h0,m,x,v,xh,vh,ngf,ngflag,opt)
        end if
c
c Recompute close encounter limits, to allow for changes in Hill radii
        call mce_hill (nbod,m,xh,vh,rce,a)
        do j = 2, nbod
          rce(j) = rce(j) * rceh(j)
        end do
c
c Check for ejections
        itmp = 2
        if (algor.eq.11.or.algor.eq.12) itmp = 3
        call mxx_ejec (time,tstart,rmax,en,am,jcen,itmp,nbod,nbig,m,xh,
     %    vh,s,stat,id,opt,ejflag,outfile(3),mem,lmem)
c
c Remove ejected objects, reset flags, calculate new Hill and physical radii
        if (ejflag.ne.0) then
          call mxx_elim (nbod,nbig,m,xh,vh,s,rho,rceh,rcrit,ngf,stat,
     %      id,mem,lmem,outfile(3),itmp)
          if (opflag.ge.0) opflag = 1
          dtflag = 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),0)
          if (algor.eq.1) then
            call mco_iden (time,jcen,nbod,nbig,h0,m,xh,vh,x,v,ngf,
     %        ngflag,opt)
          else
            call coord (time,jcen,nbod,nbig,h0,m,xh,vh,x,v,ngf,ngflag,
     %        opt)
          end if
        end if
        tfun = time
      end if
c
c Go on to the next time step
      goto 100
c
c------------------------------------------------------------------------------
c
      end
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c    MCE_BOX.FOR    (ErikSoft   30 September 2000)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c Author: John E. Chambers
c
c Given initial and final coordinates and velocities, the routine returns
c the X and Y coordinates of a box bounding the motion in between the
c end points.
c
c If the X or Y velocity changes sign, the routine performs a quadratic
c interpolation to estimate the corresponding extreme value of X or Y.
c
c------------------------------------------------------------------------------
c
      subroutine mce_box (nbod,h,x0,v0,x1,v1,xmin,xmax,ymin,ymax)
c
      implicit none
      include 'mercury.inc'
c
c Input/Output
      integer nbod
      real*8 h,x0(3,nbod), x1(3,nbod), v0(3,nbod),v1(3,nbod)
      real*8   xmin(nbod), xmax(nbod), ymin(nbod),ymax(nbod)
c
c Local
      integer j
      real*8 temp
c
c------------------------------------------------------------------------------
c
      do j = 2, nbod
        xmin(j) = min (x0(1,j), x1(1,j))
        xmax(j) = max (x0(1,j), x1(1,j))
        ymin(j) = min (x0(2,j), x1(2,j))
        ymax(j) = max (x0(2,j), x1(2,j))
c
c If velocity changes sign, do an interpolation
        if ((v0(1,j).lt.0.and.v1(1,j).gt.0).or.
     %      (v0(1,j).gt.0.and.v1(1,j).lt.0)) then
          temp = (v0(1,j)*x1(1,j) - v1(1,j)*x0(1,j) 
     %            - .5d0*h*v0(1,j)*v1(1,j)) / (v0(1,j) - v1(1,j))
          xmin(j) = min (xmin(j),temp)
          xmax(j) = max (xmax(j),temp)
        end if
c
        if ((v0(2,j).lt.0.and.v1(2,j).gt.0).or.
     %      (v0(2,j).gt.0.and.v1(2,j).lt.0)) then
          temp = (v0(2,j)*x1(2,j) - v1(2,j)*x0(2,j) 
     %            - .5d0*h*v0(2,j)*v1(2,j)) / (v0(2,j) - v1(2,j))
          ymin(j) = min (ymin(j),temp)
          ymax(j) = max (ymax(j),temp)
        end if
      end do
c
c------------------------------------------------------------------------------
c
      return
      end
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c    MCE_CENT.FOR    (ErikSoft   4 March 2001)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c Author: John E. Chambers
c
c Checks all objects with index I >= I0, to see if they have had a collision
c with the central body in a time interval H, when given the initial and 
c final coordinates and velocities. The routine uses cubic interpolation
c to estimate the minimum separations.
c
c N.B. All coordinates & velocities must be with respect to the central body!!
c ===
c------------------------------------------------------------------------------
c
      subroutine mce_cent (time,h,rcen,jcen,i0,nbod,nbig,m,x0,v0,x1,v1,
     %  nhit,jhit,thit,dhit,algor,ngf,ngflag)
c
      implicit none
      include 'mercury.inc'
c
c Input/Output
      integer i0, nbod, nbig, nhit, jhit(CMAX), algor, ngflag
      real*8 time,h,rcen,jcen(3),m(nbod),x0(3,nbod),v0(3,nbod)
      real*8 x1(3,nbod),v1(3,nbod),thit(CMAX),dhit(CMAX),ngf(4,nbod)
c
c Local
      integer j
      real*8 rcen2,mco_acsh,a,q,u0,uhit,m0,mhit,mm,r0,mcen
      real*8 hx,hy,hz,h2,p,rr0,rr1,rv0,rv1,temp,e,v2
      real*8 xu0(3,NMAX),xu1(3,NMAX),vu0(3,NMAX),vu1(3,NMAX)
c
c------------------------------------------------------------------------------
c
      if (i0.le.0) i0 = 2
      nhit = 0
      rcen2 = rcen * rcen
      mcen = m(1)
c
c If using close-binary code, convert to coords with respect to the binary
c      if (algor.eq.11) then
c        mcen = m(1) + m(2)
c        call mco_h2ub (temp,jcen,nbod,nbig,h,m,x0,v0,xu0,vu0,ngf,ngflag)
c        call mco_h2ub (temp,jcen,nbod,nbig,h,m,x1,v1,xu1,vu1,ngf,ngflag)
c      end if
c
c Check for collisions with the central body
      do j = i0, nbod
        if (algor.eq.11) then
          rr0 = xu0(1,j)*xu0(1,j) + xu0(2,j)*xu0(2,j) +xu0(3,j)*xu0(3,j)
          rr1 = xu1(1,j)*xu1(1,j) + xu1(2,j)*xu1(2,j) +xu1(3,j)*xu1(3,j)
          rv0 = vu0(1,j)*xu0(1,j) + vu0(2,j)*xu0(2,j) +vu0(3,j)*xu0(3,j)
          rv1 = vu1(1,j)*xu1(1,j) + vu1(2,j)*xu1(2,j) +vu1(3,j)*xu1(3,j)
        else
          rr0 = x0(1,j)*x0(1,j) + x0(2,j)*x0(2,j) + x0(3,j)*x0(3,j)
          rr1 = x1(1,j)*x1(1,j) + x1(2,j)*x1(2,j) + x1(3,j)*x1(3,j)
          rv0 = v0(1,j)*x0(1,j) + v0(2,j)*x0(2,j) + v0(3,j)*x0(3,j)
          rv1 = v1(1,j)*x1(1,j) + v1(2,j)*x1(2,j) + v1(3,j)*x1(3,j)
        end if
c
c If inside the central body, or passing through pericentre, use 2-body approx.
        if ((rv0*h.le.0.and.rv1*h.ge.0).or.min(rr0,rr1).le.rcen2) then
          if (algor.eq.11) then
            hx = xu0(2,j) * vu0(3,j)  -  xu0(3,j) * vu0(2,j)
            hy = xu0(3,j) * vu0(1,j)  -  xu0(1,j) * vu0(3,j)
            hz = xu0(1,j) * vu0(2,j)  -  xu0(2,j) * vu0(1,j)
            v2 = vu0(1,j)*vu0(1,j) +vu0(2,j)*vu0(2,j) +vu0(3,j)*vu0(3,j)
          else
            hx = x0(2,j) * v0(3,j)  -  x0(3,j) * v0(2,j)
            hy = x0(3,j) * v0(1,j)  -  x0(1,j) * v0(3,j)
            hz = x0(1,j) * v0(2,j)  -  x0(2,j) * v0(1,j)
            v2 = v0(1,j)*v0(1,j) + v0(2,j)*v0(2,j) + v0(3,j)*v0(3,j)
          end if
          h2 = hx*hx + hy*hy + hz*hz
          p = h2 / (mcen + m(j))
          r0 = sqrt(rr0)
          temp = 1.d0 + p*(v2/(mcen + m(j)) - 2.d0/r0)
          e = sqrt( max(temp,0.d0) )
          q = p / (1.d0 + e)
c
c If the object hit the central body
          if (q.le.rcen) then
            nhit = nhit + 1
            jhit(nhit) = j
            dhit(nhit) = rcen
c
c Time of impact relative to the end of the timestep
            if (e.lt.1) then
              a = q / (1.d0 - e)
              uhit = sign (acos((1.d0 - rcen/a)/e), -h)
              u0   = sign (acos((1.d0 - r0/a  )/e), rv0)
              mhit = mod (uhit - e*sin(uhit) + PI, TWOPI) - PI
              m0   = mod (u0   - e*sin(u0)   + PI, TWOPI) - PI
            else
              a = q / (e - 1.d0)
              uhit = sign (mco_acsh((1.d0 - rcen/a)/e), -h)
              u0   = sign (mco_acsh((1.d0 - r0/a  )/e), rv0)
              mhit = mod (uhit - e*sinh(uhit) + PI, TWOPI) - PI
              m0   = mod (u0   - e*sinh(u0)   + PI, TWOPI) - PI
            end if
            mm = sqrt((mcen + m(j)) / (a*a*a))
            thit(nhit) = (mhit - m0) / mm + time
          end if
        end if
      end do
c
c------------------------------------------------------------------------------
c
      return
      end
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c      MCE_COLL.FOR    (ErikSoft   2 October 2000)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c Author: John E. Chambers
c