cstab3.f

Go to the documentation of this file.

      SUBROUTINE cstab3(ITERM)
*
*
*       Perturbed 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)
      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)
      mb4 = mb1 + m(i4)
*
*       Form output diagnostics.
      vrel2 = 0.0d0
      vrel21 = 0.0d0
      rdot = 0.0d0
      rdot3 = 0.0d0
      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
          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))
          xc3 = (m(i3)*x(j3) + mb*xcm(k))/mb1
          vc3 = (m(i3)*v(j3) + mb*vcm(k))/mb1
          r4 = r4 + (x(j4) - xc3)**2
          v4 = v4 + (v(j4) - vc3)**2
          rdot4 = rdot4 + (x(j4) - xc3)*(v(j4) - vc3)
          xx(k,1) = x(j1)
          xx(k,2) = x(j2)
          xx(k,3) = x(j3)
          vv(k,1) = v(j1)
          vv(k,2) = v(j2)
          vv(k,3) = v(j3)
   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))
      r4 = sqrt(r4)
      semi4 = 2.0/r4 - v4/mb4
      semi4 = 1.0/semi4
      ecc4 = sqrt((1.0 - r4/semi4)**2 + rdot4**2/(semi4*mb4))
*
*       Form hierarchical stability ratio (Eggleton & Kiseleva 1995).
      ql = mb/m(i3)
      q1 = max(m(i2)/m(i1),m(i1)/m(i2))
      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)
*
*       Replace the EK criterion by the MA analytical stability formula.
      q0 = m(i3)/mb
      IF (ecc1.LT.1.0) THEN
          xfac = (1.0 + q0)*(1.0 + ecc1)/sqrt(1.0 - ecc1)
      ELSE
          xfac = 40.0*(1.0 + q0)
      END IF
      fe = 1.0
      pcrit = 2.8*fe*xfac**0.4*semi
*
*       Obtain the inclination.
      CALL inclin(xx,vv,xcm,vcm,alpha)
*
*       Include fudge factor for inclination effect.
      yfac = 1.0 - 0.3*alpha/3.14
      pcrit = yfac*pcrit
*
*       Check hierarchical stability condition for 3 + 1 configuration.
      iterm = 0
      IF (pmin.GT.pcrit.AND.semi1.GT.0.0.AND.semi4.GT.0.0.AND.
     &    rb.GT.semi) THEN
          g4 = 2.0*m(i4)/mb4*(r3/r4)**3
          pmin4 = semi4*(1.0 - ecc4)
          zfac = (1.0 + m(i4)/mb1)*(1.0 + ecc4)/sqrt(1.0 - ecc4)
          zcrit = 2.8*fe*zfac**0.4*semi4
          pz = pmin4/zcrit
*
*       Ignore outermost body for perturbation test of inner triple.
          IF (n.EQ.5) THEN
              g4 = 2.0*m(i4)/mb1*(r3/pmin4)**3
              IF (g4.LT.0.1) THEN
                  WRITE (6,19)  nstep1, ecc4, semi, pcrit, pmin4, g4
   19             FORMAT (' CSTAB3    # E4 A PCR PM4 G4 ',
     &                                i6,f8.4,1p,4e10.2)
                  iterm = -1
              END IF
              go to 40
          END IF
*
*       Continue chain integration if outer orbit unstable or large pert.
          IF (pmin4.LT.zcrit.OR.g4.GT.0.2) go to 40
*
          iterm = -1
          WRITE (6,20)  ecc, ecc1, semi, semi1, pmin, pcrit, ek, pz
   20     FORMAT (' CSTAB3    E =',f6.3,'  E1 =',f6.3,'  A =',1p,e8.1,
     &                     '  A1 =',e8.1,'  PM =',e9.2,'  PC =',e9.2,
     &                     '  EK =',e9.2,'  P/Z =',0p,f5.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