derqp4.f

Go to the documentation of this file.

      SUBROUTINE derqp4(P,Q,DP,DQ)
*
*
*       Equations of motion for chain regularization.
*       ---------------------------------------------
*
      IMPLICIT REAL*8  (a-h,o-z)
      REAL*8  m,mij,xnr(9),fnr(9),wi(9),at(3),aq(3),ap(9),fq(9),wp(9),
     &        p(12),q(12),dp(12),dq(13),upr(12),mb
      LOGICAL  switch,gtype,gtype0
      common/creg/  m(4),x(12),xd(12),pp(12),qq(12),time4,energy,epsr2,
     &              xr(9),w(9),r(6),ta(6),mij(6),cm(10),rmax4,tmax,
     &              ds,tstep,eps,nstep4,name4(4),kz15,kz27,nreg,nfn
      common/tpr/   switch,gtype,gtype0
      common/iconf/  i1,i2,i3,i4
      common/close/  rij(4,4),rcoll,qperi,SIZE(4),ecoll3,ip(4)
      common/ccoll/  qk(12),pk(12),icall,icoll,ndiss4
      common/bssave/  ep(4),dsc,facm,tfac,itfac,jc
      common/ksave/  k10,k20
*     EQUIVALENCE  (T11,TA(1)),(T22,TA(2)),(T33,TA(3)),(T12,TA(4)),
*    &             (T23,TA(5))
*
*       Copy old equivalence variables.
      t11 = ta(1)
      t22 = ta(2)
      t33 = ta(3)
      t12 = ta(4)
      t23 = ta(5)
*
      nfn = nfn + 1
      k = 0
      kp = 0
      DO 1 l = 1,3
          k1 = k + 1
          k2 = k + 2
          k3 = k + 3
          k4 = k + 4
          kp1 = kp + 1
          kp2 = kp + 2
          kp3 = kp + 3
*
*       Form physical vector & scalar distance.
          xr(kp1) = q(k1)**2 - q(k2)**2 - q(k3)**2 + q(k4)**2
          xr(kp2) = 2.d0*(q(k1)*q(k2) - q(k3)*q(k4))
          xr(kp3) = 2.d0*(q(k1)*q(k3) + q(k2)*q(k4))
          r(l) = q(k1)**2 + q(k2)**2 + q(k3)**2 + q(k4)**2
*       Set physical momentum times distance.
          wi(kp1) = q(k1)*p(k1) - q(k2)*p(k2) - q(k3)*p(k3) +
     &                                          q(k4)*p(k4)
          wi(kp2) = q(k2)*p(k1) + q(k1)*p(k2) - q(k4)*p(k3) -
     &                                          q(k3)*p(k4)
          wi(kp3) = q(k3)*p(k1) + q(k4)*p(k2) + q(k1)*p(k3) +
     &                                          q(k2)*p(k4)
          at(l) = ta(l)*(p(k1)**2 + p(k2)**2 + p(k3)**2 +
     &                                         p(k4)**2) - mij(l)
          k = k4
          kp = kp3
    1 CONTINUE
*
*       Obtain irregular vectors & distances.
      r(4) = 0.0d0
      r(5) = 0.0d0
      r(6) = 0.0d0
      DO 2 k = 1,3
          xnr(k  ) = xr(k  ) + xr(k+3)
          xnr(k+3) = xr(k+3) + xr(k+6)
          xnr(k+6) = xnr(k ) + xr(k+6)
          r(4) = r(4) + xnr(k  )**2
          r(5) = r(5) + xnr(k+3)**2
          r(6) = r(6) + xnr(k+6)**2
    2 CONTINUE
