triple.f

Go to the documentation of this file.

      SUBROUTINE triple(ISUB)
*
*
*       Three-body regularization.
*       --------------------------
*
*       Method of Aarseth & Zare, Celestial Mechanics 10, 185.
*       ......................................................
*
*
      IMPLICIT REAL*8  (a-h,m,o-z)
      parameter(ncmax=10)
      common/azreg/  time3,tmax,q(8),p(8),r1,r2,r3,energy,m(3),x3(3,3),
     &               xdot3(3,3),cm(10),c11,c12,c19,c20,c24,c25,
     &               nstep3,name3(3),kz15,kz27
      common/close/  rij(4,4),rcoll,qperi,SIZE(4),ecoll3,ip(4)
      common/azcoll/  rk(3),qk(8),pk(8),icall,icoll,ndiss3
      common/bssave/  ep(4),dsc,facm,tfac,itfac,jc
      common/ebsave/  ebs
      common/azconf/  iconf(3)
      common/clump/   bodys(ncmax,5),t0s(5),ts(5),steps(5),rmaxs(5),
     &                names(ncmax,5),isys(5)
      REAL*8  y(17)
      SAVE
*
*
*       Main variables for AZ regularization
*       ************************************
*
*       ---------------------------------------------------------------------
*       CM(1-7) Coordinates, velocities & total mass of triple system.
*       CM(8)   Total energy of N-body system (copied in routine START3).
*       CM(9)   Binary energy change (copied to SBCOLL in START3).
*       CM(10)  Relative energy error (not used).
*       C11     Inverse mass factor for DERQP3 (also C12,C19,C20,C24,C25).
*       ENERGY  Twice the initial total energy.
*       ICALL   Pericentre indicator (activated if MIN(R1,R2) < RSTAR).
*       ICOLL   Collision or capture indicator (activated in routine DERQP3).
*       IP      Polytropic index (=1: n = 3/2; =2: n = 2; =3: n = 3).
*       KZ15    Triple option copied in START3 (full output if > 1).
*       KZ27    Tidal capture & collision option (copied in START3).
*       M       Particle mass (CM(7) = M(1) + M(2) + M(3)).
*       NAME3   Particle identity (initialized to global name).
*       NREG    Number of regularization switches.
*       NSTEP3  Number of DIFSY3 calls.
*       ICONF   Three-body configuration array (initialized to 1, 2, 3).
*       P       Regularized momenta.
*       Q       Regularized coordinates.
*       QPERI   Pericentre distance of closest two-body motion.
*       R1      Distance between M(1) and M(3).
*       R2      Distance between M(2) and M(3).
*       R3      Distance between M(1) and M(2) (not regularized).
*       RCOLL   Minimum two-body separation.
*       RCRIT   Escape test criterion (RGRAV or harmonic mean with RMAXS).
*       RGRAV   Gravitational radius (sum (M(I)*M(J))/ABS(0.5*ENERGY)).
*       RIJ     Minimum pairwise separations.
*       RMAXS   Maximum of unperturbed & initial size (R1 + R2).
*       RSTAR   Pericentre check distance (= 0.2*RGRAV).
*       SIZE    Individual stellar radius (collisions or tidal capture).
*       TIME3   Local physical time in scaled units.
*       TCR     Local crossing time ((sum M(I))**2.5/ABS(ENERGY)**1.5).
*       TOL     Tolerance for DIFSY3 (1.0E-10 is recommended).
*       TMAX    Maximum integration time (based on c.m. step).
*       X3      Particle coordinates (X3(3,I) is Z-component).
*       XDOT3   Velocity components (XDOT3(3,I) is Z-component).
*       ---------------------------------------------------------------------
*
*
*       Save termination indicator and check for restart.
      iterm = isub
      IF (isub.GT.0) THEN
*       Synchronize next time interval with c.m. step unless termination.
          tmax = time3 + steps(isub)
          IF (steps(isub).LE.0.0d0) dtau3 = 0.01*dtau3
*       Copy input array for DIFSY3 (just in case).
          go to 15
      END IF
*
*       Obtain initial conditions from N-body COMMON.
      CALL start3(isub)
