deform.f

Go to the documentation of this file.

      SUBROUTINE deform(IPAIR,ECC0,ECC)
*
*
*       Deformation of elliptic 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 at peri or apocentre.
      i = n + ipair
      semi = -0.5*body(i)/h(ipair)
      IF (r(ipair).LT.semi) THEN
          efac = (1.0 - ecc)/(1.0 - ecc0)
          rnew = semi*(1.0 - ecc)
      ELSE
          efac = (1.0 + ecc)/(1.0 + ecc0)
          rnew = semi*(1.0 + ecc)
          rnew = efac*r(ipair)
      END IF
      c1 = sqrt(efac)
      v2 = 0.5*(body(i) + h(ipair)*rnew)
      IF(v2.LE.0.d0)THEN
         c2 = 1.0d-06
      ELSE
         c2 = sqrt(v2/v20)
      ENDIF
*
*       Re-scale KS variables at constant energy to new eccentricity (ECC).
      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) = c1*u(k,ipair)
          udot(k,ipair) = c2*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 motion.
      t0(2*ipair-1) = time
      IF (list(1,2*ipair-1).GT.0) THEN
          CALL resolv(ipair,1)
          CALL kspoly(ipair,1)
      ELSE
*       Determine new interval for unperturbed motion (>= TK).
          tk = twopi*semi*sqrt(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) = float(k)*min(tk,step(i))
          ELSE
              step(2*ipair-1) = step(i)
          END IF
      END IF
*
      RETURN
*
      END