qpmod3.f

Go to the documentation of this file.

      SUBROUTINE qpmod3(IM,ITERM)
*
*
*       Modification of AZ variables for tidal dissipation.
*       ---------------------------------------------------
*
      IMPLICIT REAL*8  (a-h,m,o-z)
      common/azreg/  time3,tmax,q(8),p(8),r1,r2,r3,energy,m(3),x3(3,3),
     &               xdot3(3,3),cm(10),c11,c12,c19,c20,c24,c25,
     &               nstep3,name3(3),kz15,kz27
      common/close/  rij(4,4),rcoll,qperi,SIZE(4),ecoll3,ip(4)
      common/azcoll/  rk(3),qk(8),pk(8),icall,icoll,ndiss3
      REAL*8  de(2),radius(2),p1(8)
      INTEGER is(2)
*
*
*       Obtain tidal energy loss (bodies #IM & 3 with separation QPERI).
      radius(1) = SIZE(im)
      radius(2) = SIZE(3)
      is(1) = ip(im)
      is(2) = ip(3)
*
      CALL tides(qperi,m(im),m(3),radius(1),radius(2),is,de)
*
*       Determine physical momenta of M(IM) & M(3-IM).
      DO 10 kcomp = 1,2
          k = 4*(kcomp - 1)
*
*       Form product of half Levi-Civita matrix & regularized momentum.
          p1(k+1) = q(k+1)*p(k+1) - q(k+2)*p(k+2) - q(k+3)*p(k+3) +
     &                                              q(k+4)*p(k+4)
          p1(k+2) = q(k+2)*p(k+1) + q(k+1)*p(k+2) - q(k+4)*p(k+3) -
     &                                              q(k+3)*p(k+4)
          p1(k+3) = q(k+3)*p(k+1) + q(k+4)*p(k+2) + q(k+1)*p(k+3) +
     &                                              q(k+2)*p(k+4)
          p1(k+4) = 0.0d0
*
*       Divide by the distance to obtain physical momentum (eqn (53)).
          ri = q(k+1)**2 + q(k+2)**2 + q(k+3)**2 + q(k+4)**2
          p1(k+1) = 0.5d0*p1(k+1)/ri
          p1(k+2) = 0.5d0*p1(k+2)/ri
          p1(k+3) = 0.5d0*p1(k+3)/ri
   10 CONTINUE
*
*       Define KS index for closest bodies and second interaction.
      k1 = 4*(im - 1)
      k2 = 4*(2 - im)
*
*       Evaluate kinetic energy term of least dominant interaction.
      p2 = 0.0d0
      DO 15 k = 1,3
          p2 = p2 + p1(k2+k)**2
*       Form relative momentum of dominant motion from c.m. condition.
          p1(k1+k) = 2.0d0*p1(k1+k) + p1(k2+k)
   15 CONTINUE
*
*       Set consistent third distance for energy relation.
      a21 = q(1)*q(1) - q(2)*q(2) - q(3)*q(3) + q(4)*q(4)
     &    - q(5)*q(5) + q(6)*q(6) + q(7)*q(7) - q(8)*q(8)
      a22 = q(1)*q(2) - q(3)*q(4) - q(5)*q(6) + q(7)*q(8)
      a23 = q(1)*q(3) + q(2)*q(4) - q(5)*q(7) - q(6)*q(8)
      a22 = a22 + a22
      a23 = a23 + a23
      r3 = sqrt(a21*a21 + a22*a22 + a23*a23)
*
*       Define radius, mass & reduced mass for the dominant bodies.
      rm = min(r1,r2)
      mb = m(im) + m(3)
      mu = m(im)*m(3)/mb
*
*       Obtain binding energy from total energy and perturbing function.
      mu2 = m(3-im)*mb/(m(3-im) + mb)
      vp = 0.5d0*p2/mu2 - m(1)*m(2)/r3 - m(3-im)*m(3)/max(r1,r2)
      h = (0.5d0*energy - vp)/mu
*
*       Set semi-major axis, eccentricity & pericentre (assume R' = 0).
      semi = -0.5d0*mb/h
      ecc = 1.0 - rm/semi
      peri = semi*(1.0d0 - ecc)
*
*       Determine new eccentricity from angular momentum conservation.
      dh = -(de(1) + de(2))/mu
      am0 = semi*(1.0d0 - ecc**2)
      ecc2 = ecc**2 + 2.0d0*am0*dh/mb
      IF (ecc2.GT.0.0d0) THEN
          ecc1 = sqrt(ecc2)
      ELSE
          ecc1 = 0.0
      END IF
*
*       Update binding energy and set new semi-major axis & pericentre.
      h1 = h + dh
      semi1 = -0.5d0*mb/h1
      peri1 = semi1*(1.0d0 - ecc1)
*
*       Form KS coordinate scaling factor from pericentre ratio.
      c1 = sqrt(peri1/peri)
*       Specify KS velocity scaling from angular momentum conservation.
      c2 = 1.0/c1
*
*       Modify the KS coordinates & physical momentum of dominant motion.
      ri = 0.0d0
      DO 20 k = 1,4
          q(k1+k) = c1*q(k1+k)
          ri = ri + q(k1+k)**2
          p1(k1+k) = c2**2*p1(k1+k)
          p1(k1+k) = 0.5d0*(p1(k1+k) - p1(k2+k))
   20 CONTINUE
*
*       Transform to new regularized momentum (eqn (50)).
      k = k1
      p(k+1) = 2.0d0*(+q(k+1)*p1(k+1) + q(k+2)*p1(k+2) +
     &                                  q(k+3)*p1(k+3))
      p(k+2) = 2.0d0*(-q(k+2)*p1(k+1) + q(k+1)*p1(k+2) +
     &                                  q(k+4)*p1(k+3))
      p(k+3) = 2.0d0*(-q(k+3)*p1(k+1) - q(k+4)*p1(k+2) +
     &                                  q(k+1)*p1(k+3))
      p(k+4) = 2.0d0*(+q(k+4)*p1(k+1) - q(k+3)*p1(k+2) +
     &                                  q(k+2)*p1(k+3))
*
*       Update the smallest distance (new pericentre has increased).
      IF (im.EQ.1) THEN
          r1 = ri
      ELSE
          r2 = ri
      END IF
*
*       Transform to physical variables and obtain consistent energy.
      CALL trans3(3)
      CALL trans3(0)
*
*       Update diagnostic variables (note DE > 0).
      ecoll3 = ecoll3 + (de(1) + de(2))
      ndiss3 = ndiss3 + 1
*
*       Perform stability test (ITERM < 0 denotes termination).
      CALL stabl3(iterm)
*
*       Print diagnostic if eccentricity > 0.99.
      IF (ecc.GT.0.99) THEN
          WRITE (6,30)  name3(im), name3(3), semi1, ecc, ecc1, h, qperi
   30     FORMAT (' NEW QPMOD3   NAM AF E0 EF H QP ',
     &                              2i5,1p,e10.2,0p,2f8.3,f9.1,1p,e10.2)
          CALL flush(6)
      END IF
*
      RETURN
*
      END