*
*       Specify tolerance & DSC for DIFSY3 (relative energy error < TOL).
      tol = 1.0d-10
      dsc = 1.0
*
*       Initialize diagnostic & COMMON variables.
      r12min = 100.0
      r3min = 100.0
      rcoll = 100.0
      zmi1 = 100.0
      zmi2 = 100.0
      zmi3 = 100.0
      DO 10 j = 1,3
          iconf(j) = j
          DO 5 k = 1,3
              rij(j,k) = 100.0
    5     CONTINUE
   10 CONTINUE
      icall = 0
      icoll = 0
      next = 0
      ndiss3 = 0
      ecoll3 = 0.0d0
      jc = 0
      itfac = 0
      nstep3 = 0
      itry = 0
      nreg = 0
*
*       Initialize local time & regularized time.
      time3 = 0.0d0
      tau3 = 0.0d0
      y(17) = 0.0d0
*       Specify the number of first-order equations for the integrator.
      neq = 17
*
*       Evaluate the initial total energy (superseded; done in START3).
*     CALL TRANS3(0)
*
*       Transform to regularized variables.
      CALL trans3(1)
*
*       Define gravitational radius and pericentre check distance.
      rgrav = (m(1)*m(2) + m(1)*m(3) + m(2)*m(3))/(0.5d0*abs(energy))
      rstar = 0.2*rgrav
*       Introduce average mass factor for Bulirsch-Stoer integrator.
      facm = (m(1)*m(2) + m(1)*m(3) + m(2)*m(3))/3.0
*
*       Ensure that unperturbed boundary exceeds system size.
      r1 = q(1)**2 + q(2)**2 + q(3)**2 + q(4)**2
      r2 = q(5)**2 + q(6)**2 + q(7)**2 + q(8)**2
      IF (rmaxs(isub).LT.r1 + r2) rmaxs(isub) = r1 + r2
*
*       Modify escape distance criterion in case of hierarchical system.
      rcrit = max(rgrav,sqrt(rgrav*rmaxs(isub)))
*
*       Specify local crossing time (also meaningful if ENERGY > 0).
      tcr = (m(1) + m(2) + m(3))**2.5/abs(energy)**1.5
*
*       Define a nominal crossing time in near-parabolic cases.
      rmax = max(r1,r2)
      tstar = rmax*sqrt(rmax/cm(7))
*
*       Determine the smallest two-body time-scale from parabolic orbit.
      im = 1
      rm = r1
      IF (r2.LT.r1) THEN
          im = 2
          rm = r2
      END IF
      vp2 = 2.0*(m(im) + m(3))/rm
      tp = rm/sqrt(vp2)
      tstar = min(tp,tstar)
*
*       Set step for time transformation DT/DTAU = R1*R2/(R1 + R2)**0.5.
      tpr = r1*r2/sqrt(r1 + r2)
      dtau3 = min(tcr,tstar)*tol**0.1/tpr
*
*       Form initial binding energy of the binary.
      vrel2 = (xdot3(1,2) - xdot3(1,3))**2 +
     &        (xdot3(2,2) - xdot3(2,3))**2 +
     &        (xdot3(3,2) - xdot3(3,3))**2
      eb0 = (0.5d0*vrel2/(m(2) + m(3)) - 1.0/r2)*m(2)*m(3)
*
*       Initialize input array for the integrator.
   15 DO 20 k = 1,8
          y(k) = q(k)
          y(k+8) = p(k)
   20 CONTINUE
      y(17) = time3
*
*       Advance the equations of motion by Bulirsch-Stoer integrator.
   30 CALL difsy3(neq,tol,dtau3,tau3,y)
*
*       Copy regularized coordinates, momenta & time to COMMON variables.
      DO 40 k = 1,8
          q(k) = y(k)
          p(k) = y(k+8)
   40 CONTINUE
      time0 = time3
      time3 = y(17)
*
*       Update relative distances (NB! Not quite latest value).
      r1 = q(1)**2 + q(2)**2 + q(3)**2 + q(4)**2
      r2 = q(5)**2 + q(6)**2 + q(7)**2 + q(8)**2
