Nbody6
tpert.f
Go to the documentation of this file.
1  SUBROUTINE tpert(IPAIR,GA,DT)
2 *
3 *
4 * Perturbation time scale.
5 * ------------------------
6 *
7  include 'common6.h'
8 *
9 *
10 * Set c.m. index and initialize scalars.
11  i = n + ipair
12  fmax = 0.0
13  dtin = 1.0e+20
14  jclose = 0
15  nnb1 = list(1,i) + 1
16 *
17 * Check rare case of no neighbours.
18  IF (nnb1.LE.1) THEN
19  dt = 2.0*stepr(i)
20  go to 20
21  END IF
22 *
23 * Find the most likely perturbers (first approach & maximum force).
24  DO 10 l = 2,nnb1
25  j = list(l,i)
26  rij2 = 0.0
27  rdot = 0.0
28 *
29  DO 6 k = 1,3
30  xrel = x(k,j) - x(k,i)
31  vrel = xdot(k,j) - xdot(k,i)
32  rij2 = rij2 + xrel**2
33  rdot = rdot + xrel*vrel
34  6 CONTINUE
35 *
36  vr = rdot/rij2
37  IF (vr.LT.dtin) THEN
38  dtin = vr
39 * Note DTIN is inverse travel time to include case of no RDOT < 0.
40  rcrit2 = rij2
41  jcrit = j
42  END IF
43  fij = (body(i) + body(j))/rij2
44  IF (fij.GT.fmax) THEN
45  fmax = fij
46  rjmin2 = rij2
47  jclose = j
48  END IF
49  10 CONTINUE
50 *
51 * Form radial velocity of body with shortest approach time (if any).
52  rcrit = sqrt(rcrit2)
53  rdot = rcrit*abs(dtin)
54  a1 = 2.0/(body(i)*ga)
55  semi = -0.5*body(i)/h(ipair)
56 *
57 * Use the actual apocentre for unperturbed travel time.
58  ecc2 = (1.0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(body(i)*semi)
59  ri = semi*(1.0 + sqrt(ecc2))
60 *
61 * Estimate time interval to reach tidal perturbation of GA.
62  dt = (rcrit - ri*(body(jcrit)*a1)**0.3333)/rdot
63 *
64 * Compare the travel time based on acceleration only.
65  dtmax = sqrt(2.0d0*abs(dt)*rdot*rcrit2/(body(i) + body(jcrit)))
66  dt = min(dt,dtmax)
67 *
68 * Skip dominant force test if there is only one critical body.
69  IF (jcrit.NE.jclose) THEN
70 * Form the return time of the dominant body and choose the minimum.
71  dr = sqrt(rjmin2) - ri*(body(jclose)*a1)**0.3333
72  dtmax = sqrt(2.0d0*abs(dr)/fmax)
73  dt = min(dt,dtmax)
74  END IF
75 *
76 * Apply safety test in case background force dominates c.m. motion.
77  dt = min(dt,4.0d0*step(i))
78  dt = max(dt,0.0d0)
79 *
80  20 RETURN
81 *
82  END