peri.f

Go to the documentation of this file.

      SUBROUTINE peri(U,UPR,TPR,M1,M2,Q)
*
*
*       Pericentre determination.
*       -------------------------
*
*       Routine for estimating the closest two-body approach based on the
*       KS variable U and its derivative UPR. The derivative of time, TPR,
*       is also used to convert to physical variables. This routine can be
*       used in any two-body KS, AZ and chain regularization which provides
*       a KS transformation U of the relative vector and its derivative
*       with respect to an arbitrary time variable.
*       NB!! Only reliable for the first function call with Bulirsch-Stoer.
*       Developed by Seppo Mikkola, May 1990.
*
*       U = KS coordinates of the relative position vector of bodies M1, M2.
*       UPR = derivatives of U used by equation of motion.
*       TPR = derivative of physical time used by equation of motion.
*       M1,M2 = masses of the bodies.
*       Q = osculating pericentre distance.
*
      IMPLICIT REAL*8  (a-h,m,o-z)
      REAL*8  u(4),upr(4)
*
*
*       Form the physical distance and square KS velocity.
      r = u(1)**2 + u(2)**2 + u(3)**2 + u(4)**2
      rd = upr(1)**2 + upr(2)**2 + upr(3)**2 + upr(4)**2
*
*       Evaluate the square of the KS cross-product.
      cross = 0.0d0
      DO 2 i = 1,3
          DO 1 j = i+1,4
              cross = cross + (u(i)*upr(j) - u(j)*upr(i))**2
    1     CONTINUE
    2 CONTINUE
*
*       Transform to physical velocity and include the mass factor.
      coef = (2.0d0*r/tpr)**2/(m1 + m2)
*
*       Form semi-latus rectum (note that P contains factor(s) of R).
      p = cross*coef
*
*       Set inverse semi-major axis (division by R is safe because of P).
      oa = (2.0d0 - coef*rd)/r
*
*       Obtain eccentricity & pericentre distance (avoid argument < 0).
      ecc2 = max(1.0d0 - p*oa,0.0d0)
      ecc = sqrt(ecc2)
      q = p/(1.0d0 + ecc)
*
      RETURN
*
      END