derqp.f

Go to the documentation of this file.

      SUBROUTINE derqp(Q,CMX,ENERGY,P,CMV,CHTIME,DQ,DX,DE,DP,DV,DT)
*
*
*       Derivatives of chain variables.
*       -------------------------------
*
      include 'commonc.h'
      LOGICAL  kslow,kcoll
      REAL*8  ksch
      common/slow1/   tk2(0:nmx),ejump,ksch(nmx),kslow,kcoll
      common/cpert/  rgrav,gpert,ipert,npert
      common/calls/  tpr,tkpr,step,ider,icall,nfn,nreg,iter,imcirc
      common/ccoll2/  qk(nmx4),pk(nmx4),rij(nmx,nmx),SIZE(nmx),vstar1,
     &                ecoll1,rcoll,qperi,istar(nmx),icoll,isync,ndiss1
      common/bssave/  ep(4),dsc,facm,tfac,itfac,jc
      common/ebsave/  ebs
      common/ksave/  k10,k20
      REAL*8  q(nmx4),cmx(3),p(nmx4),cmv(3)
      REAL*8  dq(nmx4),dx(3),dp(nmx4),dv(3),fcm(3),upr(nmx4)
      REAL*8  w(nmx4),ak(nmx4),dk(nmx),fnc(nmx3),fxtnl(nmx3)
      REAL*8  tp(nmx4),tq(nmx4),uq(nmx4),faux(4),aux(-2:2),xaux(3)
      SAVE
*
*
*       Mikkola & Aarseth 1990 eqs. (77) -> (80).
      uc=0.0
      rsum=0.0
*       Eq. (77).
      DO i=1,n-1
      ks = 4*(i-1)
      k1 = ks + 1
      k2 = ks + 2
      k3 = ks + 3
      k4 = ks + 4
*       Obtain W = L^T(Q)P ; L^T = transpose of L times P.
      w(k1)=(q(k1)*p(k1)-q(k2)*p(k2)-q(k3)*p(k3)+q(k4)*p(k4))
      w(k2)=(q(k2)*p(k1)+q(k1)*p(k2)-q(k4)*p(k3)-q(k3)*p(k4))
      w(k3)=(q(k3)*p(k1)+q(k4)*p(k2)+q(k1)*p(k3)+q(k2)*p(k4))
      w(k4)=(q(k4)*p(k1)-q(k3)*p(k2)+q(k2)*p(k3)-q(k1)*p(k4))
      rijl=q(k1)**2+q(k2)**2+q(k3)**2+q(k4)**2
*       Form RSUM for decision-making.
      rsum=rsum+rijl
      rinv(i)=1./rijl
      a=.5d0*rinv(i)
      uc=uc+mkk(i)*rinv(i)
      DO k=1,4
      w(ks+k)=a*w(ks+k)
      END DO
      END DO
      lri=n-1
*       Eq. (78).
      tkin=0.0
      DO i=1,n-1
      aux(-2)=0.5d0*tk2(i-1)
      aux(-1)=0.5d0*tk1(i)
      aux( 0)=tkk(i)
      aux(+1)=0.5d0*tk1(i+1)
      aux(+2)=0.5d0*tk2(i+1)
      l=4*(i-1)
      dk(i)=0.0
      DO k=1,4
      aa=0.0
      DO j=-2,2
      lj=l+4*j
      if (lj.ge.0.and.lj.le.4*n-8) then
      aa=aa+aux(j)*w(lj+k)
      end if
      END DO
      ak(l+k)=aa
*       Eq. (79).
      dk(i)=dk(i)+aa*w(l+k)
      END DO
*       Eq. (80).
      tkin=tkin+dk(i)
      END DO
