rchain.f

Go to the documentation of this file.

      SUBROUTINE rchain(IQ)
*
*
*       Regularized integration.
*       ------------------------
*
      IMPLICIT REAL*8  (a-h,o-z)
      REAL*8  m,mij,y(25),savex(12),savexd(12)
      LOGICAL  switch,gtype,gtype0
      common/creg/  m(4),x(12),xd(12),p(12),q(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/ind6/  ind(6)
      common/close/  rij(4,4),rcoll,qperi,SIZE(4),ecoll3,ip(4)
      common/ccoll/  qk(12),pk(12),icall,icoll,ndiss4
*     EQUIVALENCE  (Y(1),P(1))
      SAVE
*
*
      DO 3 i = 1,12
          y(i) = p(i)
          y(i+12) = q(i)
    3 CONTINUE
      y(25) = time4
*
      IF (time4.GT.0.0d0) go to 5
      DO 1 j = 1,12
          savex(j) = x(j)
          savexd(j) = xd(j)
    1 CONTINUE
      tsave = time4
      tstep0 = tstep
      istart = 0
      next = 0
      jc = 0
      s = 0.0d0
    2 zmi1 = 100.0
      zmi2 = 100.0
      zmi3 = 100.0
*
      gtype = .false.
      istart = istart + 1
      runav = 0.0001*rmax4
      epsi = eps
      rstar = sqrt(epsr2)
*
    5 DO 50 isteps = 1,10000
          gtype0 = gtype
          dso = ds
          time0 = time4
*
*       Advance the solution one step and order the distances.
          CALL difsy4(25,epsi,ds,s,y)
          DO 10 i = 1,12
              p(i) = y(i)
              q(i) = y(i+12)
   10     CONTINUE
          time4 = y(25)
          nstep4 = nstep4 + 1
          CALL rsort(r,switch,ind)
*
*       Check time transformation type.
          triple = (r(1) + r(2))*(r(2) + r(3))
          IF (triple.LT.epsr2) THEN
              gtype = .true.
          ELSE
              gtype = .false.
          END IF
*
*       Update smallest moment of inertia during forward integration.
          IF (time4.GT.time0.AND.jc.EQ.0) THEN
              imin = ind(1)
              zmi = r(imin)**2
              zmi1 = zmi2
              zmi2 = zmi3
              zmi3 = zmi
              runav = 0.9*runav + 0.1*r(imin)
*       Check termination during last few steps (until R(IM) > SEMI).
              IF (iq.LT.0) go to 100
          END IF
*
*       Switch on search indicator during iteration or just after pericentre.
          IF (icoll.LT.0) icall = 1
          IF (r(imin).LT.rstar.AND.nstep4.GT.next) THEN
              IF (zmi3.GT.zmi2.AND.zmi2.LT.zmi1) THEN
                  icall = 1
              END IF
          END IF
*
*       Delay next search a few steps to avoid the same pericentre.
          IF (icoll.GT.0) THEN
              next = nstep4 + 2
              go to 100
          END IF
*
*       Check the termination criterion.
          im2 = ind(5)
          IF (r(im2).GT.rmax4.OR.time4.GT.tmax) THEN
              IF (r(imin).LT.0.01*runav) go to 50
              iq = -1
              go to 100
          END IF
*
*       Save old time transformation if change of type or new chain.
          IF ((gtype.NEQV.gtype0).OR.switch) THEN
              tpr0 = r(1)*r(2)*r(3)
              IF (gtype0) tpr0 = tpr0/(r(1)*r(2) + (r(1) + r(2))*r(3))
          END IF
*
*       Check the switching condition.
          IF (switch) THEN
              CALL endreg
              CALL newreg
              nreg = nreg + 1
          END IF
*
*       Modify new step after change of type or new chain.
          IF ((gtype.NEQV.gtype0).OR.switch) THEN
              tpr = r(1)*r(2)*r(3)
              IF (gtype) tpr = tpr/(r(1)*r(2) + (r(1) + r(2))*r(3))
*       Scale regularized step by ratio of old & new time derivative.
              ds = ds*tpr0/tpr
          END IF
*
*       Check rare case of zero regularized step.
          IF (isteps.LT.3.AND.ds.EQ.0.0d0) ds = 1.0e-2*dso
          IF (ds.EQ.0.0d0) go to 120
   50 CONTINUE
*
  100 RETURN 
*
*       Make one restart with reduced tolerance.
  120 IF (istart.GT.1) RETURN
      WRITE (6,125)  nstep4, time4, r(ind(6)), rmax4
  125 FORMAT (5x,' RCHAIN RESTART:  # T R6 RMAX4 ',i5,3f12.6)
      time4 = tsave
      DO 130 j = 1,12
          x(j) = savex(j)
          xd(j) = savexd(j)
  130 CONTINUE
      switch = .false.
      CALL newreg
      epsi = 0.1*eps
      ds = 0.1*tstep0/(r(1)*r(2)*r(3))
      WRITE (6,135)  icall,jc,icoll,ds
  135 FORMAT (' RESTART:   ICALL JC ICOLL FACM DS ',3i4,1p,e9.1)
      icall = 0
      icoll = 0
      jc = 0
      go to 2
*
      END