Nbody6
peri.f
Go to the documentation of this file.
1  SUBROUTINE peri(U,UPR,TPR,M1,M2,Q)
2 *
3 *
4 * Pericentre determination.
5 * -------------------------
6 *
7 * Routine for estimating the closest two-body approach based on the
8 * KS variable U and its derivative UPR. The derivative of time, TPR,
9 * is also used to convert to physical variables. This routine can be
10 * used in any two-body KS, AZ and chain regularization which provides
11 * a KS transformation U of the relative vector and its derivative
12 * with respect to an arbitrary time variable.
13 * NB!! Only reliable for the first function call with Bulirsch-Stoer.
14 * Developed by Seppo Mikkola, May 1990.
15 *
16 * U = KS coordinates of the relative position vector of bodies M1, M2.
17 * UPR = derivatives of U used by equation of motion.
18 * TPR = derivative of physical time used by equation of motion.
19 * M1,M2 = masses of the bodies.
20 * Q = osculating pericentre distance.
21 *
22  IMPLICIT REAL*8 (a-h,m,o-z)
23  REAL*8 u(4),upr(4)
24 *
25 *
26 * Form the physical distance and square KS velocity.
27  r = u(1)**2 + u(2)**2 + u(3)**2 + u(4)**2
28  rd = upr(1)**2 + upr(2)**2 + upr(3)**2 + upr(4)**2
29 *
30 * Evaluate the square of the KS cross-product.
31  cross = 0.0d0
32  DO 2 i = 1,3
33  DO 1 j = i+1,4
34  cross = cross + (u(i)*upr(j) - u(j)*upr(i))**2
35  1 CONTINUE
36  2 CONTINUE
37 *
38 * Transform to physical velocity and include the mass factor.
39  coef = (2.0d0*r/tpr)**2/(m1 + m2)
40 *
41 * Form semi-latus rectum (note that P contains factor(s) of R).
42  p = cross*coef
43 *
44 * Set inverse semi-major axis (division by R is safe because of P).
45  oa = (2.0d0 - coef*rd)/r
46 *
47 * Obtain eccentricity & pericentre distance (avoid argument < 0).
48  ecc2 = max(1.0d0 - p*oa,0.0d0)
49  ecc = sqrt(ecc2)
50  q = p/(1.0d0 + ecc)
51 *
52  RETURN
53 *
54  END