start3.f

Go to the documentation of this file.

      SUBROUTINE start3(ISUB)
*
*
*       Initialization & restart of triple system.
*       ------------------------------------------
*
      include 'common6.h'
      REAL*8  m
      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/clump/   bodys(ncmax,5),t0s(5),ts(5),steps(5),rmaxs(5),
     &                names(ncmax,5),isys(5)
*
*
*       Decide between new run, termination or collision (= 0, > 0, < 0).
      IF (isub.NE.0) THEN
          iterm = isub
          isub = iabs(isub)
          go to 100
      END IF
*
*       Define sequence with dominant binary component as reference body.
      name3(1) = jclose
*       Perturber index from routine IMPACT (ICOMP & JCOMP set in KSTERM).
      IF (body(icomp).LT.body(jcomp)) THEN
          k = icomp
          icomp = jcomp
          jcomp = k
      END IF
      name3(2) = jcomp
      name3(3) = icomp
*
      DO 2 k = 1,7
          cm(k) = 0.0d0
    2 CONTINUE
*
*       Transform to the local c.m. reference frame.
      DO 4 l = 1,3
          j = name3(l)
          m(l) = body(j)
          cm(7) = cm(7) + m(l)
          DO 3 k = 1,3
              x3(k,l) = x(k,j)
              xdot3(k,l) = xdot(k,j)
              cm(k) = cm(k) + m(l)*x3(k,l)
              cm(k+3) = cm(k+3) + m(l)*xdot3(k,l)
    3     CONTINUE
    4 CONTINUE
*
*       Set c.m. coordinates & velocities for triple system.
      DO 5 k = 1,6
          cm(k) = cm(k)/cm(7)
    5 CONTINUE
*
*       Specify initial conditions for three-body regularization.
      DO 8 l = 1,3
          DO 7 k = 1,3
              x3(k,l) = x3(k,l) - cm(k)
              xdot3(k,l) = xdot3(k,l) - cm(k+3)
    7     CONTINUE
    8 CONTINUE
*
*       Calculate internal energy and include in total subsystem energy.
      CALL trans3(0)
      esub = esub + 0.5d0*energy
*
*       Save global indices & particle attributes of subsystem.
      DO 10 l = 1,3
          j = name3(l)
          jlist(l) = j
          name3(l) = name(j)
          SIZE(l) = radius(j)
          ip(l) = kstar(j)
   10 CONTINUE
*
*       Make perturber list for potential energy correction.
      icm = jlist(1)
      go to 200
*
*       Obtain potential energy of resolved subsystem & perturbers.
   20 CALL nbpot(3,np,pot1)
*
*       Define subsystem indicator (ISYS = 1, 2, 3 for triple, quad, chain).
      isys(nsub+1) = 1
*
*       Form ghosts and initialize c.m. motion in ICOMP (= JLIST(1) = I).
      CALL subsys(3,cm)
*
*       Include interaction of subsystem c.m. & perturbers for net effect.
      CALL nbpot(1,np,pot2)
*
*       Form square of c.m. velocity correction due to differential force.
      vi2 = x0dot(1,icomp)**2 + x0dot(2,icomp)**2 + x0dot(3,icomp)**2
      corr = 1.0 + 2.0*(pot2 - pot1)/(cm(7)*vi2)
      IF (corr.LE.0.0d0) corr = 0.0
*
*       Modify c.m. velocity by net tidal energy correction.
      DO 30 k = 1,3
          x0dot(k,icomp) = sqrt(corr)*x0dot(k,icomp)
   30 CONTINUE
*
*       Remove ghosts from perturber neighbour lists.
      CALL nbrem(icm,3,np)
*
*       Set maximum integration interval equal to c.m. step.
      tmax = step(icomp)
*
*       Copy total energy and output & capture option for routine TRIPLE.
      cm(8) = be(3)
      kz15 = kz(15)
      kz27 = kz(27)
*
*       Assign new subsystem index and begin triple regularization.
      isub = nsub
      ntrip = ntrip + 1
      go to 180
