zare.f

Go to the documentation of this file.

      SUBROUTINE zare(I1,I2,SP)
*
*
*       Zare stability parameter.
*       -------------------------
*
*       Computes the stability parameter (SP) of a 3-body hierarchical
*       binary system with inner binary (I1,I2) and outer body JCOMP.
*       The system is stable to exchange if SP > 1 is returned.
*       Reference:- K. Zare, Celestial Mechanics 16, 35 (1977).
*       Developed by Murray Alexander 1984.
*
      include 'common6.h'
      REAL*8  m1,m2,m3,m(3),am(3),a(6),xcm(3),vcm(3),x3(3,3),xdot3(3,3)
      INTEGER  ic(3,2),iname(3)
      DATA  ic/2,3,1,3,1,2/
      DATA  err/1.0d-6/
*
*
*       Copy masses and global indices I1, I2 & JCOMP (set in IMPACT).
      m1 = body(i1)
      m2 = body(i2)
      m3 = body(jcomp)
      iname(1) = i1
      iname(2) = i2
      iname(3) = jcomp
*
*       Transform to local centre of mass frame.
      DO 1 k = 1,3
          xcm(k) = 0.0
          vcm(k) = 0.0
          i = iname(k)
          m(k) = body(i)
    1 CONTINUE
      DO 3 l = 1,3
          i = iname(l)
          DO 2 k = 1,3
              xcm(k) = xcm(k) + m(l)*x(k,i)
              vcm(k) = vcm(k) + m(l)*xdot(k,i)
    2     CONTINUE
    3 CONTINUE
      DO 5 l = 1,3
          i = iname(l)
          DO 4 k = 1,3
              x3(k,l) = x(k,i) - xcm(k)/(m1 + m2 + m3)
              xdot3(k,l) = xdot(k,i) - vcm(k)/(m1 + m2 + m3)
    4     CONTINUE
    5 CONTINUE
*
*       Evaluate the total angular momentum, kinetic & potential energy.
      zk3 = 0.0
      DO 10 ia = 1,3
          am(ia) = 0.0d0
          zk3 = zk3 + m(ia)*(xdot3(1,ia)**2 + xdot3(2,ia)**2
     &                                      + xdot3(3,ia)**2)
   10 CONTINUE
      pot3 = 0.0
      DO 30 ia = 1,3
          ja = ic(ia,1)
          ka = ic(ia,2)
          rij2 = 0.0
          DO 20 i = 1,3
              am(ia) = am(ia) + m(i)*(x3(ja,i)*xdot3(ka,i) -
     &                                x3(ka,i)*xdot3(ja,i))
              rij2 = rij2 + (x3(i,ja) - x3(i,ka))**2
   20     CONTINUE
          pot3 = pot3 - m(ja)*m(ka)/sqrt(rij2)
   30 CONTINUE
*
*       Set the bifurcation parameter h*c**2.
      energy = 0.5*zk3 + pot3
      c1 = 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)
*       Note bug fix by Douglas Heggie 1/11/2000.
      a(4) = - (3.0d0*m2 + m1)
      a(5) = - (3.0d0*m2 + 2.0d0*m1)
      a(6) = - (m2 + m1)
*
*       Solve quintic equation for the central collinear configuration.
      icount = 0
      s = 1.0d0
*  50 F1 = A(1)
*     FP1 = F1*5.0D0
*     DO 60 I = 2,5
*         F1 = F1*S + A(I)
*         FP1 = FP1*S + (6-I)*A(I)
*  60 CONTINUE
*
*       Replace by iteration of f(s)/s**2 = 0 for faster convergence (DCH).
   50 f1 = ((a(1)*s + a(2))*s + a(3))*s + a(4) + a(5)/s + a(6)/s**2
      fp1 = (3.0*a(1)*s + 2.0*a(2))*s + a(3) - (2.0*a(6)/s + a(5))/s**2
*
*     F1 = F1*S + A(6)
      dx = -f1/fp1
      IF (abs(dx).GT.err) THEN
          s = s + dx
          icount = icount + 1
          IF (icount.LT.20) go to 50
          WRITE (6,70)  s
   70     FORMAT (' WARNING!    ZARE    NO CONVERGENCE    S =',f8.4)
      END IF
*
*       Compute critical value CCRIT of C1.
      y = 1.0d0 + s
      fx = m1*m3 + m2*(m3/y + m1/s)
      gx = m1*m3 + m2*(m3*y**2 + m1*s**2)
      ccrit = -0.5d0*fx**2*gx/(m1 + m2 + m3)
*
*       Form stability parameter SP.
      sp = c1/ccrit
*
      RETURN
*
      END