stabl4.f

Go to the documentation of this file.

      SUBROUTINE stabl4(ITERM)
*
*
*       Stability test of four-body system.
*       -----------------------------------
*
      IMPLICIT REAL*8  (a-h,m,o-z)
      LOGICAL  switch,gtype,gtype0,aux
      common/creg/  m(4),x(12),xd(12),p(12),q(12),time4,energy,epsr2,
     &              xr(9),w(9),r(6),ta(6),mij(6),cm(10),rmax4,tmax,
     &              ds,tstep,eps,nstep4,name4(4),kz15,kz27,nreg,nfn
      common/tpr/   switch,gtype,gtype0
      common/ksave/  k1,k2
      REAL*8  rc(3),vc(3),rc0(3),vc0(3)
*
*
*       Transform to physical variables (retain value of SWITCH).
      aux = switch
      switch = .false.
      CALL endreg
      switch = aux
*
*       Specify indices of two least dominant bodies (denoted K3 & K4).
      k3 = 0
      DO 1 l = 1,4
          IF (l.EQ.k1.OR.l.EQ.k2) go to 1
          IF (k3.EQ.0) THEN
              k3 = l
          ELSE
              k4 = l
          END IF
    1 CONTINUE
*
*       Initialize scalars for orbital elements.
      rb0 = 0.0d0
      rb = 0.0d0
      rb1 = 0.0d0
      rb2 = 0.0d0
      rdot = 0.0d0
      rdot1 = 0.0d0
      rdot2 = 0.0d0
      vrel2 = 0.0d0
      vrel20 = 0.0d0
      vrel21 = 0.0d0
      vrel22 = 0.0d0
*
*       Define binary masses of smallest & widest pair (K1 & K2 and K3 & K4).
      mb0 = m(k1) + m(k2)
      mb = m(k3) + m(k4)
*
*       Form separations & velocities of MB0, MB and their relative orbit.
      DO 10 k = 1,3
          j1 = 3*(k1-1) + k
          j2 = 3*(k2-1) + k
          j3 = 3*(k3-1) + k
          j4 = 3*(k4-1) + k
          rc0(k) = (m(k1)*x(j1) + m(k2)*x(j2))/mb0
          rc(k) = (m(k3)*x(j3) + m(k4)*x(j4))/mb
          vc0(k) = (m(k1)*xd(j1) + m(k2)*xd(j2))/mb0
          vc(k) = (m(k3)*xd(j3) + m(k4)*xd(j4))/mb
          rb0 = rb0 + (x(j1) - x(j2))**2
          rb = rb + (x(j3) - x(j4))**2
          rb1 = rb1 + (rc(k) - rc0(k))**2
          rb2 = rb2 + (rc0(k) - x(j3))**2
          rdot = rdot + (x(j3) - x(j4))*(xd(j3) - xd(j4))
          rdot1 = rdot1 + (rc(k) - rc0(k))*(vc(k) - vc0(k))
          rdot2 = rdot2 + (rc0(k) - x(j3))*(vc0(k) - xd(j3))
          vrel2 = vrel2 + (xd(j3) - xd(j4))**2
          vrel20 = vrel20 + (xd(j1) - xd(j2))**2
          vrel21 = vrel21 + (vc(k) - vc0(k))**2
          vrel22 = vrel22 + (vc0(k) - xd(j3))**2
   10 CONTINUE
*
*       Determine semi-major axis of inner binary.
      rb0 = sqrt(rb0)
      semi0 = 2.0/rb0 - vrel20/mb0
      semi0 = 1.0/semi0
      e0 = 1.0 - rb0/semi0
*
*       Form semi-major axis & eccentricity of outer pair.
      rb = sqrt(rb)
      semi = 2.0d0/rb - vrel2/mb
      semi = 1.0/semi
