ksperi.f

Go to the documentation of this file.

      SUBROUTINE ksperi(IPAIR)
*
*
*       Pericentre KS variables.
*       ------------------------
*
      include 'common6.h'
      common/slow0/  range,islow(10)
*
*
*       Save current time and initialize integration indicator.
      time0 = time
      itime = 0
      i1 = 2*ipair - 1
      icm = n + ipair
*
*       Re-calculate TDOT2 which may have misleading value.
    1 td2 = 0.0d0
      DO 5 k = 1,4
          td2 = td2 + u0(k,ipair)*udot(k,ipair)
    5 CONTINUE
      td2 = 2.0*td2
*
*       Define auxiliary quantities and obtain the eccentricity.
      semi = -0.5d0*body(icm)/h(ipair)
      zeta = 1.0 - r(ipair)/semi
      psi = td2/sqrt(body(icm))
      ecc = sqrt(zeta**2 + psi**2/semi)
*
*       Avoid nearly circular orbits (undefined pericentre).
      IF (ecc.LT.0.0001) go to 30
*
*       Distinguish between near-collision, elliptic & hyperbolic case.
      q = (psi/ecc)**2
      IF (zeta.GT.0.0d0.AND.abs(q/semi).LT.0.1) THEN
*       Expansion of THETA/(SIN(THETA)*COS(THETA)) - 1 for Q = SIN(THETA)**2.
          sn1 = 1.0
          s = 0.0d0
          q = q/semi
          DO 8 it = 1,10
              sn1 = q*sn1*float(2*it)/float(2*it + 1)
              s = s + sn1
    8     CONTINUE
          s = s + sn1*q*0.9d0/(1.0 - q)
          dt = (r(ipair)*zeta - psi**2 + zeta*s*semi)*psi/
     &                                          (ecc**2*sqrt(body(icm)))
      ELSE IF (semi.GT.0.0d0) THEN
*       Determine the eccentric anomaly with respect to pericentre (0,PI).
          theta = atan2(abs(psi)/sqrt(semi),zeta)
*       Obtain pericentre time interval from Kepler's equation.
          dt = semi*sqrt(semi/body(icm))*(theta - abs(psi)/sqrt(semi))
      ELSE IF (semi.LT.0.0d0) THEN
          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/body(icm))*(abs(psi)/sqrt(a0) - theta)
      END IF
*
*       Re-define current time (NB! not quantized; T0(I1) may be << TIME).
      time = time - dt
*
*       Integrate backwards for perturbed motion (reflection gives errors).
      IF (list(1,i1).GT.0.AND.itime.EQ.0.AND.dt.GT.step(i1)) THEN
          time = time0
          imod = kslow(ipair)
          zmod = float(islow(imod))
          iph = iphase
          iphase = -1
*
*       Integrate step by step if interval is too large (note IPHASE < 0).
   10     IF (dt.GT.step(i1)) THEN
              time = time - step(i1)
              dt = dt - step(i1)
              h0(ipair) = h(ipair)
              z = -0.5d0*h(ipair)*dtau(ipair)**2
              CALL stumpf(ipair,z)
              dtau(ipair) = -abs(dtau(ipair))
              CALL ksint(i1)
              dtu = dtau(ipair)
*       Use negative DTU and treat STEP as positive (not used elsewhere).
              step(i1) = ((one6*tdot3(ipair)*dtu + 0.5*tdot2(ipair))*dtu
     &                                                   + r(ipair))*dtu
              step(i1) = -zmod*step(i1)
              itime = itime + 1
              IF (itime.LT.200) go to 10
          END IF
*
          itime = itime + 1
          dtu = dt/(r(ipair)*zmod)
          dtu0 = dtau(ipair)
          iter = 0
*       Determine the regularized step by Newton-Raphson iteration (DT > 0).
   20     y0 = dt - zmod*((one6*tdot3(ipair)*dtu +
     &                             0.5*tdot2(ipair))*dtu + r(ipair))*dtu
          ypr = -((0.5*tdot3(ipair)*dtu + tdot2(ipair))*dtu + r(ipair))
          ypr = zmod*ypr
          dtu = dtu - y0/ypr
          dt1 = ((one6*tdot3(ipair)*dtu + 0.5*tdot2(ipair))*dtu +
     &                                                     r(ipair))*dtu
          dt1 = zmod*dt1
          iter = iter + 1
          IF (abs(dt - dt1).GT.1.0d-10*step(i1).AND.iter.LT.5) go to 20
*
*       Integrate back to pericentre using temporary indicator < 0 for exit.
          time = time - dt
          dtau(ipair) = -dtu
*       Re-initialize Stumpff functions.
          h0(ipair) = h(ipair)
          z = -0.5d0*h(ipair)*dtau(ipair)**2
          CALL stumpf(ipair,z)
          CALL ksint(i1)
*       Update energy and set positive step in case no KS initialization.
          h0(ipair) = h(ipair)
          dtau(ipair) = abs(dtu0)
          step(i1) = zmod*dt
          iphase = iph
*       Use reflection procedure to improve provisional pericentre.
          go to 1
*       Note: typically 2 iterations and final TDOT2 > 0.
      END IF
*
*       Specify transformation coefficients (Seppo Mikkola's procedure).
      IF (zeta.GE.0.0) THEN
          xc = sqrt(0.5d0 + 0.5d0*zeta/ecc)
          ys = psi/(ecc*xc*sqrt(body(icm)))
      ELSE
*       Employ well behaved expressions for R > A (SM 29/5/97).
          xc = 0.5*abs(psi)/(sqrt(semi)*ecc)/sqrt(0.5d0-0.5d0*zeta/ecc)
*       Avoid division by small XC near apocentre (ZETA < 0 only).
          ys = 2.0*sqrt(semi/body(icm)*(0.5d0 - 0.5d0*zeta/ecc))
          IF (psi.LT.0.0) ys = -ys
      END IF
*
      zz = body(icm)/(4.0*semi)
      r(ipair) = 0.0d0
      tdot2(ipair) = 0.0d0
*
*       Generate analytical solutions for U & UDOT using old U0 & UDOT.
      t0(i1) = time
      DO 25 k = 1,4
          u(k,ipair) = u0(k,ipair)*xc - udot(k,ipair)*ys
          udot(k,ipair) = u0(k,ipair)*ys*zz + udot(k,ipair)*xc
          u0(k,ipair) = u(k,ipair)
          r(ipair) = r(ipair) + u(k,ipair)**2
          tdot2(ipair) = tdot2(ipair) + 2.0d0*u(k,ipair)*udot(k,ipair)
   25 CONTINUE
*
*       Do not allow R' < 0 for repeated pericentre in routine KSINT.
      IF (tdot2(ipair).LT.0.0d0) THEN
          tdot2(ipair) = 0.0d0
      END IF
*
*       Predict c.m. coordinates & velocities.
      IF (abs(dt).LT.step(icm).AND.time.GT.0.0) THEN
          CALL xvpred(icm,0)
      END IF
*
   30 RETURN
*
      END