histab.f

Go to the documentation of this file.

      SUBROUTINE histab(IPAIR,J,PMIN,RSTAB)
*
*
*       Hierarchical stability condition.
*       ---------------------------------
*
      include 'common6.h'
      REAL*8  xx(3,3),vv(3,3)
*
*
*       Define KS & c.m. indices and save IPAIR for RETURN.
      i1 = 2*ipair - 1
      i2 = i1 + 1
      i = n + ipair
      jpair = ipair
*
*       Evaluate semi-major axis & eccentricity of inner binary.
      semi = -0.5d0*body(i)/h(ipair)
      ecc2 = (1.0d0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(body(i)*semi)
      ecc = sqrt(ecc2)
*
*       Form stability ratio (Eggleton & Kiseleva, Ap.J. 455, 640, 1995).
*     Q = BODY(I)/BODY(J)
*     Q1 = MAX(BODY(I2)/BODY(I1),BODY(I1)/BODY(I2))
*     Q3 = Q**0.33333
*     Q13 = Q1**0.33333
*       Note sign error in third term of Eq. (2).
*     AR = 1.0 + 3.7/Q3 - 2.2/(1.0 + Q3) + 1.4/Q13*(Q3 - 1.0)/(Q3 + 1.0)
*
*       Save hierarchical stability separation in COMMON variable.
*     PCRIT = AR*SEMI*(1.0D0 + ECC)
*
*       Determine the outer eccentricity.
      rij2 = 0.0
      vij2 = 0.0
      rdot = 0.0
      DO 5 k = 1,3
          rij2 = rij2 + (x(k,i) - x(k,j))**2
          vij2 = vij2 + (xdot(k,i) - xdot(k,j))**2
          rdot = rdot + (x(k,i) - x(k,j))*(xdot(k,i) - xdot(k,j))
    5 CONTINUE
      rij = sqrt(rij2)
      a1 = 2.0/rij - vij2/(body(i) + body(j))
*
*       Exit on hyperbolic orbit.
      IF (a1.LT.0.0) THEN
          pmin = 1.0
          rstab = 0.0
          go to 20
      END IF
*
      semi1 = 1.0/a1
      ecc2 = (1.0 - rij/semi1)**2 + rdot**2/(semi1*(body(i) + body(j)))
      ecc1 = sqrt(ecc2)
      pmin = semi1*(1.0 - ecc1)
*
*       Evaluate the basic stability condition without fudge factor.
      q = body(j)/body(i)
      IF (ecc1.LT.1.0) THEN
          xfac = (1.0 + q)*(1.0 + ecc1)/sqrt(1.0 - ecc1)
      ELSE
          xfac = 40.0*(1.0 + q)
      END IF
      pcrit = 2.8*xfac**0.4*semi
*
*       Exit if stability value falls outside practical limits.
      IF ((pcrit.GT.1.5*pmin.OR.pcrit.LT.0.5*pmin).AND.j.LE.n) THEN
          rstab = pcrit
          go to 20
      END IF
*
*       Choose the most active triple in case of two binaries.
      jj = j
      IF (j.GT.n) THEN
          semi2 = -0.5d0*body(j)/h(j-n)
*       Adopt 10% fudge factor with linear dependence on smallest ratio.
          yfac = 1.0 + 0.1*min(semi2/semi,semi/semi2)
          IF (semi2.GT.semi) THEN
              ecc2 = (1.0 - r(j-n)/semi2)**2 +
     &                             tdot2(j-n)**2/(body(j)*semi2)
              ecc = sqrt(ecc2)
              semiz = semi2
              semi2 = semi
              semi = semiz
              ipair = j - n
              i1 = 2*ipair - 1
              i2 = i1 + 1
              jj = i
              i = j
          END IF
      ELSE
          yfac = 1.0
      END IF
*
*       Resolve weakly perturbed binary (prevent X(K,I1) = X(K,I2)).
      IF (gamma(ipair).LT.gmin.OR.x(1,i1).EQ.x(1,i2)) THEN
          CALL resolv(ipair,1)
      END IF
*
*       Determine inclination for triple configuration (NB! radians).
      DO 10 k = 1,3
          xx(k,1) = x(k,i1)
          xx(k,2) = x(k,i2)
          xx(k,3) = x(k,jj)
          vv(k,1) = xdot(k,i1)
          vv(k,2) = xdot(k,i2)
          vv(k,3) = xdot(k,jj)
  10  CONTINUE
      CALL inclin(xx,vv,x(1,i),xdot(1,i),angle)
*
*       Employ the improved stability criterion for doubtful cases.
      rstab = stability(body(i1),body(i2),body(jj),ecc,ecc1,angle)*semi
*       Note: the present stability routine includes inclination!
*     RSTAB = YFAC*RSTAB
      pcrit = rstab
      ipair = jpair
*
   20 RETURN
*
      END