expand.f

Go to the documentation of this file.

      SUBROUTINE expand(IPAIR,SEMI0)
*
*
*       Expansion (contraction) of KS orbit.
*       ------------------------------------
*
      include 'common6.h'
*
*
*       Evaluate square regularized velocity.
      v20 = 0.0
      DO 10 k = 1,4
          v20 = v20 + udot(k,ipair)**2
   10 CONTINUE
*
*       Form KS coordinate & velocity scaling factors (general point is OK).
      i = n + ipair
      IF (semi0.GE.0.0d0) THEN
          semi = -0.5d0*body(i)/h(ipair)
      ELSE
          semi = semi0
      END IF
      c2 = sqrt(semi/semi0)
      v2 = 0.5*(body(i) + h(ipair)*r(ipair)*(semi/semi0))
      c1 = sqrt(v2/v20)
*
*       Re-scale KS variables to new energy (H < 0: constant eccentricity).
      r(ipair) = 0.0d0
*       Retain sign of radial velocity for unperturbed KS (apo or peri).
*     TDOT2(IPAIR) = 0.0D0
      DO 20 k = 1,4
          u(k,ipair) = c2*u(k,ipair)
          udot(k,ipair) = c1*udot(k,ipair)
          u0(k,ipair) = u(k,ipair)
          r(ipair) = r(ipair) + u(k,ipair)**2
*         TDOT2(IPAIR) = TDOT2(IPAIR) + 2.0*U(K,IPAIR)*UDOT(K,IPAIR)
   20 CONTINUE
*
*       Re-initialize KS polynomials for perturbed or hyperbolic motion.
      t0(2*ipair-1) = time
      IF (list(1,2*ipair-1).GT.0.OR.h(ipair).GT.0.0) THEN
          CALL resolv(ipair,1)
          CALL kspoly(ipair,1)
      ELSE
*       Determine new interval for unperturbed motion (>= TK).
          tk = twopi*semi*sqrt(abs(semi)/body(i))
          CALL tpert(ipair,gmin,dt)
*       Adopt c.m. step instead if integer argument exceeds 10**9.
          IF (dt.LT.2.0e+09*tk) THEN
              k = 1 + int(0.5d0*dt/tk)
*       Restrict Kepler period to c.m. step (case of very wide orbit).
              step(2*ipair-1) = min(float(k)*tk,step(i))
          ELSE
              step(2*ipair-1) = step(i)
          END IF
      END IF
*
      RETURN
*
      END