*
*       Evaluate irregular physical forces.
      unr = 0.0d0
      DO 5 i = 4,6
          dist = sqrt(r(i))
          fc = mij(i)/dist
          unr = unr + fc
          fc = fc/r(i)
          r(i) = dist
          ik = 3*(i - 4)
          DO 3 k = 1,3
              fnr(ik+k) = fc*xnr(ik+k)
    3     CONTINUE
    5 CONTINUE
      w1w2 = wi(1)*wi(4) + wi(2)*wi(5) + wi(3)*wi(6)
      w2w3 = wi(4)*wi(7) + wi(5)*wi(8) + wi(6)*wi(9)
      unre = unr + energy
      r12 = r(1)*r(2)
      r13 = r(1)*r(3)
      r23 = r(2)*r(3)
      aq(1) = at(2)*r(3) + at(3)*r(2) + t23*w2w3-r23*unre
      aq(2) = at(1)*r(3) + at(3)*r(1) - r13*unre
      aq(3) = at(2)*r(1) + at(1)*r(2) + t12*w1w2-r12*unre
      r123 = r12*r(3)
      ap(1) = 2.0d0*t11*r23
      ap(2) = 2.0d0*t22*r13
      ap(3) = 2.0d0*t33*r12
      wk12 = t12*r(3)
      wk23 = t23*r(1)
      DO 6 k = 1,3
          wp(k  ) = wk12*wi(k+3)
          wp(k+3) = wk12*wi(k  ) + wk23*wi(k+6)
          wp(k+6) = wk23*wi(k+3)
          fq(k  ) = fnr(k  ) + fnr(k+6)
          fq(k+3) = fnr(k  ) + fnr(k+3) + fnr(k+6)
          fq(k+6) = fnr(k+3) + fnr(k+6)
    6 CONTINUE
*
*       Form regularized derivatives.
      kq = 0
      l = 0
      DO 10 i = 1,3
          k1 = kq + 1
          k2 = kq + 2
          k3 = kq + 3
          k4 = kq + 4
          l1 = l + 1
          l2 = l + 2
          l3 = l + 3
*
          f1 = +fq(l1)*q(k1) + fq(l2)*q(k2) + fq(l3)*q(k3)
          f2 = -fq(l1)*q(k2) + fq(l2)*q(k1) + fq(l3)*q(k4)
          f3 = -fq(l1)*q(k3) - fq(l2)*q(k4) + fq(l3)*q(k1)
          f4 = +fq(l1)*q(k4) - fq(l2)*q(k3) + fq(l3)*q(k2)
*
          g1 = +wp(l1)*p(k1) + wp(l2)*p(k2) + wp(l3)*p(k3)
          g2 = -wp(l1)*p(k2) + wp(l2)*p(k1) + wp(l3)*p(k4)
          g3 = -wp(l1)*p(k3) - wp(l2)*p(k4) + wp(l3)*p(k1)
          g4 = +wp(l1)*p(k4) - wp(l2)*p(k3) + wp(l3)*p(k2)
*
          rrr = r123 + r123
          dp(k1) = -(2.d0*aq(i)*q(k1) + g1+rrr*f1)
          dp(k2) = -(2.d0*aq(i)*q(k2) + g2+rrr*f2)
          dp(k3) = -(2.d0*aq(i)*q(k3) + g3+rrr*f3)
          dp(k4) = -(2.d0*aq(i)*q(k4) + g4+rrr*f4)
*
          g1 = +wp(l1)*q(k1) + wp(l2)*q(k2) + wp(l3)*q(k3)
          g2 = -wp(l1)*q(k2) + wp(l2)*q(k1) + wp(l3)*q(k4)
          g3 = -wp(l1)*q(k3) - wp(l2)*q(k4) + wp(l3)*q(k1)
          g4 = +wp(l1)*q(k4) - wp(l2)*q(k3) + wp(l3)*q(k2)
*
          dq(k1) = ap(i)*p(k1) + g1
          dq(k2) = ap(i)*p(k2) + g2
          dq(k3) = ap(i)*p(k3) + g3
          dq(k4) = ap(i)*p(k4) + g4
*
          kq = kq + 4
          l = l + 3
   10 CONTINUE
*
*       Set the time derivative and check tolerance scaling (ITFAC > 1).
      dq(13) = r123
      IF (itfac.GT.0) THEN
          tfac = facm*(r(1)*r(2) + r(1)*r(3) + r(2)*r(3))
          itfac = -1
      END IF
*
*       Procedure for obtaining GAMMA (only after first call).
*     IF (IDIAG.EQ.1) THEN
*         GAMMA  =  AT(1)*R23 + AT(2)*R13 + AT(3)*R12 +
*    &              (R(3)*T12*W1W2 + T23*R(1)*W2W3) - R123*UNRE
*         IDIAG  =  0
*       NB! IDIAG must be in COMMON and set > 0 by RCHAIN.
*     END IF
*
      IF (gtype) THEN