*
*       Prepare KS regularization and direct integration of third body.
  100 imin = 1
      IF (r2.LT.r1) imin = 2
*
*       Specify global names of the KS candidates and least dominant body.
      nam1 = name3(3)
      nam2 = name3(imin)
      nam3 = name3(3-imin)
*
*       Identify current global indices by searching all single particles.
      i = 0
      DO 102 j = ifirst,n
          IF (name(j).EQ.nam1) icomp = j
          IF (name(j).EQ.nam2) jcomp = j
          IF (name(j).EQ.nam3) i = j
  102 CONTINUE
*
*       Identify global index for c.m. body.
      icm = i
      IF (body(icomp).GT.0.0d0) icm = icomp
      IF (body(jcomp).GT.0.0d0) icm = jcomp
*
*       Quantize the elapsed interval since last step.
      time2 = t0s(isub) + time3 - tprev
      dt8 = (tblock - tprev)/8.0d0
*
*       Adopt the nearest truncated step (at most 8 subdivisions).
      dt2 = time2
      IF (time2.GT.0.0d0) THEN
          CALL stepk(dt2,dtn2)
          dtn = nint(dtn2/dt8)*dt8
      ELSE
*       Choose negative step if pericentre time < TPREV (cf. iteration).
          dt2 = -dt2
          CALL stepk(dt2,dtn2)
          dtn = -nint(dtn2/dt8)*dt8
      END IF
*
*       Update time for new polynomial initializations (also for CMBODY).
      time = tprev + dtn
      time = min(tblock,time)
*
*       Predict current coordinates & velocities before termination.
      CALL xvpred(icm,0)
*
*       Copy c.m. coordinates & velocities.
      DO 105 k = 1,3
          cm(k) = x(k,icm)
          cm(k+3) = xdot(k,icm)
  105 CONTINUE
*
*       Re-determine the perturber list.
      go to 200
*
*       Obtain potential energy for the c.m. subsystem & JPERT(NP).
  110 jlist(1) = icm
      CALL nbpot(1,np,pot1)
*
*       Save global indices of three-body system.
      jlist(1) = i
      jlist(2) = icomp
      jlist(3) = jcomp
*
*       Set configuration pointers for escaper & KS candidates.
      jlist(4) = 3 - imin
      jlist(5) = 3
      jlist(6) = imin
*
*       Place new coordinates in the original locations.
      DO 120 l = 1,3
          j = jlist(l)
*       Compare global name & subsystem name to restore the mass.
          DO 112 k = 1,3
              IF (name(j).EQ.names(k,isub)) THEN
                  body(j) = bodys(k,isub)
                  t0(j) = time
              END IF
  112     CONTINUE
          ll = jlist(l+3)
          DO 115 k = 1,3
              x(k,j) = x3(k,ll) + cm(k)
  115     CONTINUE
  120 CONTINUE
*
*       Obtain potential energy of subsystem & perturbers at the end.
      CALL nbpot(3,np,pot2)
*
*       Form square of c.m. velocity correction due to differential force.
      vi2 = cm(4)**2 + cm(5)**2 + cm(6)**2
      corr = 1.0 + 2.0*(pot2 - pot1)/(cm(7)*vi2)
      IF (corr.LE.0.0d0) corr = 0.0
*
*       Modify c.m. velocity by net tidal energy correction.
      DO 122 k = 1,3
          cm(k+3) = sqrt(corr)*cm(k+3)
  122 CONTINUE
*
*       Transform to global velocities using corrected c.m. velocity.
      DO 130 l = 1,3
          j = jlist(l)
          ll = jlist(l+3)
          DO 125 k = 1,3
              xdot(k,j) = xdot3(k,ll) + cm(k+3)
              x0dot(k,j) = xdot(k,j)
  125     CONTINUE
  130 CONTINUE
*
*       Predict coordinates & velocities of perturbers to order FDOT.
      DO 140 l = 1,np
          j = jpert(l)
          CALL xvpred(j,0)
  140 CONTINUE