*
*       Check minimum two-body separations and increase step counter.
      r3min = min(r3min,r3)
      rm = min(r1,r2)
      r12min = min(r12min,rm)
      nstep3 = nstep3 + 1
*
*       Check minimum pairwise separations (also in DERQP3).
      rk(1) = r1
      rk(2) = r2
      rk(3) = r3
*
*       Consider pairs 1-2, 1-3 & 2-3 (initial names = 1, 2, 3).
      DO 44 k = 1,3
          DO 42 l = k+1,3
              i = iconf(k)
              j = iconf(l)
*       Use cyclic loop index (3,1,2) for distances R3, R1 & R2.
              kk = k - 1
              IF (kk.EQ.0) kk = 3
              rij(i,j) = min(rij(i,j),rk(kk))
              rij(j,i) = min(rij(j,i),rk(kk))
   42     CONTINUE
   44 CONTINUE
*
*       Update smallest moment of inertia during normal forward integration.
      IF (time3.GT.time0.AND.jc.EQ.0) THEN
          zmi = rm**2
          zmi1 = zmi2
          zmi2 = zmi3
          zmi3 = zmi
          dtau0 = dtau3
      END IF
*
*       Switch on search indicator during iteration or just after pericentre.
      IF (icoll.LT.0) icall = 1
      IF (rm.LT.rstar.AND.nstep3.GT.next) THEN
          IF (zmi3.GT.zmi2.AND.zmi2.LT.zmi1) THEN
              icall = 1
          END IF
      END IF
*
*       Check for tidal dissipation or collision during last step.
      IF (icoll.LE.0) go to 48
*
*       Restore the minimum configuration.
      DO 45 k = 1,8
          q(k) = qk(k)
          p(k) = pk(k)
   45 CONTINUE
*
*       Set two-body separations.
      r1 = q(1)**2 + q(2)**2 + q(3)**2 + q(4)**2
      r2 = q(5)**2 + q(6)**2 + q(7)**2 + q(8)**2
*
*       Delay next search a few steps to avoid the same pericentre.
      next = nstep3 + 2 
*
*       Distinguish between tidal energy loss and collision (ICOLL holds IM).
      im = icoll
      icoll = 0
      IF (qperi.LT.4.0*max(SIZE(im),SIZE(3))) THEN
          IF (qperi.LT.0.75*(SIZE(im) + SIZE(3))) THEN
*
*       Transform to physical variables.
              CALL trans3(2)
*
*       Obtain global coordinates & velocities (ITERM < 0: termination).
              iterm = -1
              isub = -isub
              CALL start3(isub)
*
*       Combine internal energy and external c.m. kinetic energy.
              h3 = 0.5d0*energy + 0.5d0*cm(7)*(cm(4)**2 + cm(5)**2 +
     &                                                    cm(6)**2)
*
*       Evaluate the two-body energy for diagnostic purposes.
              CALL erel3(im,ebs,semi)
              dminc = min(rcoll,dminc)
*
*       Form composite body and begin KS regularization of new pair.
              CALL cmbody(h3,3)
              go to 100
          ELSE
*
*       Modify variables due to tidal effects and check stability parameter.
              CALL qpmod3(im,iterm)
*       Ensure positive step from pericentre and switch off search indicator.
              dtau3 = abs(dtau0)
              icall = 0
*       Update relative distances (NB! Not quite latest value).
      r1 = q(1)**2 + q(2)**2 + q(3)**2 + q(4)**2
      r2 = q(5)**2 + q(6)**2 + q(7)**2 + q(8)**2
              IF (iterm.LT.0) go to 90
*       Initialize input array and continue integration.
              go to 15
          END IF
      END IF
*
*       See whether switching of reference body is desirable.
   48 IF (r3.GT.rm) go to 70
*
*       Use a simple distance test to determine new reference body IMIN.
      imin = 1
      IF (r2.LT.r1) imin = 2
*
*       Transform to physical variables and rename the exchanged particles.
      CALL trans3(2)
