invert.f

Go to the documentation of this file.

      SUBROUTINE invert(DT,DTAU)
*
*
*       Inversion of physical time interval.
*       ------------------------------------
*
      include 'commonc.h'
      include 'common2.h'
      LOGICAL  kslow,kcoll
      REAL*8  vi(nmx3),vc(nmx3),ksch,c(5)
      common/slow1/   tk2(0:nmx),ejump,ksch(nmx),kslow,kcoll
      common/ccoll2/  qk(nmx4),pk(nmx4),rij(nmx,nmx),SIZE(nmx),vstar1,
     &                ecoll1,rcoll,qperi,istar(nmx),icoll,isync,ndiss1
      common/chreg/  timec,tmax,rmaxc,cm(10),namec(6),nstep1,kz27,kz30
      common/calls/  tpr,tkpr,step,ider,icall,nfn,nreg,iter,imcirc
      common/clump/   bodys(nmx,5),t0s(5),ts(5),steps(5),rmaxs(5),
     &                names(nmx,5),isys(5)
      common/ksave/  k1,k2
      SAVE
      DATA  iwarn/0/
*
*
*   Transform to current chain coordinates and set RINV.
      DO i=1,n-1
          l1=3*(i-1)+1
          ks1=4*(i-1)+1
          CALL ksphys(q(ks1),p(ks1),xc(l1),wc(l1))
          rk = q(ks1)**2 + q(ks1+1)**2 + q(ks1+2)**2 + q(ks1+3)**2
          rinv(i) = 1.0/rk
      END DO
*
*       Obtain physical velocities (first absolute, then relative).
      l = 3*(n-2)
      DO k = 1,3
          vi(k) = -wc(k)/mc(1)
          vi(l+k+3) = wc(l+k)/mc(n)
      END DO
      DO i = 2,n-1
          l = 3*(i-1)
          DO k = 1,3
              vi(l+k) = (wc(l+k-3) - wc(l+k))/mc(i)
          END DO
      END DO  
      do j = 1,3*(n-1)
          vc(j) = vi(j+3) - vi(j)
      end do
*
*       Determine the smallest two-body distance and chain index.
      rm = 0.0
      do i = 1,n-1
          if (rinv(i).gt.rm) then
              rm = rinv(i)
              i1 = i
          end if
      end do
      rm = 1.0/rm
*
*       Form semi-major axes for the two closest distances.
      dm = 0.0
      do i = 1,n-1
          l = 3*(i - 1)
          if (i.eq.i1) then
              w2 = vc(l+1)**2 + vc(l+2)**2 + vc(l+3)**2
              mb = mc(i) + mc(i+1)
              amax = 2.0*rinv(i) - w2/mb
              a1 = 1.0/amax
              v21 = w2
*       Obtain the largest relative perturbation from nearest neighbour.
              if (i.eq.1) then
                  g1 = 2.0*mc(3)/mb*(rinv(2)/rinv(1))**3
              else
                  g1 = 2.0*mc(i-1)/mb*(rinv(i-1)/rinv(i))**3
                  if (i.lt.n-1) then
                      g1b = 2.0*mc(i+2)/mb*(rinv(i+1)/rinv(i))**3
                      g1 = max(g1b,g1)
                  end if
              end if
*       Treat the second smallest distance similarly.
          else if (rinv(i).gt.dm) then
              w2 = vc(l+1)**2 + vc(l+2)**2 + vc(l+3)**2
              mb = mc(i) + mc(i+1)
              amax = 2.0*rinv(i) - w2/mb
              a2 = 1.0/amax
              i2 = i
              dm = rinv(i)
              if (i.eq.1) then
                  g2 = 2.0*mc(3)/mb*(rinv(2)/rinv(1))**3
              else
                  g2 = 2.0*mc(i-1)/mb*(rinv(i-1)/rinv(i))**3
                  if (i.lt.n-1) then
                      g2b = 2.0*mc(i+2)/mb*(rinv(i+1)/rinv(i))**3
                      g2 = max(g2b,g2)
                  end if
              end if
          end if
      end do
      dm = 1.0/dm
*
*       Obtain reliable semi-major axis for small pericentre or large EB.
      eb = -0.5*mc(i1)*mc(i1+1)/a1
      eb1 = eb/energy
      IF (rm.lt.0.2*a1.or.eb1.gt.0.9) then
*       Copy current configuration for EREL & TRANSK.
          DO 2 i = 1,n-1
              ks = 4*(i - 1)
              DO 1 j = 1,4
                  qk(ks+j) = q(ks+j)
                  pk(ks+j) = p(ks+j)
    1         CONTINUE
    2     CONTINUE
*
*       Evaluate semi-major axis of K1 & K2 from non-singular variables.
          k1 = iname(i1)
          k2 = iname(i1+1)
          CALL erel(i1,eb,semi)
*
          err = (a1 - semi)/semi
          IF (abs(err).GT.1.0e-05.AND.iwarn.LT.10) THEN
              iwarn = iwarn + 1
              zk = 1.0
              DO i = 1,n-1
                  zk = max(ksch(i),zk)
              END DO
              WRITE (6,4)  nstep1, zk, g1, rm, semi, err
    4         FORMAT (' WARNING!    INVERT    # K g1 RM SEMI DA/A ',
     &                                        i6,f7.2,1p,4e10.2)
          END IF
