tpert.f

Go to the documentation of this file.

      SUBROUTINE tpert(IPAIR,GA,DT)
*
*
*       Perturbation time scale.
*       ------------------------
*
      include 'common6.h'
*
*
*       Set c.m. index and initialize scalars.
      i = n + ipair
      fmax = 0.0
      dtin = 1.0e+20
      jclose = 0
      nnb1 = list(1,i) + 1
*
*       Check rare case of no neighbours.
      IF (nnb1.LE.1) THEN
          dt = 2.0*stepr(i)
          go to 20
      END IF
*
*       Find the most likely perturbers (first approach & maximum force).
      DO 10 l = 2,nnb1
          j = list(l,i)
          rij2 = 0.0
          rdot = 0.0
*
          DO 6 k = 1,3
              xrel = x(k,j) - x(k,i)
              vrel = xdot(k,j) - xdot(k,i)
              rij2 = rij2 + xrel**2
              rdot = rdot + xrel*vrel
    6     CONTINUE
*
          vr = rdot/rij2
          IF (vr.LT.dtin) THEN
              dtin = vr
*       Note DTIN is inverse travel time to include case of no RDOT < 0.
              rcrit2 = rij2
              jcrit = j
          END IF
          fij = (body(i) + body(j))/rij2
          IF (fij.GT.fmax) THEN
              fmax = fij
              rjmin2 = rij2
              jclose = j
          END IF
   10 CONTINUE
*
*       Form radial velocity of body with shortest approach time (if any).
      rcrit = sqrt(rcrit2)
      rdot = rcrit*abs(dtin)
      a1 = 2.0/(body(i)*ga)
      semi = -0.5*body(i)/h(ipair)
*
*       Use the actual apocentre for unperturbed travel time.
      ecc2 = (1.0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(body(i)*semi)
      ri = semi*(1.0 + sqrt(ecc2))
*
*       Estimate time interval to reach tidal perturbation of GA.
      dt = (rcrit - ri*(body(jcrit)*a1)**0.3333)/rdot
*
*       Compare the travel time based on acceleration only.
      dtmax = sqrt(2.0d0*abs(dt)*rdot*rcrit2/(body(i) + body(jcrit)))
      dt = min(dt,dtmax)
*
*       Skip dominant force test if there is only one critical body.
      IF (jcrit.NE.jclose) THEN
*       Form the return time of the dominant body and choose the minimum.
          dr = sqrt(rjmin2) - ri*(body(jclose)*a1)**0.3333
          dtmax = sqrt(2.0d0*abs(dr)/fmax)
          dt = min(dt,dtmax)
      END IF
*
*       Apply safety test in case background force dominates c.m. motion.
      dt = min(dt,4.0d0*step(i))
      dt = max(dt,0.0d0)
*
   20 RETURN
*
      END