*       Use modified time transformation for critical case.
          gamma = at(1)*r23 + at(2)*r13 + at(3)*r12 +
     &            r(3)*t12*w1w2 + t23*r(1)*w2w3 - r123*unre
          sigin = 1.d0/(r12 + r23 + r13)
          gsigin = 2.d0*gamma*sigin
          sumr = r(1) + r(2) + r(3)
          DO 15 i = 1,3
              si = (sumr - r(i))*gsigin
              kq = 4*(i - 1)
              DO 12 k = kq+1,kq+4
                  dq(k) = sigin*dq(k)
                  dp(k) = sigin*(dp(k) + si*q(k))
   12         CONTINUE
   15     CONTINUE
          dq(13) = dq(13)*sigin
          IF (itfac.NE.0) THEN
              tfac = tfac*sigin
              itfac = 0
          END IF
      END IF
*
*       Check osculating pericentre of closest pair (first call only).
      IF (icall.EQ.0) go to 50
*
*       Examine all close pair indices K1 & K2 in chain vector /ICONF/.
      qpmin = 100.0
      im0 = 1
      DO 20 im = 1,3
          IF (im.EQ.1) THEN
              k1 = i1
              k2 = i2
          ELSE IF (im.EQ.2) THEN
              k1 = i2
              k2 = i3
          ELSE
              k1 = i3
              k2 = i4
          END IF
*
*       Set distance of nearest perturber.
          IF (im.EQ.1.OR.im.EQ.3) THEN
              rp = r(2)
          ELSE
              rp = min(r(1),r(3))
          END IF
*
*       Obtain pericentre for small perturbations (ignore mass effect).
          gi = (r(im)/rp)**3
          IF (gi.LT.0.005) THEN
              k = 4*(im - 1) + 1
              CALL peri(q(k),dq(k),dq(13),m(k1),m(k2),qperi)
          ELSE
              qperi = r(im)
          END IF
*
*       Compare pericentre with previous mutual distances.
          rij(k1,k2) = min(rij(k1,k2),qperi)
          rij(k2,k1) = min(rij(k2,k1),qperi)
*
*       Save indices for smallest pericentre and switch off indicator.
          IF (im.EQ.1.OR.r(im).LE.r(im0)) THEN
              qpmin = qperi
              rp0 = rp
              k10 = k1
              k20 = k2
              im0 = im
              icall = 0
          END IF
   20 CONTINUE
*
*       Check smallest pericentre and reset corresponding indices.
      rcoll = min(rcoll,qpmin)
      qperi = qpmin
      im = im0
      k = 4*(im - 1) + 1
      k1 = k10
      k2 = k20
*
*       Save minimum configuration (for EREL4 & TRANS4).
      DO 25 j = 1,12
          qk(j) = q(j)
          pk(j) = p(j)
   25 CONTINUE
*
*       Check for tidal two-body interaction or stellar collision.
      IF (qpmin.LT.4.0*max(SIZE(k1),SIZE(k2))) THEN
*
*       Exit if collision distance is not the smallest.
          IF (r(im).GT.rp0) go to 50
*
*       Convert Q' to standard KS with T' = R and form radial velocity R'.
          rpr = 0.0d0
          DO 30 j = k,k+3
              upr(j) = dq(j)*r(im)/dq(13)
              rpr = rpr + 2.0d0*q(j)*upr(j)
   30     CONTINUE
*
*       Determine small semi-major axis from non-singular expressions.
          CALL erel4(im,eb,semi)
*
*       Obtain pericentre time interval from elements & Kepler's equation.
          mb = m(k1) + m(k2)
          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-06*sqrt(2.0d0*mb*r(im)).OR.
     &           (abs(dsc).LT.1.0e-06)) THEN
*       Reset B-S step selector and set collision indicator to CHAIN index.
                  jc = 0
                  icoll = im
              END IF
          END IF
*
*       Set regularized step for DIFSY4 using accelerated iteration rate.
          IF (r(im).GT.2.0*qpmin) THEN
              dsc = 2.0*abs(tp)/dq(13)
          ELSE
              dsc = abs(tp)/dq(13)
          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
*
*       Begin iteration if R(IM) is smallest or wide pair is past pericentre.
          IF (r(im).LT.rp0.OR.rpr.GT.0.0d0) THEN
              dsfac = 0.2*abs(ds/dsc)
              IF (r(im).LT.rp0) dsfac = 1.0
              ds = dsfac*dsc
          ELSE
*       Switch off indicators and carry on normally until next check.
              icoll = 0
              jc = 0
          END IF
      END IF
*
   50 RETURN
*
      END