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