stablc.f

Go to the documentation of this file.

      SUBROUTINE stablc(IM,ITERM,SEMI)
*
*
*       Chain stability test.
*       ---------------------
*
      include 'commonc.h'
      include 'common2.h'
      common/chreg/  timec,tmax,rmaxc,cm(10),namec(6),nstep1,kz27,kz30
      common/clump/   bodys(nmx,5),t0s(5),ts(5),steps(5),rmaxs(5),
     &                names(nmx,5),isys(5)
      common/ccoll2/  qk(nmx4),pk(nmx4),rij(nmx,nmx),SIZE(nmx),vstar1,
     &                ecoll1,rcoll,qperi,istar(nmx),icoll,isync,ndiss1
      common/incond/  x4(3,nmx),xdot4(3,nmx)
      REAL*8  xcm(3),vcm(3)
*
*
*       Check hierarchical stability of triple using simple criterion.
      IF (n.EQ.3) THEN
          CALL chstab(iterm)
*
*       Examine general case of N > 3.
      ELSE IF (n.GE.4) THEN
          iesc = 0
          jesc = 0
*       Distinguish between synchronous binary and large distance ratio.
          IF (isync.NE.0) THEN
*       Specify removal of binary or single body (IM = 1, N-1; 1 < IM < N-1).
              IF (im.EQ.1) THEN
                  iesc = iname(1)
                  jesc = iname(2)
              ELSE IF (im.EQ.n - 1) THEN
                  iesc = iname(n-1)
                  jesc = iname(n)
              ELSE
                  IF (rinv(1).LT.rinv(n-1)) THEN
                      iesc = iname(1)
                  ELSE
                      iesc = iname(n)
                  END IF
              END IF
          ELSE
*       Check for distant binary or single body.
              IF (im.EQ.1.OR.im.EQ.n-1) THEN
                  jm = min(im+1,n-2)
*       Reduce chain if close binary is hierarchical and at beginning or end.
                  IF (1.0/rinv(jm).GT.50.0*semi) THEN
                      iesc = iname(im)
                      jesc = iname(im+1)
                  END IF
              ELSE IF (im.EQ.2.OR.im.EQ.3) THEN
                  jm = 1
*       Identify index of the maximum end separation.
                  IF (1.0/rinv(jm).LT.1.0/rinv(n-1)) jm = n - 1
                  IF (1.0/rinv(jm).GT.50.0*semi) THEN
                      iesc = iname(jm)
                  END IF
              END IF   
          END IF
*
*       Copy chain variables to standard form for escape test.
          IF (iesc.GT.0) THEN
              lk = 0
              DO 10 l = 1,n
                  DO 5 k = 1,3
                      lk = lk + 1
                      x4(k,l) = x(lk)
                      xdot4(k,l) = v(lk)
    5             CONTINUE
   10         CONTINUE
          END IF
*
*       Check binary escape.
          IF (iesc.GT.0.AND.jesc.GT.0) THEN
*       Form coordinates & velocities of escaping binary (local c.m. frame).
              bcm = m(iesc) + m(jesc)
              ri2 = 0.0
              rdot = 0.0
              DO 20 k = 1,3
                  xcm(k) = (m(iesc)*x4(k,iesc) + m(jesc)*x4(k,jesc))/bcm
                  vcm(k) = (m(iesc)*xdot4(k,iesc) +
     &                      m(jesc)*xdot4(k,jesc))/bcm
                  ri2 = ri2 + xcm(k)**2
                  rdot = rdot + xcm(k)*vcm(k)
   20         CONTINUE
*
*       Convert to relative distance & radial velocity w.r. to inner part.
              fac = mass/(mass - bcm)
              ri = sqrt(ri2)
              rdot = fac*rdot/ri
              ri = fac*ri
              resc = 0.5*rsum
*
*       Employ parabolic escape criterion (reduce if RI > RSUM/2 & RDOT > 0).
              IF (ri.GT.resc.AND.rdot.GT.0.0) THEN
                  IF (rdot**2.LT.2.0*mass/ri) THEN
                      iesc = 0
                  END IF
              ELSE
                  iesc = 0
              END IF
*       Check single body escape.
          ELSE IF (iesc.GT.0.AND.jesc.EQ.0) THEN
*       Form relative distance and radial velocity for single particle.
              ri = sqrt(x4(1,iesc)**2 + x4(2,iesc)**2 + x4(3,iesc)**2)
              rdot = x4(1,iesc)*xdot4(1,iesc) + x4(2,iesc)*xdot4(2,iesc)
     &                                        + x4(3,iesc)*xdot4(3,iesc)
              fac = mass/(mass - m(iesc))
              rdot = fac*rdot/ri
              ri = fac*ri
              resc = 0.5*rsum
*       Ensure that escaper is not close to another body.
              rm = min(1.0/rinv(im),ri)
*
*       Check approximate escape criterion outside RESC (allow RX > 2*RESC).
              IF (rm.GT.resc.AND.rdot.GT.0.0) THEN
                  IF (rdot**2.LT.2.0*mass/ri) THEN
                      er = 0.5*rdot**2 - mass/ri
                      rx = -mass/er
                      IF (er.LT.0.0.AND.rx.LT.2.0*resc) THEN
                          iesc = 0
                      END IF
                  END IF
              ELSE
                  iesc = 0
              END IF
          END IF
*
*       Remove identified binary or single escaper.
          IF (iesc.GT.0) THEN
              WRITE (6,40)  iesc, jesc, namec(iesc), ri, rdot**2,
     &                      2.0*mass/ri
   40         FORMAT (' STABLC:    IESC JESC NAM RI RD2 2*M/R ',
     &                             2i4,i6,1p,3e9.1)
              isub = isys(5)
              CALL reduce(iesc,jesc,isub)
              iterm = -2
          END IF
      END IF
*
      RETURN
*
      END