*     E = SQRT((1.0D0 - RB/SEMI)**2 + RDOT**2/(SEMI*MB))
*
*       Evaluate orbital elements of relative c.m. motion.
      rb1 = sqrt(rb1)
      semi1 = 2.0d0/rb1 - vrel21/cm(7)
      semi1 = 1.0/semi1
      e1 = sqrt((1.0d0 - rb1/semi1)**2 + rdot1**2/(semi1*cm(7)))
*
*       Consider the inner triple.
      mb2 = mb0 + m(k3)
      rb2 = sqrt(rb2)
      semi2 = 2.0d0/rb2 - vrel22/mb2
      semi2 = 1.0/semi2
      e2 = sqrt((1.0d0 - rb2/semi2)**2 + rdot2**2/(semi2*mb2))
*
*       Obtain standard stability ratio (outer pericentre / inner apocentre).
      ratio = semi1*(1.0d0 - e1)/(semi0*(1.0d0 + e0))
*
*       Form coefficients for stability test (Valtonen, Vistas Ast 32, 1988).
*     AM = (2.65 + E0)*(1.0 + MB0/MB)**0.3333
*     FM = (2.0*MB0 - MB)/(3.0*MB)
*
*       Expand natural logarithm for small arguments.
*     IF (ABS(FM).LT.0.67) THEN
*         BM = FM*(1.0 - (0.5 - 0.3333*FM)*FM)
*     ELSE
*         BM = LOG(1.0D0 + FM)
*     END IF
*
*       Adopt mass dependent criterion of Harrington (A.J. 80) & Bailyn.
*     PCRIT = AM*(1.0 + 0.7*BM)*SEMI0
*
*       Form hierarchical stability ratio (Kiseleva & Eggleton 1995).
*     Q0 = MB/MB0
*     Q1 = MAX(M(K2)/M(K1),M(K1)/M(K2))
*     Q3 = Q0**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)
*     PCRIT = AR*SEMI0*(1.0D0 + E0)
*
*       Check stability (AM 1997; inner triple or well separated quadruple).
      iterm = 0
      IF (rb1.GT.5.0*rb2.AND.e2.LT.1.0) THEN
          q1 = m(k3)/mb0
          xfac = (1.0 + q1)*(1.0 + e2)/sqrt(1.0 - e2)
          pcrit = 2.8*xfac**0.4*semi0
          pmin = semi2*(1.0 - e2)
          IF (pcrit.LT.pmin) THEN
              iterm = -1
              ratio = semi2*(1.0d0 - e2)/(semi0*(1.0d0 + e0))
              WRITE (6,15)  semi0, semi2, e0, e2, ratio, rb0, rb2,
     &                      pcrit, pmin
   15         FORMAT ('  STABT:    A0 A2 E0 E2 RATIO R0 R2 PCR PM ',
     &                             1p,2e10.2,0p,2f7.3,f6.2,1p,4e9.1)
          END IF
      ELSE IF (rb1.GT.5.0*max(rb0,rb).AND.e1.LT.1.0.AND.
     &         min(semi0,semi).GT.0.0) THEN
*       Choose smallest binary as third body and ignore fudge factor.
          IF (semi.GT.semi0) THEN
              q1 = mb0/mb
              ain = semi
          ELSE
              q1 = mb/mb0
              ain = semi0
          END IF
          xfac = (1.0 + q1)*(1.0 + e1)/sqrt(1.0 - e1)
          pcrit = 2.8*xfac**0.4*ain
          pmin = semi1*(1.0 - e1)
          IF (pcrit.LT.pmin) THEN
              iterm = -1
              WRITE (6,20)  ain, semi1, e0, e1, ratio, rb0, rb1,
     &                      pcrit, pmin
   20         FORMAT ('  STABQ:    AIN A1 E0 E1 RATIO R0 R1 PCR PM ',
     &                             1p,2e10.2,0p,2f7.3,f6.2,1p,4e9.1)
          END IF
      END IF
*
      RETURN
*
      END