tstab.f

Go to the documentation of this file.

      SUBROUTINE tstab(I,ECC1,SEMI1,PMIN1,YFAC,ITERM)
*
*
*       Hierarchical stability time.
*       ----------------------------
*
      include 'common6.h'
      common/binary/  cm(4,mmax),xrel(3,mmax),vrel(3,mmax),
     &                hm(mmax),um(4,mmax),umdot(4,mmax),tmdis(mmax),
     &                namem(mmax),nameg(mmax),kstarm(mmax),iflagm(mmax)
*
*
*       Distinguish between absolute and approximate stability.
      iterm = 0
      IF (pmin1.GT.yfac*pcrit.OR.kz(19).EQ.0.OR.kz(27).EQ.0) THEN
          tmdis(nmerge+1) = 1.0d+10
      ELSE IF (name(i).LT.0.OR.name(jcomp).LT.0) THEN
*       Skip double hierarchy (RESET2 uses standard stability criterion).
          iterm = 1
      ELSE IF (pmin1.GT.0.8*yfac*pcrit.OR.
     &        (ecc1.GT.0.99.AND.ecc1.LT.1.0)) THEN
*       Estimate time-scale for long-lived hierarchy.
          nk = 1 + 10.0*ecc1/(1.0 - ecc1)
*       Define termination time in terms of outer periods.
          tk = twopi*semi1*sqrt(abs(semi1)/(body(i) + body(jcomp)))
*       Set disruption time in new merger location for routine KSINT.
          tmdis(nmerge+1) = time + abs(nk*tk)
*         WRITE (6,3)  YFAC,PMIN1,PCRIT,0.8*YFAC*PCRIT,NK*TK
*   3     FORMAT (' TSTAB:    YFAC PMIN1 PCR PTEST N*TK ',F6.2,1P,4E9.2)
*       Note that large NK implies termination by other means.
          IF (ecc1.LT.1.0.AND.nk.LE.1) iterm = 1
*       Modify PCRIT so it becomes < SEMI1*(1 - ECC1)*(1 - 2*PERT).
          IF (pcrit.LT.pmin1) pcrit = 0.999*pmin1/yfac
      ELSE IF (ecc1.LT.1.0) THEN
*       Specify termination in all other cases (but allow hyperbolic orbit).
          iterm = 1
      END IF
*
      RETURN
*
      END