cstab2.f

Go to the documentation of this file.

      SUBROUTINE cstab2(ITERM)
*
*
*       Degenerate triple chain stability test.
*       ---------------------------------------
*
      include 'commonc.h'
      include 'common2.h'
      common/chreg/  timec,tmax,rmaxc,cm(10),namec(6),nstep1,kz27,kz30
      REAL*8  m,mb,mb1,r2(nmx,nmx),xcm(3),vcm(3),xx(3,3),vv(3,3),
     &        xcm3(3),vcm3(3)
      INTEGER  ij(nmx)
*
*
*       Sort particle separations (I1 & I2 form closest pair).
      CALL r2sort(ij,r2)
      i1 = ij(1)
      i2 = ij(2)
      i3 = ij(3)
      i4 = ij(4)
      mb = m(i1) + m(i2)
      mb1 = mb + m(i3)
      mb2 = mb1 + m(i4)
*
*       Form output diagnostics with smallest binary as one body.
      vrel2 = 0.0d0
      vrel21 = 0.0d0
      rdot = 0.0d0
      rdot3 = 0.0d0
      rb0 = 0.0
      ri2 = 0.0d0
      r4 = 0.0d0
      v4 = 0.0
      rdot4 = 0.0d0
      DO 10 k = 1,3
          j1 = 3*(i1-1) + k
          j2 = 3*(i2-1) + k
          j3 = 3*(i3-1) + k
          j4 = 3*(i4-1) + k
          rb0 = rb0 + (x(j1) - x(j2))**2
          vrel2 = vrel2 + (v(j1) - v(j2))**2
          rdot = rdot + (x(j1) - x(j2))*(v(j1) - v(j2))
          xcm(k) = (m(i1)*x(j1) + m(i2)*x(j2))/mb
          vcm(k) = (m(i1)*v(j1) + m(i2)*v(j2))/mb
          ri2 = ri2 + (x(j3) - xcm(k))**2
          vrel21 = vrel21 + (v(j3) - vcm(k))**2
          rdot3 = rdot3 + (x(j3) - xcm(k))*(v(j3) - vcm(k))
          xcm3(k) = (m(i3)*x(j3) + mb*xcm(k))/mb1
          vcm3(k) = (m(i3)*v(j3) + mb*vcm(k))/mb1
          r4 = r4 + (x(j4) - xcm3(k))**2
          v4 = v4 + (v(j4) - vcm3(k))**2
          rdot4 = rdot4 + (x(j4) - xcm3(k))*(v(j4) - vcm3(k))
          xx(k,1) = xcm(k)
          xx(k,2) = xcm3(k)
          xx(k,3) = x(j4) - xcm3(k)
          vv(k,1) = vcm(k)
          vv(k,2) = vcm3(k)
          vv(k,3) = v(j4) - vcm3(k)
   10 CONTINUE
*
*       Evaluate orbital elements for inner and outer motion.
      rb = sqrt(ri2)
      r4 = sqrt(r4)
      semi = 2.0d0/rb - vrel21/mb1
      semi = 1.0/semi
      ecc = sqrt((1.0d0 - rb/semi)**2 + rdot3**2/(semi*mb1))
      semi1 = 2.0/r4 - v4/mb2
      semi1 = 1.0/semi1
      ecc1 = sqrt((1.0d0 - r4/semi1)**2 + rdot4**2/(semi1*mb2))
*
*       Skip if innermost triple is not stable (including ECC > 1).
      rb0 = sqrt(rb0)
      semi0 = 2.0/rb0 - vrel2/mb
      semi0 = 1.0/semi0
      IF (semi0.LT.0.0.OR.ecc.GT.1.0) go to 40
      ecc0 = sqrt((1.0 - rb0/semi0)**2 + rdot**2/(semi0*mb))
      pcrit0 = stability(m(i1),m(i2),m(i3),ecc0,ecc,0.0d0)*semi0
*       Allow 10% extra with MA 1999 (proper test done for outer triple).
      pmin0 = 1.1*semi*(1.0 - ecc)
      IF (pmin0.LT.pcrit0) go to 40
*
*       Form hierarchical stability ratio (Eggleton & Kiseleva 1995).
*     QL = MB1/M(I4)
*     Q1 = MAX(MB/MB1,MB1/MB)
*     Q3 = QL**0.33333
*     Q13 = Q1**0.33333
*     AR = 1.0 + 3.7/Q3 - 2.2/(1.0 + Q3) + 1.4/Q13*(Q3 - 1.0)/(Q3 + 1.0)
*
*     EK = AR*SEMI*(1.0D0 + ECC)
      pmin = semi1*(1.0d0 - ecc1)
*
*       Obtain the inclination (in radians).
      CALL inclin(xx,vv,xcm3,vcm3,alpha)
*
*       Replace the EK criterion by the MA 1999 stability formula.
      pc99 = stability(mb,m(i3),m(i4),ecc,ecc1,alpha)*semi
*
*       Evaluate the general stability function.
      IF (ecc1.LT.1.0) THEN
          nst = nstab(semi,semi1,ecc,ecc1,alpha,mb,m(i3),m(i4))
          IF (nst.EQ.0) THEN
              pcrit = 0.99*pmin
          ELSE
              pcrit = 1.01*pmin
          END IF
      ELSE
          pcrit = 1.01*pmin
      END IF
*
*       Check hierarchical stability condition (SEMI1 > 0 => ECC1 < 1).
      iterm = 0
      IF (pmin.GT.pcrit.AND.semi.GT.0.0.AND.semi1.GT.0.0) THEN
          iterm = -1
          WRITE (6,20)  ecc, ecc1, semi, semi1, pmin, pcrit, pc99
   20     FORMAT (' CSTAB2    E =',f6.3,'  E1 =',f6.3,'  A =',1p,e8.1,
     &                     '  A1 =',e8.1,'  PM =',e9.2,'  PC =',e9.2,
     &                     '  PC99 =',e9.2)
          WRITE (6,25)  namec(i1), namec(i2), ecc0, semi0, pmin0, pcrit0
   25     FORMAT (' INNER TRIPLE    NAM E0 A0 PM0 PC0 ',
     &                              2i6,f7.3,1p,3e10.2)
          ri = sqrt(cm(1)**2 + cm(2)**2 + cm(3)**2)
          emax = 0.0
          WRITE (81,30)  timec, ri, namec(i3), ql, q1, ecc, ecc1,
     &                   semi, semi1, pcrit/pmin, 180.*alpha/3.14, emax
   30     FORMAT (f9.5,f5.1,i6,2f6.2,2f6.3,1p,2e10.2,0p,f5.2,f6.1,f6.3)
          CALL flush(81)
      END IF
*
   40 RETURN
*
      END