*
      DO 50 k = 1,3
          temp1 = x3(k,3)
          temp2 = xdot3(k,3)
          x3(k,3) = x3(k,imin)
          xdot3(k,3) = xdot3(k,imin)
          x3(k,imin) = temp1
          xdot3(k,imin) = temp2
   50 CONTINUE
*
      temp1 = m(3)
      m(3) = m(imin)
      m(imin) = temp1
      temp2 = SIZE(3)
      SIZE(3) = SIZE(imin)
      SIZE(imin) = temp2
      name33 = name3(3)
      name3(3) = name3(imin)
      name3(imin) = name33
      i3 = iconf(3)
      iconf(3) = iconf(imin)
      iconf(imin) = i3
      i3 = ip(3)
      ip(3) = ip(imin)
      ip(imin) = i3
*
*       Transform back to regularized variables and initialize input array.
      CALL trans3(4)
      DO 60 k = 1,8
          y(k) = q(k)
          y(k+8) = p(k)
   60 CONTINUE
*
*       Update regularization counter at the end of switching procedure.
      nreg = nreg + 1
*
*       Set consistent relative distances and minimum separation.
      r1 = q(1)**2 + q(2)**2 + q(3)**2 + q(4)**2
      r2 = q(5)**2 + q(6)**2 + q(7)**2 + q(8)**2
      rm = min(r1,r2)
*
*       Terminate triple integration if R3 > specified perturber limit.
   70 IF (r3.GT.rmaxs(isub)) go to 90
      IF (iterm.LT.0) go to 74
*
*       Check for hierarchical stability in case of tidal capture.
      IF (ndiss3.GT.0) THEN
          iskip = 10
          an = int(nstep3/float(iskip)) - float(nstep3)/float(iskip)
*       Evaluate stability parameter every ISKIP step.
          IF (abs(an).LT.0.01) THEN
              CALL stabl3(iterm)
              IF (iterm.LT.0) go to 90
          END IF
      END IF
*
*       Check temporary termination time for return to routine INTGRT.
      IF (time3.GT.tmax) THEN
          IF (steps(isub).LE.0.0d0) go to 90
          go to 100
      END IF
*
*       See if the configuration permits testing of escape condition.
      IF (r1 + r2 + r3.LT.3.0*rcrit) go to 30
*
*       Obtain current physical variables.
      CALL trans3(2)
*
*       Identify index of second binary component & escaper.
      imin = 1
      IF (r2.LT.r1) imin = 2
      i = 3 - imin
*
*       Skip escape test if distant body is approaching the system c.m.
      ridot = x3(1,i)*xdot3(1,i) + x3(2,i)*xdot3(2,i) +
     &                             x3(3,i)*xdot3(3,i)
      IF (ridot.LT.0.0) go to 30
*
*       Set distance & radial velocity of body #I with respect to binary.
      mb = m(imin) + m(3)
      fac = cm(7)/mb
      ri = sqrt(x3(1,i)**2 + x3(2,i)**2 + x3(3,i)**2)
      ridot = fac*ridot/ri
      ri = fac*ri
*
*       Check the escape criterion due to Standish (Celes. Mech. 4, 44).
      ratio = rgrav/(mb*ri)
      vcrit2 = 2.0*cm(7)*(1.0/ri + m(3)*m(imin)*ratio**2/(ri - rgrav))
      IF (ridot**2.LT.vcrit2.OR.ri.LT.rgrav) go to 30
*
*       Define basic variables for termination of tidal capture event.
   74 IF (iterm.LT.0.AND.ndiss3.GT.0) THEN
          CALL trans3(2)
          im = 1
          IF (r2.LT.r1) im = 2
          i = 3 - im
          ri = r3
          rm = min(r1,r2)
          CALL erel3(im,ebs,semi)
          mb = m(imin) + m(3)
          fac = cm(7)/mb
      END IF
*
*       Evaluate orbital elements.
      vrel2 = 0.0d0
      rdot = 0.0d0
      vi2 = 0.0d0
      DO 75 k = 1,3
          rdot = rdot + 
     &               (x3(k,3) - x3(k,imin))*(xdot3(k,3) - xdot3(k,imin))
          vrel2 = vrel2 + (xdot3(k,3) - xdot3(k,imin))**2
          vi2 = vi2 + xdot3(k,i)**2
   75 CONTINUE
