quad.f

Go to the documentation of this file.

      SUBROUTINE quad(ISUB)
*
*
*       Four-body regularization. 
*       -------------------------
*
*       Method of Mikkola & Aarseth, Celestial Mechanics 47, 375.
*       .........................................................
*
*
      IMPLICIT REAL*8  (a-h,o-z)
      parameter(ncmax=10)
      REAL*8  m,mij,mb,mb1,savex(12),savexd(12)
      LOGICAL  switch,gtype,gtype0
      common/creg/  m(4),x(12),xd(12),p(12),q(12),time4,energy,epsr2,
     &              xr(9),w(9),r(6),ta(6),mij(6),cm(10),rmax4,tmax,
     &              ds,tstep,eps,nstep4,name4(4),kz15,kz27,nreg,nfn
      common/tpr/   switch,gtype,gtype0
      common/config/  r2(4,4),i1,i2,i3,i4
      common/close/  rij(4,4),rcoll,qperi,SIZE(4),ecoll3,ip(4)
      common/ccoll/  qk(12),pk(12),icall,icoll,ndiss4
      common/bssave/  ep(4),dsc,facm,tfac,itfac,jc
      common/ksave/  k1,k2
      common/ebsave/  ebs
      common/clump/   bodys(ncmax,5),t0s(5),ts(5),steps(5),rmaxs(5),
     &                names(ncmax,5),isys(5)
      SAVE
*
*
*       Main variables for four-body regularization
*       *******************************************
*
*       -------------------------------------------------------------------
*       CM(1-7) Coordinates, velocities & total mass of four-body system.
*       CM(8)   Total energy of N-body system (copied in routine START4).
*       CM(9)   Binary energy change (copied to BBCOL in START4).
*       DS      Regularized time-step (rescaled after switching & new TPR).
*       ENERGY  Total energy.
*       EPS     Tolerance for DIFSY4 (1.0E-08 is recommended).
*       EPSR2   Distance criterion for stabilized time transformation.
*       ICALL   Pericentre indicator (activated if R(IMIN) < EPSR2**0.5).
*       ICOLL   Collision indicator (activated in routine DERQP4).
*       IND     Sorting index (R(IND(6)) is largest distance).
*       I1-I4   Indices of final configuration (I1 & I2 is closest binary).
*       KZ15    Diagnostic option (copied in START4; full output if > 1).
*       M       Particle mass (CM(7) = M(1) + M(2) + M(3) + M(4)).
*       NAME4   Particle identity (initialized to global name).
*       NFN     Number of function calls.
*       NREG    Number of regularization switches.
*       NSTEP4  Number of DIFSY4 calls.
*       P       Regularized momenta.
*       Q       Regularized coordinates.
*       QPERI   Minimum two-body separation.
*       R       Mutual distances (R(1), R(2), R(3) are regularized).
*       RCOLL   Minimum two-body separation.
*       RGRAV   Gravitational radius ((sum M(I)*M(J))/ABS(ENERGY)).
*       RIJ     Minimum pairwise separations.
*       RMAXS   Maximum size of unperturbed configuration.
*       RMAX4   Maximum of unperturbed & initial size (+ 20 per cent).
*       R2      Square mutual distances (R2(I1,I2) is closest pair).
*       TCR     Local crossing time ((sum M(I))**2.5/ABS(2*ENERGY)**1.5).
*       TIME4   Local physical time in scaled units.
*       TMAX    Maximum integration time (based on c.m. step).
*       TSTEP   Initial time-step for regularization (physical units).
*       X       Particle coordinates (X(1), X(2), X(3) is X, Y, Z).
*       XD      Velocity components.
*       -------------------------------------------------------------------
*
*
*       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 = time4 + steps(isub)
          IF (steps(isub).LE.0.0d0) ds = 0.01*ds
          go to 30
      END IF
*
*       Obtain initial conditions from N-body COMMON.
      CALL start4(isub)
*
*       Specify tolerance & DSC for DIFSY4 (relative energy error < TOL).
      eps = 1.0d-08
      eps = 1.0d-12
      dsc = 1.0
*
*       Initialize diagnostic & COMMON variables.
      time4 = 0.0d0
      rcoll = 100.0
      DO 10 j = 1,4
          DO 5 k = 1,4
              rij(j,k) = 100.0
    5     CONTINUE
   10 CONTINUE
      icall = 0
      icoll = 0
      ndiss4 = 0
      jc = 0
      itfac = 0
      nstep4 = 0
      iq = 0
      ecoll3 = 0.0d0
      nreg = 0
      nfn = 0
      itry = 0
      isave = 0
      ratio = 0.0
      switch = .false.
*
*       Evaluate the initial total energy (superseded; done in START4).
*     CALL NEWSYS(X,XD,M,4,ENERGY,GAM)
*
*       Initialize regularized coordinates, momenta & particle indices.
      CALL newreg
