derqp3.f

Go to the documentation of this file.

      SUBROUTINE derqp3(Q,P,QPR,PPR,TPR)
*
*
*       Equations of motion for AZ regularization.
*       ------------------------------------------
*
      IMPLICIT REAL*8  (a-h,m,o-z)
      common/azreg/  time3,tmax,v(8),w(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
*       Note COMMON locations of Q & P are replaced by dummy variables.
      common/close/  rij(4,4),rcoll,qperi,SIZE(4),ecoll3,ip(4)
      common/azcoll/  rk(3),qk(8),pk(8),icall,icoll,ndiss3
      common/bssave/  ep(4),dsc,facm,tfac,itfac,jc
      common/azconf/  iconf(3)
      REAL*8  q(8),p(8),qpr(8),ppr(8),s2(8),s5(8),s8(8),upr(8)
*
*
*       Form scalar distances & coefficients.
      r1 = q(1)*q(1) + q(2)*q(2) + q(3)*q(3) + q(4)*q(4)
      r2 = q(5)*q(5) + q(6)*q(6) + q(7)*q(7) + q(8)*q(8)
      c3 = q(1)*p(1) - q(2)*p(2) - q(3)*p(3) + q(4)*p(4)
      c4 = q(5)*p(5) - q(6)*p(6) - q(7)*p(7) + q(8)*p(8)
      c5 = q(2)*p(1) + q(1)*p(2) - q(4)*p(3) - q(3)*p(4)
      c6 = q(6)*p(5) + q(5)*p(6) - q(8)*p(7) - q(7)*p(8)
      c7 = q(3)*p(1) + q(4)*p(2) + q(1)*p(3) + q(2)*p(4)
      c8 = q(7)*p(5) + q(8)*p(6) + q(5)*p(7) + q(6)*p(8)
      c9 = p(1)*p(1) + p(2)*p(2) + p(3)*p(3) + p(4)*p(4)
      c10 = p(5)*p(5) + p(6)*p(6) + p(7)*p(7) + p(8)*p(8)
      c13 = c11*r2
      c14 = c12*r1
      c15 = c12*c10
      c16 = c11*c9
      c17 = r2*energy
      c18 = r1*energy
*       Note that twice the energy is saved in COMMON.
      c21 = 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)
      c22 = q(1)*q(2) - q(3)*q(4) - q(5)*q(6) + q(7)*q(8)
      c23 = q(1)*q(3) + q(2)*q(4) - q(5)*q(7) - q(6)*q(8)
      c22 = c22 + c22
      c23 = c23 + c23
      rr = c21*c21 + c22*c22 + c23*c23
      r3 = sqrt(rr)
      a = c25/r3
*
*       Set first derivative of the physical time (standard transformation).
      tpr = r1*r2
      b = a*tpr/rr
*
*       Combine vectorial components with relevant coefficients.
      s2(1) = q(1)*c4 + q(2)*c6 + q(3)*c8
      s2(2) =-q(2)*c4 + q(1)*c6 + q(4)*c8
      s2(3) =-q(3)*c4 - q(4)*c6 + q(1)*c8
      s2(4) = q(4)*c4 - q(3)*c6 + q(2)*c8
      s2(5) = q(5)*c3 + q(6)*c5 + q(7)*c7
      s2(6) =-q(6)*c3 + q(5)*c5 + q(8)*c7
      s2(7) =-q(7)*c3 - q(8)*c5 + q(5)*c7
      s2(8) = q(8)*c3 - q(7)*c5 + q(6)*c7
      s5(1) = p(1)*c4 + p(2)*c6 + p(3)*c8
      s5(2) =-p(2)*c4 + p(1)*c6 + p(4)*c8
      s5(3) =-p(3)*c4 - p(4)*c6 + p(1)*c8
      s5(4) = p(4)*c4 - p(3)*c6 + p(2)*c8
      s5(5) = p(5)*c3 + p(6)*c5 + p(7)*c7
      s5(6) =-p(6)*c3 + p(5)*c5 + p(8)*c7
      s5(7) =-p(7)*c3 - p(8)*c5 + p(5)*c7
      s5(8) = p(8)*c3 - p(7)*c5 + p(6)*c7
      s8(1) = q(1)*c21 + q(2)*c22 + q(3)*c23
      s8(2) =-q(2)*c21 + q(1)*c22 + q(4)*c23
      s8(3) =-q(3)*c21 - q(4)*c22 + q(1)*c23
      s8(4) = q(4)*c21 - q(3)*c22 + q(2)*c23
      s8(5) =-q(5)*c21 - q(6)*c22 - q(7)*c23
      s8(6) = q(6)*c21 - q(5)*c22 - q(8)*c23
      s8(7) = q(7)*c21 + q(8)*c22 - q(5)*c23
      s8(8) =-q(8)*c21 + q(7)*c22 - q(6)*c23
      c1 = c17 - c15 + c19 + a*r2
      c2 = c18 - c16 + c20 + a*r1
