extend.f

Go to the documentation of this file.

      SUBROUTINE extend(ISUB)
*
*
*       Size of (un)perturbed subsystem.
*       --------------------------------
*
      include 'common6.h'
      common/chainc/  xc(3,ncmax),uc(3,ncmax),bodyc(ncmax),ich,
     &                listc(lmax)
      common/clump/   bodys(ncmax,5),t0s(5),ts(5),steps(5),rmaxs(5),
     &                names(ncmax,5),isys(5)
*
*
*       Set global index for chain c.m. body.
      ncm = names(1,isub)
      IF (isys(isub).GE.3) THEN
          icm = ich
          IF (name(ich).EQ.0) go to 20
          go to 12
      END IF
*
*       Find global index for c.m. body (unperturbed triple or quad).
      DO 10 i = ifirst,n
          IF (ncm.EQ.name(i)) THEN
              icm = i
              go to 20
          END IF
   10 CONTINUE
*
*       Include safety procedure in case c.m. body not identified.
   12 WRITE (6,15)  isub, ncm
   15 FORMAT (5x,'WARNING!   SUBSYSTEM TERMINATION   ISUB =',i3,
     &                                                     '  NCM =',i5)
      steps(isub) = 0.0d0
      WRITE (6,60)  isub, nch, ich, name(ich)
   60 FORMAT (' EXTEND:   ISUB NCH ICH NAME(ICH) ',4i5)
      go to 40
*
*       Determine the largest perturbing force (M/R**3).
   20 pmax = 0.0
      nnb = list(1,icm)
      DO 30 l = 2,nnb+1
          j = list(l,icm)
          rij2 = (x(1,j) - x(1,icm))**2 + (x(2,j) - x(2,icm))**2 +
     &                                    (x(3,j) - x(3,icm))**2
          pij = body(j)/(rij2*sqrt(rij2))
          IF (pij.GT.pmax) THEN
              pmax = pij
          END IF
   30 CONTINUE
*
*       Choose maximum of dominant perturber and all neighbours at boundary.
      pmax = max(pmax,float(nnb)*bodym/rs(icm)**3)
*
*       Specify maximum size of (un)perturbed system (limit of 100*GSTAR).
      IF (isys(isub).LE.2) THEN
          gstar = gmin
      ELSE
          gstar = 0.01*gmax
      END IF
      rmaxs(isub) = (100.0*gstar*body(icm)/(2.0*pmax))**0.3333
*
*       Update time limit interval unless termination has been signalled.
      IF (steps(isub).GT.0.0d0) THEN
          steps(isub) = step(icm)
      ELSE
*       Set phase indicator for termination (step reduction in STEPS).
          IF (isys(isub).LE.2) iphase = 4
      END IF
*
*     WRITE (6,35)  ICM,NAME(ICM),RMAXS(ISUB),STEPS(ISUB),STEP(ICM)
*  35 FORMAT (' EXTEND:  ICM NM RMAX STEPS STEP ',2I5,1P,3E10.2)
   40 RETURN
*
      END