stabl3.f

Go to the documentation of this file.

      SUBROUTINE stabl3(ITERM)
*
*
*       Stability test of three-body system.
*       ------------------------------------
*
      IMPLICIT REAL*8  (a-h,m,o-z)
      common/azreg/  time3,tmax,q(8),p(8),r1,r2,r3,energy,m(3),x3(3,3),
     &               xdot3(3,3),cm(10),c11,c12,c19,c20,c24,c25,
     &               nstep3,name3(3),kz15,kz27
      common/azout/  mb,rb,ri,semi,semi1,e,e1
*
*
*       Transform to physical variables.
      CALL trans3(3)
*
*       Identify index of second binary component and the third body.
      im = 1
      IF (r2.LT.r1) im = 2
      i = 3 - im
*
*       Form scalar product of distance & velocity for body #I.
      ridot = x3(1,i)*xdot3(1,i) + x3(2,i)*xdot3(2,i) +
     &                             x3(3,i)*xdot3(3,i)
*
*       Set distance & radial velocity of body #I with respect to binary.
      mb = m(im) + m(3)
      fac = cm(7)/mb
      ri = sqrt(x3(1,i)**2 + x3(2,i)**2 + x3(3,i)**2)
      ridot = fac*ridot/ri
      ri = fac*ri
*
*       Evaluate orbital elements.
      rdot = 0.0d0
      vrel2 = 0.0d0
      vi2 = 0.0d0
      DO 10 k = 1,3
          rdot = rdot + 
     &               (x3(k,3) - x3(k,im))*(xdot3(k,3) - xdot3(k,im))
          vrel2 = vrel2 + (xdot3(k,3) - xdot3(k,im))**2
          vi2 = vi2 + xdot3(k,i)**2
   10 CONTINUE
*
*       Determine semi-major axis & eccentricity of inner binary.
      rb = min(r1,r2)
      semi = 2.0d0/rb - vrel2/mb
      semi = 1.0/semi
      e = sqrt((1.0d0 - rb/semi)**2 + rdot**2/(semi*mb))
*
*       Form outer semi-major axis & eccentricity.
      vi2 = vi2*fac**2
      semi1 = 2.0d0/ri - vi2/cm(7)
      semi1 = 1.0/semi1
      e1 = sqrt((1.0d0 - ri/semi1)**2 + (ri*ridot)**2/(semi1*cm(7)))
*
*       Obtain standard stability ratio (outer pericentre / inner apocentre).
      ratio = semi1*(1.0d0 - e1)/(semi*(1.0d0 + e))
*
*       Form coefficients for stability test (Valtonen, Vistas Ast 32, 1988).
*     AM = (2.65 + E)*(1.0 + M(I)/MB)**0.3333
*     FM = (2.0*M(I) - MB)/(3.0*MB)
*
*       Expand natural logarithm for small arguments.
*     IF (ABS(FM).LT.0.67) THEN
*         BM = FM*(1.0 - (0.5 - ONE3*FM)*FM)
*     ELSE
*         BM = LOG(1.0D0 + FM)
*     END IF
*       Define mass dependent criterion of Harrington (A.J. 80) & Bailyn.
*     PCRIT = AM*(1.0 + 0.7*BM)*SEMI
*
*       Form hierarchical stability ratio (Kiseleva & Eggleton 1995).
*     Q0 = MB/M(I)
*     Q1 = MAX(M(3)/M(IM),M(IM)/M(3))
*     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*SEMI*(1.0D0 + E)
*
*       Adopt the semi-analytical stability criterion (MA 1997).
      q1 = m(i)/mb
      IF (e1.LT.1.0) THEN
          xfac = (1.0 + q1)*(1.0 + e1)/sqrt(1.0 - e1)
      ELSE
          xfac = 40.0*(1.0 + q1)
      END IF
      pcrit = 2.8*xfac**0.4*semi
      pmin = semi1*(1.0d0 - e1)
*
*       Set negative termination index if system is stable and RB > SEMI.
      IF (pcrit.LT.pmin.AND.e1.LT.1.0.AND.1./rb.LT.1./semi) THEN
          iterm = -1
*       Obtain Zare's stability parameter (valid for small inclinations).
          m1 = m(im)
          m2 = m(3)
          m3 = m(i)
          CALL stablz(m1,m2,m3,sp)
          WRITE (6,20)  semi, semi1, e, e1, ri, ratio, sp, pcrit, pmin
   20     FORMAT ('  STABT:    A A1 E E1 RI RATIO SP PCR PM ',
     &                           1p,2e10.2,0p,2f6.2,f9.5,2f6.2,1p,2e9.1)
*       Terminate if escaper is outside 3*SEMI (> 0).
      ELSE IF (e1.GT.1.0.AND.1./ri.LT.1./(3.0*semi)) THEN
          iterm = -1
      ELSE
          iterm = 0
      END IF
*
      RETURN
*
      END