Nbody6
tperi.f
Go to the documentation of this file.
1  SUBROUTINE tperi(SEMI,Q,UPR,MB,DT)
2 *
3 *
4 * Pericentre time for KS motion.
5 * ------------------------------
6 *
7 * Routine for calculating the pericentre time
8 * from the KS coordinates Q and derivatives UPR.
9 * SEMI = semi-major axis of dominant motion.
10 * Q = KS coordinates of relative position vector.
11 * UPR = derivatives of Q in standard KS formulation.
12 * MB = combined mass of the two bodies.
13 * DT = pericentre time interval (DT < 0 before peri).
14 * Algorithm due to Seppo Mikkola (1991).
15 *
16  IMPLICIT REAL*8 (a-h,m,o-z)
17  REAL*8 q(4),upr(4)
18 *
19 *
20 * Set scalar KS distance, velocity and radial velocity.
21  r = q(1)**2 + q(2)**2 + q(3)**2 + q(4)**2
22  v2 = upr(1)**2 + upr(2)**2 + upr(3)**2 + upr(4)**2
23  vr = 2.0d0*(q(1)*upr(1) + q(2)*upr(2) + q(3)*upr(3) + q(4)*upr(4))
24 *
25 * Define auxiliary variables.
26  psi = vr/sqrt(mb)
27  zeta = 4.0d0*v2/mb - 1.0d0
28 *
29 * Form semi-major axis & eccentricity.
30 * ALPHA = 2.0D0*(1.0D0 - 2.0D0*V2/MB)/R
31 * SEMI = 1.0D0/ALPHA
32 *
33 * Employ regular value of semi-major axis determined by routine EREL.
34  alpha = 1.0d0/semi
35  ecc = sqrt(zeta*zeta + alpha*psi*psi)
36 *
37 * Distinguish between near-collision, elliptic & hyperbolic case.
38  q1 = (psi/ecc)**2
39  IF (zeta.GT.0.0d0.AND.abs(q1/semi).LT.0.1) THEN
40 * Expansion of THETA/(SIN(THETA)*COS(THETA)) - 1 for Q = SIN(THETA)**2.
41  sn1 = 1.0
42  s = 0.0d0
43  q1 = q1/semi
44  DO 1 it = 1,10
45  sn1 = q1*sn1*float(2*it)/float(2*it + 1)
46  s = s + sn1
47  1 CONTINUE
48  s = s + sn1*q1*0.9d0/(1.0 - q1)
49  dt = (r*zeta - psi**2 + zeta*s*semi)*psi/(ecc**2*sqrt(mb))
50  ELSE IF (semi.GT.0.0d0) THEN
51 * Determine the eccentric anomaly with respect to pericentre (-PI,PI).
52  theta = datan2(psi/sqrt(semi),zeta)
53 * Obtain pericentre time interval from Kepler's equation.
54  dt = semi*sqrt(semi/mb)*(theta - psi/sqrt(semi))
55  ELSE IF (semi.LT.0.0d0) THEN
56 * Hyperbolic case.
57  a1 = psi/(ecc*sqrt(abs(semi)))
58 * Use EXP(F) = SINH(F) + COSH(F) to obtain the eccentric anomaly THETA.
59  a2 = abs(a1) + sqrt(a1**2 + 1.0d0)
60  theta = log(a2)
61  IF (a1.LT.0.0d0) theta = -theta
62  a0 = abs(semi)
63  dt = a0*sqrt(a0/mb)*(psi/sqrt(a0) - theta)
64  END IF
65 *
66  RETURN
67 *
68  END