*
*       Define gravitational radius & average mass factor for routine DERQP4.
      rgrav = (mij(1) + mij(2) + mij(3) + mij(4) + mij(5) + mij(6))/
     &                                                      abs(energy)
      facm = rgrav*abs(energy)/6.0
*       Set small distance criterion for stabilized time transformation.
      epsr2 = (0.2*rgrav)**2
*
*       Ensure that unperturbed boundary exceeds system size (+ 20 per cent).
      rmax = max(sqrt(r2(i1,i3)),sqrt(r2(i2,i4)))
      rmax4 = max(rmaxs(isub),1.2*rmax)
*
*       Specify local crossing time.
      tcr = (m(1) + m(2) + m(3) + m(4))**2.5/abs(2.0d0*energy)**1.5
*
*       Define a nominal crossing time in near-parabolic cases.
      tstar = rmax*sqrt(rmax/cm(7))
*
*       Determine the smallest two-body time-scale from parabolic orbit.
      rm = sqrt(r2(i1,i2))
      vp2 = 2.0*(m(i1) + m(i2))/rm
      tp = rm/sqrt(vp2)
      tstar = min(tp,tstar)
*
*       Set initial time-step in physical units (also used by RCHAIN).
      tstep = min(tcr,tstar)*eps**0.1
*
*       Convert physical step to regularized time.
      ds = tstep/(r(1)*r(2)*r(3))
*
*       Form initial binary energies (scaled by internal energy).
      vrel2 = 0.0d0
      vrel21 = 0.0d0
      DO 25 k = 1,3
          j1 = 3*(i1-1) + k
          j2 = 3*(i2-1) + k
          j3 = 3*(i3-1) + k
          j4 = 3*(i4-1) + k
          vrel2 = vrel2 + (xd(j1) - xd(j2))**2
          vrel21 = vrel21 + (xd(j3) - xd(j4))**2
   25 CONTINUE
      eb0 = 2.0d0/sqrt(r2(i1,i2)) - vrel2/(m(i1) + m(i2))
      eb10 = 2.0d0/sqrt(r2(i3,i4)) - vrel21/(m(i3) + m(i4))
      eb0 = -0.5d0*m(i1)*m(i2)*eb0/energy
      eb10 = -0.5d0*m(i3)*m(i4)*eb10/energy
*
*       Perform regularized integration until termination is reached.
   30 CALL rchain(iq)
*
*       See whether temporary or actual termination (otherwise R > RMAX4).
      IF (time4.GT.tmax) THEN
         IF (steps(isub).GT.0.0d0) go to 100
      END IF
*
*       Transform to physical variables and order particle indices.
      switch = .false.
      CALL endreg
      CALL status(x,j1,j2,j3,j4)
*
*       Save the current configuration.
      DO 32 j = 1,12
          savex(j) = x(j)
          savexd(j) = xd(j)
   32 CONTINUE
      tsave = time4
      ds0 = ds
*
*       Check for collision during last step.
      IF (icoll.LE.0) go to 40
*
*       Restore the minimum configuration.
      DO 35 k = 1,12
          q(k) = qk(k)
          p(k) = pk(k)
   35 CONTINUE
*
*       Transform to physical variables.
      switch = .false.
      CALL endreg
*
*       Order particle indices of the collision configuration.
*     CALL STATUS(X,J1,J2,J3,J4)
*
*       Distinguish between tidal energy loss and collision (ICOLL holds IM).
      im = icoll
      icoll = 0
      IF (qperi.LT.4.0*max(SIZE(k1),SIZE(k2))) THEN
          IF (qperi.LT.0.75*(SIZE(k1) + SIZE(k2))) THEN
*
*       Obtain global coordinates & velocities (ITERM < 0 denotes collision).
              iterm = -1
              isub = -isub
              CALL start4(isub)
*
*       Combine internal energy and external c.m. kinetic energy.
              h4 = energy + 0.5d0*cm(7)*(cm(4)**2 + cm(5)**2 + cm(6)**2)
*
*       Determine dominant two-body energy from non-singular expressions.
              CALL erel4(im,ebs,semi)
              dminc = min(rcoll,dminc)
*
*       Form composite body and begin KS regularization of closest pair.
              CALL cmbody(h4,4)
              go to 100
          ELSE
*
*       Modify regularized variables and check stability (ITERM < 0).
              CALL qpmod4(im,iterm)
              switch = .true.
              CALL newreg
              icall = 0
              iq = iterm
              IF (iterm.EQ.-1) go to 100
              go to 30
          END IF
      END IF
*
*       Form output diagnostics (only hard binary printed if KZ15 < 2).
   40 vrel2 = 0.0d0
      vrel21 = 0.0d0
      rdot = 0.0d0