*
*       Update subsystem COMMON variables unless last or only case.
      IF (isub.LT.nsub) THEN
          DO 150 l = isub,nsub
              DO 145 k = 1,6
                  bodys(k,l) = bodys(k,l+1)
                  names(k,l) = names(k,l+1)
  145         CONTINUE
              t0s(l) = t0s(l+1)
              ts(l) = ts(l+1)
              steps(l) = steps(l+1)
              rmaxs(l) = rmaxs(l+1)
              isys(l) = isys(l+1)
  150     CONTINUE
      END IF
*
*       Reduce subsystem counter and subtract internal binding energy.
      nsub = nsub - 1
      esub = esub - 0.5d0*energy - ecoll3
*
*       Ensure ICOMP < JCOMP for KS treatment.
      k = icomp
      icomp = min(icomp,jcomp)
      jcomp = max(k,jcomp)
*
*       Select consistent neighbours for single body and new c.m.
      rs0 = rs(icm)
      CALL nblist(i,rs0)
      CALL nblist(icomp,rs0)
*
*       Add any other ghosts to perturber list for replacement of #ICM.
      DO 160 jsub = 1,nsub
          DO 155 l = 1,4
              IF (names(l,jsub).EQ.0) go to 155
              DO 154 j = 1,n
              IF (name(j).EQ.names(l,jsub).AND.body(j).LE.0.0d0) THEN
                  np = np + 1
                  jpert(np) = j
                  go to 155
              END IF
  154         CONTINUE
  155     CONTINUE
  160 CONTINUE
*
*       Replace ICM in relevant neighbour lists by all subsystem members.
      CALL nbrest(icm,3,np)
*
*       Check for stellar collision (only needs coordinates & velocities).
      IF (iterm.LT.0) THEN
          jlist(1) = icomp
          jlist(2) = jcomp
          jlist(3) = i
*       Initialize fictitious fourth body for general treatment.
          jlist(4) = 0
          go to 170
      END IF
*
*       Exclude the dominant interaction for c.m. approximation (large FDOT).
      IF (r3**2.GT.cmsep2*min(r1,r2)**2) THEN
          nnb = 1
      ELSE
          nnb = 3
      END IF
*
*       Set dominant F & FDOT on ICOMP & JCOMP for body #I in FPOLY2.
      CALL fclose(icomp,nnb)
      CALL fclose(jcomp,nnb)
*
*       Initialize force polynomials & time-steps for body #I.
      CALL fpoly1(i,i,0)
      CALL fpoly2(i,i,0)
*
*       Perform KS regularization of dominant components.
      CALL ksreg
*
*       Check minimum two-body distance.
      dmin3 = min(dmin3,rcoll)
*
*       Update net binary energy change.
      sbcoll = sbcoll + cm(9)
*
*       Update number of DIFSY calls, tidal dissipations & collision energy.
  170 nstept = nstept + nstep3
      ndiss = ndiss + ndiss3
      ecoll = ecoll + ecoll3
      e(10) = e(10) + ecoll3
*
*       Check for subsystem at last COMMON dump (no restart with NSUB > 0).
      IF (nsub.EQ.0.AND.kz(2).GE.1) THEN
          IF (time - tdump.LT.time3) THEN
              tdump = time
              CALL mydump(1,2)
          END IF
      END IF
*
*       Set phase indicator = -1 to ensure new time-step list in INTGRT.
  180 iphase = -1
*
      RETURN
*
*       Form the current perturber list.
  200 rp2 = rs(icm)**2
      np = 0
      nnb2 = list(1,icm) + 1
*
*       Loop over all single particles & c.m. but skip subsystem members.
      DO 210 l = 2,nnb2
          j = list(l,icm)
          rij2 = (x(1,j) - cm(1))**2 + (x(2,j) - cm(2))**2 +
     &                                 (x(3,j) - cm(3))**2
          IF (rij2.LT.rp2) THEN
              IF (j.EQ.i) go to 210
              IF (j.EQ.icomp.OR.j.EQ.jcomp) go to 210
              np = np + 1
              jpert(np) = j
          END IF
  210 CONTINUE
*
      IF (isub.EQ.0) go to 20
      go to 110
*
      END