*
*       Form outer semi-major axis (not used).
      vi2 = vi2*fac**2
      semi1 = 2.0d0/ri - vi2/cm(7)
      semi1 = 1.0/semi1
*
*       Determine semi-major axis & eccentricity of inner binary.
      rb = rm
      semi = 2.0d0/rb - vrel2/mb
      semi = 1.0/semi
      e = sqrt((1.0d0 - rb/semi)**2 + rdot**2/(semi*mb))
*
*       Delay new regularization if binary is near small pericentre.
      itry = itry + 1
      IF (rb.LT.semi.AND.itry.LE.10) go to 30
*
*       Set binary energy ratio for triple system and whole N-body system.
      eb = -m(imin)*m(3)/(semi*energy)
      et = -m(imin)*m(3)/(2.0*semi*cm(8))
*       Form relative perturbation on inner binary due to body #I.
      gb = 2.0*m(i)/mb*(rb/ri)**3
      db = (eb*energy - eb0)/eb0
*
*       Print final configuration for significant energy increase.
      IF (db.GT.0.1) THEN
          WRITE (6,80)  name3(imin), name3(3), mb, semi, e, eb, gb, rb,
     &                  m(i), ri, et
   80     FORMAT (/,' TRIPLE BINARY ',2i5,'  MB =',f7.4,'  A =',1p,e8.1,
     &     '  E =',0p,f5.2,'  EB =',f5.2,'  GB =',1p,e8.1,'  RB =',e8.1,
     &               '  MI =',0p,f7.4,'  RI =',1p,e8.1,'  ET =',0p,f6.3)
      END IF
*
*       Postpone termination by one small step unless restart case.
      IF (steps(isub).GT.0.0d0) THEN
          steps(isub) = 0.0d0
          go to 100
      END IF
*
*       Transform to physical variables for termination.
   90 CALL trans3(3)
*
*       Identify second component and obtain osculating elements.
      imin = 1
      IF (r2.LT.r1) imin = 2
      rb = min(r1,r2)
      vrel2 = (xdot3(1,imin) - xdot3(1,3))**2 +
     &        (xdot3(2,imin) - xdot3(2,3))**2 +
     &        (xdot3(3,imin) - xdot3(3,3))**2
*
*       Form the semi-major axis & energy of closest pair.
      semi = 2.0d0/rb - vrel2/(m(imin) + m(3))
      semi = 1.0/semi
      eb = -0.5d0*m(imin)*m(3)/semi
*
*       Avoid new KS regularization inside half semi-major axis.
      itry = itry + 1
      IF (rb.LT.abs(semi).AND.itry.LE.10) THEN
          next = nstep3 + 10
          go to 15
      END IF
*
      IF (steps(isub).GT.0.0d0) THEN
          steps(isub) = 0.0d0
          go to 100
      END IF
*
*       Check optional print diagnostics of triple integration.
      IF (kz15.GT.1.OR.db.GT.0.1) THEN
          tc = time3/tcr
*       Print relative energy change of binary (including exchange).
          db = (eb - eb0)/eb0
          WRITE (6,95)  nreg, r12min, r3min, r3, rgrav, tc, nstep3, db,
     &                  name3(3-imin)
   95     FORMAT (/,' END TRIPLE    NREG =',i3,'  MIN(RB&R3) =',1p,e8.1,
     &                e9.1,'  R3 =',e8.1,'  RG =',e8.1,'  TC =',0p,f5.1,
     &                      '  STEPS =',i4,'  DB =',f6.2,'  NMESC =',i5)
      END IF
*
*       Form binary energy change (set in SBCOLL; routine START3).
      cm(9) = eb - eb0
*
*       Transform to global variables and begin new KS & single body #I.
      CALL start3(isub)
*
*       Activate termination index for routine INTGRT.
      iterm = -1
*
*       Update current time unless termination and set subsystem index.
  100 IF (iterm.GE.0) ts(isub) = t0s(isub) + time3
      isub = iterm
*
      RETURN
*
      END