*
      DO 50 k = 1,3
          j1 = 3*(i1-1) + k
          j2 = 3*(i2-1) + k
          j3 = 3*(i3-1) + k
          j4 = 3*(i4-1) + k
          vrel2 = vrel2 + (xd(j1) - xd(j2))**2
          vrel21 = vrel21 + (xd(j3) - xd(j4))**2
          rdot = rdot + (x(j1) - x(j2))*(xd(j1) - xd(j2))
   50 CONTINUE
*
*       Evaluate orbital elements.
      rb = sqrt(r2(i1,i2))
      rb1 = sqrt(r2(i3,i4))
      r13 = sqrt(r2(i1,i3))
      r24 = sqrt(r2(i2,i4))
      mb = m(i1) + m(i2)
      semi = 2.0d0/rb - vrel2/mb
      semi = 1.0/semi
      ecc = sqrt((1.0d0 - rb/semi)**2 + rdot**2/(semi*mb))
      eb = -m(i1)*m(i2)/(2.0d0*semi*energy)
      et = -m(i1)*m(i2)/(2.0d0*semi*cm(8))
      mb1 = m(i3) + m(i4)
      semi1 = 2.0d0/rb1 - vrel21/mb1
      semi1 = 1.0/semi1
      eb1 = -m(i3)*m(i4)/(2.0d0*semi1*energy)
*       Scale binding energy of relative motion by initial value.
      e1 = (1.0 - eb - eb1)/eb0
*
*       Determine nearest and most distant binary perturber.
      IF (r13.LT.r24) THEN
          im = i3
          rm = r13
          imax = i4
          rmax = r24
      ELSE
          im = i4
          rm = r24
          imax = i3
          rmax = r13
      END IF
*
*       Estimate perturbations on smallest binary.
      gb = 2.0*mb1/mb*(rb/rm)**3
      IF (eb1.LT.0.0) gb = gb*m(im)/mb1
      g4 = 2.0*m(imax)/(mb + m(im))*(rm/rmax)**3
      db = (eb - eb0)/eb0
      IF (itry.LT.0) go to 70
*
*       Find the optimum configuration for two binaries or compact triple.
      itry = itry + 1
      IF (rm/semi.LT.ratio) go to 60
      IF (rb.LT.semi.OR.rb1.LT.semi1) go to 60
*       Avoid new regularization of any binary near small pericentre.
      ratio = rm/semi
      isave = itry
   60 IF (ratio.LT.5.0.AND.itry.LE.5.
     & or.ratio.LT.3.0.AND.itry.LE.10) THEN
          go to 30
      END IF
*
      IF (isave.EQ.itry) go to 70
      itry = -1
*       Restore the maximum triple configuration and obtain indices.
      DO 65 j = 1,12
          x(j) = savex(j)
          xd(j) = savexd(j)
   65 CONTINUE
      time4 = tsave
      CALL newreg
      ds = 0.5*ds0
      go to 40
*
*       Postpone temination by one small step unless restart case.
   70 IF (steps(isub).GT.0.0d0) THEN
          steps(isub) = 0.0d0
          go to 100
      END IF
*
*       Print final binary if relative energy increase exceeds 0.1.
      IF (db.GT.0.1) THEN 
          WRITE (6,80)  mb, semi, ecc, eb, gb, g4, rb1, eb1, e1, et
   80     FORMAT (/,' QUAD BINARY','  MB =',f7.4,'  A =',1p,e8.1,
     &     '  E =',0p,f5.2,'  EB =',f7.4,'  GB =',1p,e8.1,'  G4 =',e8.1,
     &     '  RB1 =',e8.1,'  EB1 =',0p,f5.2,'  E1 =',f5.2,'  ET =',f6.3)
      END IF
*
*       Check optional print diagnostics of four-body regularization.
      IF (kz15.GT.1) THEN
          tc = time4/tcr
          ec = energy/cm(8)
          WRITE (6,90)  i1, i2, i3, i4, rb, r13, r24, rgrav, tc, nstep4,
     &                  nreg, db, ec
   90     FORMAT (/,' END QUAD   ',4i3,'  RB =',1p,e8.1,'  R13 =',e8.1,
     &              '  R24 =',e8.1,'  RG =',e8.1,'  TC =',0p,f5.1,'  #',
     &                 i5,i4,'  DB =',f5.2,'  EC =',f6.3)
      END IF
*
*       Form binary energy change (set in BBCOLL; routine START4).
      cm(9) = (eb - eb0 - eb10)*energy
      IF (eb1.GT.0.0) cm(9) = cm(9) + eb1*energy
*
*       Transform to global variables and begin new KS (I1 & I2), I3 & I4.
      CALL start4(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) + time4
      isub = iterm
*
      RETURN
*
      END