stablz.f

Go to the documentation of this file.

      SUBROUTINE stablz(M1,M2,M3,SP)
*
*
*       Zare stability parameter.
*       -------------------------
*
*       Computes the stability parameter (SP) of a 3-body hierarchical
*       binary system with (M1,M2) constituting the inner binary.
*       The system is stable if SP > 1 is returned.
*       Reference:- K. Zare, Celestial Mechanics 16, 35 (1977).
*       Developed by Murray Alexander 1984.
*
      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
      REAL*8  m1,m2,m3,am(3),a(6)
      INTEGER  ic(3,2)
      DATA  ic/2,3,1,3,1,2/
      DATA  err/1.0d-6/
*
*
*       Evaluate the total angular momentum.
      DO 10 ia = 1,3
          am(ia) = 0.0d0
   10 CONTINUE
      DO 30 ia = 1,3
          ja = ic(ia,1)
          ka = ic(ia,2)
          DO 20 i = 1,3
              am(ia) = am(ia) + m(i)*(x3(ja,i)*xdot3(ka,i) -
     &                                x3(ka,i)*xdot3(ja,i))
   20     CONTINUE
   30 CONTINUE
*
*       Set the bifurcation parameter h*c**2 (note ENERGY = 2*h).
      c1 = 0.5d0*energy*(am(1)**2 + am(2)**2 + am(3)**2)
      a(1) = m3 + m1
      a(2) = 3.0d0*m3 + 2.0d0*m1
      a(3) = 3.0d0*m3 + m1
      a(4) = - a(3)
      a(5) = - (3.0d0*m2 + 2.0d0*m1)
      a(6) = - (m2 + m1)
*
*       Solve quintic equation for the central collinear configuration.
      icount = 0
      x = 1.0d0
   50 f1 = a(1)
      fp1 = f1*5.0d0
      DO 60 i = 2,5
          f1 = f1*x + a(i)
          fp1 = fp1*x + (6-i)*a(i)
   60 CONTINUE
*
      f1 = f1*x + a(6)
      dx = -f1/fp1
      IF (abs(dx).GT.err) THEN
          x = x + dx
          icount = icount + 1
          IF (icount.LT.20) go to 50
          WRITE (6,70)  x
   70     FORMAT (5x,'WARNING!   NO CONVERGENCE IN STABLZ   X =',f8.4)
      END IF
*
*       Compute critical value CCRIT of C1.
      y = 1.0d0 + x
      fx = m1*m3 + m2*(m3/y + m1/x)
      gx = m1*m3 + m2*(m3*y**2 + m1*x**2)
      ccrit = -0.5d0*fx**2*gx/(m1 + m2 + m3)
*
*       Form stability parameter SP.
      sp = c1/ccrit
*
      RETURN
*
      END