Nbody6
stablz.f
Go to the documentation of this file.
1  SUBROUTINE stablz(M1,M2,M3,SP)
2 *
3 *
4 * Zare stability parameter.
5 * -------------------------
6 *
7 * Computes the stability parameter (SP) of a 3-body hierarchical
8 * binary system with (M1,M2) constituting the inner binary.
9 * The system is stable if SP > 1 is returned.
10 * Reference:- K. Zare, Celestial Mechanics 16, 35 (1977).
11 * Developed by Murray Alexander 1984.
12 *
13  IMPLICIT REAL*8 (a-h,m,o-z)
14  common/azreg/ time3,tmax,q(8),p(8),r1,r2,r3,energy,m(3),x3(3,3),
15  & xdot3(3,3),cm(10),c11,c12,c19,c20,c24,c25,
16  & nstep3,name3(3),kz15,kz27
17  REAL*8 m1,m2,m3,am(3),a(6)
18  INTEGER ic(3,2)
19  DATA ic/2,3,1,3,1,2/
20  DATA err/1.0d-6/
21 *
22 *
23 * Evaluate the total angular momentum.
24  DO 10 ia = 1,3
25  am(ia) = 0.0d0
26  10 CONTINUE
27  DO 30 ia = 1,3
28  ja = ic(ia,1)
29  ka = ic(ia,2)
30  DO 20 i = 1,3
31  am(ia) = am(ia) + m(i)*(x3(ja,i)*xdot3(ka,i) -
32  & x3(ka,i)*xdot3(ja,i))
33  20 CONTINUE
34  30 CONTINUE
35 *
36 * Set the bifurcation parameter h*c**2 (note ENERGY = 2*h).
37  c1 = 0.5d0*energy*(am(1)**2 + am(2)**2 + am(3)**2)
38  a(1) = m3 + m1
39  a(2) = 3.0d0*m3 + 2.0d0*m1
40  a(3) = 3.0d0*m3 + m1
41  a(4) = - a(3)
42  a(5) = - (3.0d0*m2 + 2.0d0*m1)
43  a(6) = - (m2 + m1)
44 *
45 * Solve quintic equation for the central collinear configuration.
46  icount = 0
47  x = 1.0d0
48  50 f1 = a(1)
49  fp1 = f1*5.0d0
50  DO 60 i = 2,5
51  f1 = f1*x + a(i)
52  fp1 = fp1*x + (6-i)*a(i)
53  60 CONTINUE
54 *
55  f1 = f1*x + a(6)
56  dx = -f1/fp1
57  IF (abs(dx).GT.err) THEN
58  x = x + dx
59  icount = icount + 1
60  IF (icount.LT.20) go to 50
61  WRITE (6,70) x
62  70 FORMAT (5x,'WARNING! NO CONVERGENCE IN STABLZ X =',f8.4)
63  END IF
64 *
65 * Compute critical value CCRIT of C1.
66  y = 1.0d0 + x
67  fx = m1*m3 + m2*(m3/y + m1/x)
68  gx = m1*m3 + m2*(m3*y**2 + m1*x**2)
69  ccrit = -0.5d0*fx**2*gx/(m1 + m2 + m3)
70 *
71 * Form stability parameter SP.
72  sp = c1/ccrit
73 *
74  RETURN
75 *
76  END