cstab4.f

Go to the documentation of this file.

      SUBROUTINE cstab4(ITERM)
*
*
*       Four-body 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,mb2,r2(nmx,nmx),xcm(3),vcm(3),xx(3,3),vv(3,3),
     &        xcm2(3),vcm2(3)
      INTEGER  ij(nmx)
*
*
*       Sort particle separations (K1 & K2 form closest pair).
      CALL r2sort(ij,r2)
*
*       Save indices with I1 & I2 as wide binary and I3 & I4 as outer body.
      i1 = ij(3)
      i2 = ij(4)
      i3 = ij(1)
      i4 = ij(2)
      mb = m(i1) + m(i2)
      mb2 = m(i3) + m(i4)
      mb1 = mb + mb2
*
*       Use perturbed triple stability instead for small middle distance.
      IF (1.0/rinv(2).LT.0.5*rsum) THEN
          CALL cstab3(iterm)
          go to 40
      END IF
*
*       Form orbital parameters with c.m. of I3 & I4 as third body.
      vrel2 = 0.0d0
      vrel21 = 0.0d0
      vrel34 = 0.0d0
      rdot = 0.0d0
      rdot3 = 0.0d0
      rb0 = 0.0
      ri2 = 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
          xcm2(k) = (m(i3)*x(j3) + m(i4)*x(j4))/mb2
          vcm2(k) = (m(i3)*v(j3) + m(i4)*v(j4))/mb2
          ri2 = ri2 + (xcm2(k) - xcm(k))**2
          vrel21 = vrel21 + (vcm2(k) - vcm(k))**2
          vrel34 = vrel34 + (v(j3) - v(j4))**2
          rdot3 = rdot3 + (xcm2(k) - xcm(k))*(vcm2(k) - vcm(k))
          xx(k,1) = x(j1)
          xx(k,2) = x(j2)
          xx(k,3) = xcm2(k)
          vv(k,1) = v(j1)
          vv(k,2) = v(j2)
          vv(k,3) = vcm2(k)
   10 CONTINUE
*
*       Evaluate orbital elements for inner and outer motion.
      rb = sqrt(r2(i1,i2))
      r3 = sqrt(ri2)
      semi = 2.0d0/rb - vrel2/mb
      semi = 1.0/semi
      ecc = sqrt((1.0d0 - rb/semi)**2 + rdot**2/(semi*mb))
      semi1 = 2.0/r3 - vrel21/mb1
      semi1 = 1.0/semi1
      ecc1 = sqrt((1.0d0 - r3/semi1)**2 + rdot3**2/(semi1*mb1))
      rb2 = sqrt(r2(i3,i4))
      semi2 = 2.0/rb2 - vrel34/mb2
      semi2 = 1.0/semi2
*
*       Skip if innermost triple is not stable (including ECC > 1).
      IF (semi.LT.0.0.OR.ecc.GT.1.0) THEN
          iterm = 1
          go to 40
      END IF
      pcrit0 = stability(m(i1),m(i2),m(i3),ecc,ecc1,0.0d0)*semi
      pmin0 = semi*(1.0 - ecc)
      IF (pmin0.LT.pcrit0) THEN
          iterm = 1
          go to 40
      END IF
*
*       Obtain the inclination.
      CALL inclin(xx,vv,xcm,vcm,alpha)
*
*       Use the general stability formula for the widest binary.
      IF (semi.GT.semi2) THEN
          pcrit = stability(m(i1),m(i2),mb2,ecc,ecc1,alpha)
          ain = semi
      ELSE IF (semi2.GT.0.0) THEN
          pcrit = stability(m(i3),m(i4),mb,0.0d0,ecc1,alpha)
          ain = semi2
      ELSE
          ain = 0.0
          pcrit = 0.0
      END IF
      pcrit = pcrit*(1.0 + 0.1*abs(semi/semi2))*ain
*
*       Check hierarchical stability condition for bound close pair.
      iterm = 0
      pmin = semi1*(1.0d0 - ecc1)
      IF (pmin.GT.pcrit.AND.ain.GT.0.0.AND.semi1.GT.0.0.AND.
     &    rb.GT.semi) THEN
          iterm = -1
          alpha = 180.0*alpha/3.14
          WRITE (6,20)  ecc, ecc1, semi, semi1, pmin, pcrit, alpha
   20     FORMAT (' QUAD HIARCH   E =',f6.3,'  E1 =',f6.3,
     &                     '  A =',1p,e8.1,'  A1 =',e8.1,'  PM =',e9.2,
     &                     '  PC =',e9.2,'  IN =',0p,f6.1)
          ri = sqrt(cm(1)**2 + cm(2)**2 + cm(3)**2)
          emax = 0.0
          WRITE (81,30)  timec, ri, namec(i3), ecc, ecc1, semi, semi1,
     &                   pcrit/pmin, alpha, emax
   30     FORMAT (f9.5,f5.1,i6,2f6.3,1p,2e10.2,0p,f5.2,f6.1,f6.3)
          CALL flush(81)
      END IF
*
   40 RETURN
*
      END