tperi.f

Go to the documentation of this file.

      SUBROUTINE tperi(SEMI,Q,UPR,MB,DT)
*
*
*       Pericentre time for KS motion.
*       ------------------------------
*
*       Routine for calculating the pericentre time
*       from the KS coordinates Q and derivatives UPR.
*       SEMI = semi-major axis of dominant motion.
*       Q = KS coordinates of relative position vector.
*       UPR = derivatives of Q in standard KS formulation.
*       MB = combined mass of the two bodies.
*       DT = pericentre time interval (DT < 0 before peri).
*       Algorithm due to Seppo Mikkola (1991).
*
      IMPLICIT REAL*8  (a-h,m,o-z)
      REAL*8  q(4),upr(4)
*
*
*       Set scalar KS distance, velocity and radial velocity.
      r = q(1)**2 + q(2)**2 + q(3)**2 + q(4)**2
      v2 = upr(1)**2 + upr(2)**2 + upr(3)**2 + upr(4)**2
      vr = 2.0d0*(q(1)*upr(1) + q(2)*upr(2) + q(3)*upr(3) + q(4)*upr(4))
*
*       Define auxiliary variables.
      psi = vr/sqrt(mb)
      zeta = 4.0d0*v2/mb - 1.0d0
*
*       Form semi-major axis & eccentricity.
*     ALPHA = 2.0D0*(1.0D0 - 2.0D0*V2/MB)/R
*     SEMI = 1.0D0/ALPHA
*
*       Employ regular value of semi-major axis determined by routine EREL.
      alpha = 1.0d0/semi
      ecc = sqrt(zeta*zeta + alpha*psi*psi)
*
*       Distinguish between near-collision, elliptic & hyperbolic case.
      q1 = (psi/ecc)**2
      IF (zeta.GT.0.0d0.AND.abs(q1/semi).LT.0.1) THEN
*       Expansion of THETA/(SIN(THETA)*COS(THETA)) - 1 for Q = SIN(THETA)**2.
          sn1 = 1.0
          s = 0.0d0
          q1 = q1/semi
          DO 1 it = 1,10
              sn1 = q1*sn1*float(2*it)/float(2*it + 1)
              s = s + sn1
    1     CONTINUE
          s = s + sn1*q1*0.9d0/(1.0 - q1)
          dt = (r*zeta - psi**2 + zeta*s*semi)*psi/(ecc**2*sqrt(mb))
      ELSE IF (semi.GT.0.0d0) THEN
*       Determine the eccentric anomaly with respect to pericentre (-PI,PI).
          theta = datan2(psi/sqrt(semi),zeta)
*       Obtain pericentre time interval from Kepler's equation.
          dt = semi*sqrt(semi/mb)*(theta - psi/sqrt(semi))
      ELSE IF (semi.LT.0.0d0) THEN
*       Hyperbolic case.
          a1 = psi/(ecc*sqrt(abs(semi)))
*       Use EXP(F) = SINH(F) + COSH(F) to obtain the eccentric anomaly THETA.
          a2 = abs(a1) + sqrt(a1**2 + 1.0d0)
          theta = log(a2)
          IF (a1.LT.0.0d0) theta = -theta
          a0 = abs(semi)
          dt = a0*sqrt(a0/mb)*(psi/sqrt(a0) - theta)
      END IF
*
      RETURN
*
      END