*
*       Form derivatives for standard equations of motion.
      DO 10 i = 1,4
          k = i + 4
          qpr(i) = c13*p(i) + c24*s2(i)
          qpr(k) = c14*p(k) + c24*s2(k)
          ppr(i) = c1*q(i) - c24*s5(i) - b*s8(i)
          ppr(k) = c2*q(k) - c24*s5(k) - b*s8(k)
   10 CONTINUE
*
*       Check the Bulirsch-Stoer tolerance scaling TFAC = TPR*U (ITFAC > 1).
      IF (itfac.GT.0) THEN
          tfac = facm*(r1 + r2 + tpr/r3)
          itfac = -1
      END IF
*
*       Evaluate the regularized Hamiltonian (denoted GAMMA in 1974 paper).
      gamma = 0.5d0*(-r1*c17 + r2*c11*c9 + r1*c12*c10 - c20*r2 - c19*r1
     &                            - a*tpr) + c24*(c3*c4 + c5*c6 + c7*c8)
      g2 = 1.0/sqrt(r1 + r2)
*       Time transformation of 1976 paper DT/DTAU = R1*R2/(R1 + R2)**0.5.
*
*       ------------------------------------------------------------------
*     U = M(1)*M(3)*R2 + M(2)*M(3)*R1 + 0.5D0*A*TPR
*     G2 = 1.0/U
*       Alternative transformation DT/DTAU = R1*R2/U allows explicit time.
*       ------------------------------------------------------------------
*
*       Set the explicit time derivative and check tolerance scaling.
      tpr = tpr*g2
      IF (itfac.NE.0) THEN
          tfac = tfac*g2
          itfac = 0
      END IF
*
*       Include stabilization terms due to modified time transformation.
      gamma = gamma*g2
*     S6 = C19 + A*R2
      DO 20 i = 1,8
          qpr(i) = g2*qpr(i)
          ppr(i) = g2*(ppr(i) + gamma*q(i))
*         PPR(I) = G2*(PPR(I) + GAMMA*(S6*Q(I) - B*S8(I)))
*         IF (I.EQ.4) S6 = C20 + A*R1
   20 CONTINUE
*
*       Check osculating pericentre of closest pair (first call only).
      IF (icall.EQ.0) go to 50
*
*       Find index & distance of closest pair and set corresponding KS index.
      IF (r1.LT.r2) THEN
          im = 1
          rm = r1
      ELSE
          im = 2
          rm = r2
      END IF    
      k = 1 + 4*(im - 1)
*
*       Obtain pericentre by Mikkola's algorithm (perturbation < 0.005).
      gi = 2.0*m(3-im)*(rm/r3)**3/(m(im) + m(3))
      IF (gi.LT.0.005) THEN
          CALL peri(q(k),qpr(k),tpr,m(im),m(3),qperi)
      ELSE
          qperi = rm
      END IF
*
*       Identify the close bodies and check mutual distances.
      i = iconf(im)
      j = iconf(3)
      rij(i,j) = min(rij(i,j),qperi)
      rij(j,i) = min(rij(j,i),qperi)
*
*       Ckeck minimum two-body distance and switch off indicator.
      rcoll = min(rcoll,qperi)
      icall = 0
*
*       Check for tidal two-body interaction or stellar collision.
      IF (qperi.LT.4.0*max(SIZE(im),SIZE(3))) THEN
*
*       Convert QPR to standard KS with T' = R and form radial velocity R'.
          rpr = 0.0d0
          DO 25 j = k,k+3
              upr(j) = qpr(j)*rm/tpr
              rpr = rpr + 2.0d0*q(j)*upr(j)
   25     CONTINUE
*
*       Save minimum configuration.
          DO 30 j = 1,8
              qk(j) = q(j)
              pk(j) = p(j)
   30     CONTINUE
*
*       Determine small semi-major axis from non-singular expressions.
          CALL erel3(im,eb,semi)
*
*       Obtain pericentre time interval from elements & Kepler's equation.
          mb = m(im) + m(3)
          CALL tperi(semi,q(k),upr(k),mb,tp)
*
*       Activate collision indicator & B-S step selector (first time).
          IF (icoll.EQ.0) THEN
              icoll = -1
              jc = 1
          ELSE
*
*       Check convergence: radial velocity < 1.0E-06 parabolic velocity.
              IF (abs(rpr).LT.1.0e-08*sqrt(2.0d0*mb*rm).OR.
     &           (abs(dsc).LT.1.0e-08)) THEN
*       Reset B-S step selector and set collision indicator to AZ pair index.
                  jc = 0
                  icoll = im
              END IF
          END IF
*
*       Set regularized step for DIFSY3 using accelerated iteration rate.
          IF (rm.GT.2.0*qperi) THEN
              dsc = 2.0*abs(tp)/tpr
          ELSE
              dsc = abs(tp)/tpr
          END IF       
*
*       Ensure negative step beyond pericentre (restore step at end).
          IF (jc.GT.0.AND.rpr.GT.0.0d0) THEN
              dsc = -dsc
          END IF
          IF (jc.EQ.0) dsc = 1.0
      END IF
*
   50 RETURN
*
      END