*       Replace direct value by improved determination.
          a1 = semi
      END IF
*
*       Restrict the interval to a fractional period or crossing time.
      IF (a1.GT.0.0) THEN
          tks = 6.28*ksch(i1)*a1*sqrt(a1/(mc(i1) + mc(i1+1)))
*       Adopt upper limit of TKS/2 but allow genuine small DT.
          dt = min(0.5*tks,dt)
      ELSE
*       Use crossing time for hyperbolic velocity (TPR compensates small R).
          dt = min(rm/sqrt(v21),dt)
      END IF
*
*       Check second closest pair if bound (small period; a1 < 0 possible).
      if (a2.gt.0.0) then
          tks = 6.28*ksch(i2)*a2*sqrt(a2/(mc(i2) + mc(i2+1)))
          dt = min(0.5*tks,dt)
      end if
*
*       Initialize time variables.
      dt0 = dt
      sum1 = 1.0/tpr
      sum2 = 0.0
*
*       Skip iteration for significant perturbations.
      IF (g1.GT.0.01) THEN
          dtau = sum1*dt
          go to 30
      END IF
*
*       Begin the loop (second time depends on conditions).
    5 dt = dt0
      mb = mc(i1) + mc(i1+1)
*       Form useful quantities (including EB and eccentricity if needed).
      beta = mb/a1
      li = 3*(i1 - 1)
      eta = xc(li+1)*vc(li+1) + xc(li+2)*vc(li+2) + xc(li+3)*vc(li+3)
      zeta = mb - beta*rm
*     EB = -0.5*MC(I1)*MC(I1+1)/a1
*     ECC = SQRT((1.0 - RM/a1)**2 + ETA**2/(MB*a1))
*
*       Sum the constant part of the integral.
      sum1 = sum1 - 2.0*mc(i1)*mc(i1+1)/(ksch(i1)*rm)
*       Reduce physical time by slow-down factor for the two-body integral.
      dt = dt/ksch(i1)
*
*       Obtain improved initial guess by factor of 2 bisecting (Chris Tout).
      it = 0
*       Choose a conservative starting value in case of small RM.
      y = dt*min(1.0/rm,0.5/abs(a1))
    8 z = beta*y**2
      CALL cfuncs(z,c)
      y0 = dt - ((zeta*c(3)*y + eta*c(2))*y + rm)*y
      it = it + 1
      IF (y0.GT.0.0) THEN
          y = 2.0*y
      ELSE
          y = 0.5*y
      END IF
      IF (it.LT.4) go to 8
*
*       Perform Newton-Raphson iteration using first and second derivatives.
   10 z = beta*y**2
*
*       Obtain new c-functions for argument Z = -2h*Y^2.
      CALL cfuncs(z,c)
*
*       Define the function Y0 and its derivatives (cf. Book eqns. 12.9).
      y0 = dt - ((zeta*c(3)*y + eta*c(2))*y + rm)*y
      ypr = -((zeta*c(2)*y + eta*c(1))*y + rm)
      y2pr = -eta*(1.0 - z*c(2)) - zeta*y*c(1)
*       Adopt safety measure to avoid negative argument (Seppo Mikkola).
      d = y0*y2pr/ypr**2
      y = y - y0/ypr/sqrt(0.01 + abs(0.99d0 - d))
      dt1 = ((zeta*c(3)*y + eta*c(2))*y + rm)*y
      it = it + 1
      IF (abs(dt - dt1).GT.1.0d-06*abs(dt).AND.it.LT.25) go to 10
*
*       Accumulate the integral of dt/R.
      sum2 = sum2 + 2.0*mc(i1)*mc(i1+1)*y
*
      IF (it.GT.15.AND.iwarn.LT.50) THEN
          iwarn = iwarn + 1
          ecc = sqrt((1.0 - rm/a1)**2 + eta**2/(mb*a1))
          WRITE (6,15)  it, ksch(i1), ecc, rm, a1, dt-dt1, dtau
   15     FORMAT (' WARNING!    INVERT   IT K E R A ERR DTAU ',
     &                                   i4,f7.2,f7.3,1p,5e9.1)
      END IF
*
*       Consider the second close interaction for small perturbation.
      IF (i2.ne.i1.and.g2.lt.0.01) THEN
          i1 = i2
          a1 = a2
          rm = dm
          go to 5
      END IF
*
*       Combine the two terms and ensure positive value.
      dtau = sum1*dt0 + sum2
      IF (dtau.LT.0.0) THEN
          dtau = abs(sum1)*dt0
          IF (iwarn.LT.100) THEN
              iwarn = iwarn + 1
              WRITE (6,16)  ksch(i1), rm/a1, sum1*dt0, sum2, dtau
   16         FORMAT (' INVERT NEGATIVE!    KSCH R/A S1*DT0 S2 DTAU ',
     &                                      f7.2,1p,4e10.2)
          END IF
      END IF
*
*     WRITE (6,20)  IT, a1, DT, DT-DT1, DTAU, STEP
*  20 FORMAT (' INVERT:    IT A DT ERR DTAU STEP ',I4,1P,5E10.2)
*
   30 RETURN
*
      END