*
*       Obtain physical coordinates.
      DO k=1,3
      xi(k)=0.0
      END DO
      DO i=1,n-1
      l=3*(i-1)
      ks=4*(i-1)
      xc(l+1)=q(ks+1)**2-q(ks+2)**2-q(ks+3)**2+q(ks+4)**2
      xc(l+2)=2.d0*(q(ks+1)*q(ks+2)-q(ks+3)*q(ks+4))
      xc(l+3)=2.d0*(q(ks+1)*q(ks+3)+q(ks+2)*q(ks+4))
      DO k=1,3
      xi(l+3+k)=xi(l+k)+xc(l+k)
      END DO
      END DO
*
*       Evaluate external force (FXTNL = W', not the 'usual' force!).
      IF (ipert.GT.0) THEN
      iskip=0
      CALL xtf(fxtnl,fcm,cmx,chtime)
*       Note that perturbation limit must be consistent with chlist.f.
      IF (gpert.LT.1.0d-06) THEN
      ipert=0
      ELSE
      ipert=1
      END IF
      ELSE
      iskip=1
      DO i=1,3*(n-1)
      fxtnl(i)=0.0d0
      END DO
      DO i=1,3
      fcm(i)=0.0d0
      END DO
      END IF
*
*   Form non-chained contributions.
      unc=0.0
      DO i=1,3*(n-1)
      fnc(i)=fxtnl(i)
      END DO
*
      DO i=1,n-2
      li=3*(i-1)
      DO j=i+2,n
      lj=3*(j-1)
      rij2=0.0
      IF(j.GT.i+2)THEN
      DO k=1,3
      xaux(k)=xi(lj+k)-xi(li+k)
      rij2=rij2+xaux(k)**2
      END DO
      ELSE
      DO k=1,3
      xaux(k)=xc(li+k)+xc(li+k+3)
      rij2=rij2+xaux(k)**2
      END DO
      END IF
      rij2in=1./rij2
*       Introduce the inverse distances.
      lri=lri+1
      rinv(lri)=sqrt(rij2in)
*
      fm=mij(i,j)*rinv(lri)
      unc=unc+fm
      fm=fm*rij2in
*       Fij attraction.
      DO k=1,3
      faux(k)=-fm*xaux(k)
      END DO
*       Add the contribution to interactions depending on Rij.
      DO ik=i,j-1
      l=3*(ik-1)
      DO k=1,3
      fnc(l+k)=fnc(l+k)+faux(k)
      END DO
      END DO
      END DO
      END DO
*
*   Evaluate UQ & TP.
      DO i=1,n-1
      l1=3*(i-1)+1
      ks=4*(i-1)
      ks1=ks+1
      CALL qforce(q(ks1),fnc(l1),uq(ks1))
      CALL vector(q(ks1),ak(ks1),tp(ks1))
*       The * operation of eq. (84).
      ak(ks+4)=-ak(ks+4)
      CALL vector(p(ks1),ak(ks1),tq(ks1))
*
      DO k=1,4
      uq(ks+k)=uq(ks+k)-2.0d0*mkk(i)*q(ks+k)*rinv(i)**2
      tq(ks+k)=tq(ks+k)-4.d0*dk(i)*q(ks+k)
      END DO
      END DO
*   NOTE: The division by R above (in TP & TQ) is delayed.
*
*   Proceed to final evaluation of derivatives (90)->(94).
      upot=uc+unc
      g=1./(tkin+upot)
      h=tkin-upot
*       Reset EJUMP after significant slow-down change (Seppo's procedure).
*     IF (KJUMP) THEN
*         EJUMP = H - ENERGY
*         KJUMP = .false.
*     END IF
      gamma=(h-(energy+ejump))*g
*
      gt= (1.-gamma)*g
      gu=-(1.+gamma)*g
*
      DO i=1,n-1
      ks=4*(i-1)
*       Apply the division by R here (to TP & TQ).
      gtoverr=gt*rinv(i)
*       NOTE: TP & TQ never get 'correct' values (thus TP = R*TPtrue).
      DO k=1,4
      dq(ks+k)=gtoverr*tp(ks+k)
      dp(ks+k)=-gtoverr*tq(ks+k)-gu*uq(ks+k)
      END DO
      END DO
      dt=g
      tpr=g
      DO k=1,3
      dx(k)=cmv(k)*g
      dv(k)=fcm(k)*g
      END DO
*   Evaluate E'.
      de=0.0
      IF (iskip.EQ.0) THEN
      DO i=1,n-1
      l1=3*(i-1)+1
      ks=4*(i-1)
      ks1=ks+1
      CALL qforce(q(ks1),fxtnl(l1),faux)
      DO k=1,4
      de=de+dq(ks+k)*faux(k)
      END DO
      END DO
      END IF
*
*       Copy the time derivative for step control and reset indicator.
*     IF (IDER.GT.0) THEN
*         TPR = G
*         IDER = 0
*     END IF
*
*       Check osculating pericentre of closest pair (first call only).
      IF (icall.EQ.0) go to 50
      IF (jc.GT.0) go to 10
*
*       Perform a fast pericentre search (saves unnecessary complications).
      rm = 0.0
      DO 5 i = 1,n-1
          IF (rinv(i).GT.rm) THEN
              rm = rinv(i)
              im = i
          END IF
    5 CONTINUE
*
      k1 = iname(im)
      k2 = iname(im+1)
      k = 4*(im - 1) + 1
      CALL peri(q(k),dq(k),tpr/ksch(im),m(k1),m(k2),qperi)
*
*       Switch off indicator and exit for large pericentre.
      IF (qperi.GT.4.0*max(SIZE(k1),SIZE(k2))) THEN
          icall = 0
          go to 50
      END IF
*
*       Examine all close pair indices K1 & K2 in chain vector INAME.
   10 qpmin = 100.0
      rpmin = 100.0
      im0 = 1
      DO 20 im = 1,n-1
          k1 = iname(im)
          k2 = iname(im+1)
*
*       Set distance of nearest perturber (note extra test for N > 4).
          IF (im.EQ.1) THEN
              rp = 1.0/rinv(2)
          ELSE IF (im.EQ.2.AND.n.GT.3) THEN
              rp = min(1.0/rinv(1),1.0/rinv(3))
          ELSE IF (im.EQ.3.AND.n.GT.4) THEN
              rp = min(1.0/rinv(2),1.0/rinv(4))
          ELSE
              rp = 1.0/rinv(im-1)
          END IF
          rpmin = min(rp,rpmin)
*
*       Determine pericentre for small perturbations by Mikkola's algorithm.
          gi = (1.0/(rinv(im)*rp))**3
          IF (gi.LT.0.005) THEN
              k = 4*(im - 1) + 1
              CALL peri(q(k),dq(k),tpr/ksch(im),m(k1),m(k2),qperi)
          ELSE
              qperi = 1.0/rinv(im)
          END IF
*
*       Compare pericentre with previous mutual distances (note symmetry).
          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 (qperi.LT.qpmin.AND.im.NE.imcirc) THEN
              g0 = gi
              qpmin = qperi
              rp0 = rp
              k10 = k1
              k20 = k2
              im0 = im
              icall = 0
          END IF
   20 CONTINUE
*
*       Check smallest pericentre of current chain and copy latest value.
      rcoll = min(rcoll,qpmin)
      qperi = qpmin
*
*       Specify KS index and closest pair indices.
      im = im0
      ks1 = 4*(im - 1) + 1
      k1 = k10
      k2 = k20
*
*       Save TPR and current configuration (for EREL, TRANSK & CHAIN).
      tkpr = tpr
      DO 25 i = 1,n-1
          ks = 4*(i - 1)
          DO 24 j = 1,4
              qk(ks+j) = q(ks+j)
              pk(ks+j) = p(ks+j)
   24     CONTINUE
   25 CONTINUE
*
*       Check for tidal two-body interaction or stellar collision.
      IF (qpmin.LT.4.0*max(SIZE(k1),SIZE(k2)).OR.iter.GT.0) THEN
*
*       Terminate after 25 iterations (convergence problem).
          iter = iter + 1
          IF (iter.GE.25) THEN
              WRITE (6,26)  iter, im, g0, qpmin, ecc
   26         FORMAT (' WARNING!    NO CONVERGENCE    # IM GI QP ECC ',
     &                                                i5,i4,1p,3e10.2)
              jc = 0
              kcoll = .false.
              icoll = -1
              iter = 0
              imcirc = 0
              go to 50
          END IF
*
*       Exit if collision candidate distance is not the smallest.
          rb = 1.0/rinv(im)
*         IF (RB.GT.1.001*RPMIN) THEN
*             JC = 0
*             DSC = 1.0
*             GO TO 50
*         END IF
*
*       Convert Q' to standard KS with T' = R and form radial velocity R'.
          rpr = 0.0d0
          DO 30 j = ks1,ks1+3
              upr(j) = dq(j)*rb*ksch(im)/tpr
              rpr = rpr + 2.0d0*q(j)*upr(j)
   30     CONTINUE
*
*       Determine small semi-major axis from non-singular expressions.
          CALL erel(im,eb,semi)
*
*       Exclude circularized binaries and consider second smallest peri.
          ecc = 1.0 - qpmin/semi
*       Skip exclusion in GR case (denoted by large VSTAR1).
          IF (ecc.LT.0.002.AND.imcirc.EQ.0.AND.vstar1.LT.100.0) THEN
              iter = 0
              imcirc = im
              go to 10
          END IF
*
*       Temporary exit because SEMI < 0 does not converge (tested OK 3/99).
*         IF (SEMI.LT.0.0) THEN
*             JC = 0
*             GO TO 50
*         END IF
*
*       Obtain pericentre time interval from elements & Kepler's equation.
          mb = m(k1) + m(k2)
          CALL tperi(semi,q(ks1),upr(ks1),mb,tper)
*
*       Activate collision indicator & B-S step selector (first time).
          IF (icoll.EQ.0) THEN
              icoll = -1
              jc = 1
              kcoll = .true.
          ELSE
*
*       Check convergence: radial velocity < 1.0E-09 parabolic velocity.
              IF (abs(rpr).LT.1.0e-09*sqrt(2.0d0*mb*rb)) THEN
*       Reset B-S step selector and copy chain index to collision indicator.
                  jc = 0
                  icoll = im
                  ebs = eb
                  iter = 0
                  imcirc = 0
              END IF
          END IF
*
*       Set regularized step for DIFSY1 using accelerated iteration rate.
*         IF (RB.GT.2.0*QPMIN) THEN
*             DSC = 2.0*ABS(TPER)/TPR
*         ELSE
*             DSC = ABS(TPER)/TPR
*         END IF       
*
*       Evaluate regularized pericentre time (Stiefel & Scheifele, p. 85).
          hi = -0.5*mb/semi
*       Note use of Seppo's sign convention (TPER > 0 after peri).
          dsc = 2.0d0*(hi*tper - 0.5d0*rpr)/mb
*       Scale from KS to chain time derivative.
          dsc = dsc*rb/tpr
*
*       Ensure negative step beyond pericentre (case RPR > 0 works OK).
          IF (jc.GT.0.AND.rpr.LT.0.0d0) THEN
              dsc = abs(dsc)
          END IF
*       Restore step to dummy value at the end (not used).
          IF (jc.EQ.0) dsc = 1.0
*
*       Switch off iteration on large perturbation after five tries.
          IF (g0.GT.0.005.AND.iter.GT.5) THEN
              jc = 0
              kcoll = .false.
              icoll = -1
              iter = 0
              imcirc = 0
*       Avoid apocentre region of secondary binary (algorithmic confusion).
          ELSE IF (rb.GT.semi.AND.imcirc.GT.0) THEN
              jc = 0
              iter = 0
              imcirc = 0
              icoll = 0
              kcoll = .false.
          END IF
      ELSE
          jc = 0
*       Enforce restart from saved variables on failed iteration.
          kcoll = .false.
          icoll = -1
          icall = 0
          iter = 0
      END IF
*
*       Increase function call counter.
   50 nfn = nfn + 1
*
      RETURN
*
      END