Nbody6
zare.f
Go to the documentation of this file.
1  SUBROUTINE zare(I1,I2,SP)
2 *
3 *
4 * Zare stability parameter.
5 * -------------------------
6 *
7 * Computes the stability parameter (SP) of a 3-body hierarchical
8 * binary system with inner binary (I1,I2) and outer body JCOMP.
9 * The system is stable to exchange if SP > 1 is returned.
10 * Reference:- K. Zare, Celestial Mechanics 16, 35 (1977).
11 * Developed by Murray Alexander 1984.
12 *
13  include 'common6.h'
14  REAL*8 m1,m2,m3,m(3),am(3),a(6),xcm(3),vcm(3),x3(3,3),xdot3(3,3)
15  INTEGER ic(3,2),iname(3)
16  DATA ic/2,3,1,3,1,2/
17  DATA err/1.0d-6/
18 *
19 *
20 * Copy masses and global indices I1, I2 & JCOMP (set in IMPACT).
21  m1 = body(i1)
22  m2 = body(i2)
23  m3 = body(jcomp)
24  iname(1) = i1
25  iname(2) = i2
26  iname(3) = jcomp
27 *
28 * Transform to local centre of mass frame.
29  DO 1 k = 1,3
30  xcm(k) = 0.0
31  vcm(k) = 0.0
32  i = iname(k)
33  m(k) = body(i)
34  1 CONTINUE
35  DO 3 l = 1,3
36  i = iname(l)
37  DO 2 k = 1,3
38  xcm(k) = xcm(k) + m(l)*x(k,i)
39  vcm(k) = vcm(k) + m(l)*xdot(k,i)
40  2 CONTINUE
41  3 CONTINUE
42  DO 5 l = 1,3
43  i = iname(l)
44  DO 4 k = 1,3
45  x3(k,l) = x(k,i) - xcm(k)/(m1 + m2 + m3)
46  xdot3(k,l) = xdot(k,i) - vcm(k)/(m1 + m2 + m3)
47  4 CONTINUE
48  5 CONTINUE
49 *
50 * Evaluate the total angular momentum, kinetic & potential energy.
51  zk3 = 0.0
52  DO 10 ia = 1,3
53  am(ia) = 0.0d0
54  zk3 = zk3 + m(ia)*(xdot3(1,ia)**2 + xdot3(2,ia)**2
55  & + xdot3(3,ia)**2)
56  10 CONTINUE
57  pot3 = 0.0
58  DO 30 ia = 1,3
59  ja = ic(ia,1)
60  ka = ic(ia,2)
61  rij2 = 0.0
62  DO 20 i = 1,3
63  am(ia) = am(ia) + m(i)*(x3(ja,i)*xdot3(ka,i) -
64  & x3(ka,i)*xdot3(ja,i))
65  rij2 = rij2 + (x3(i,ja) - x3(i,ka))**2
66  20 CONTINUE
67  pot3 = pot3 - m(ja)*m(ka)/sqrt(rij2)
68  30 CONTINUE
69 *
70 * Set the bifurcation parameter h*c**2.
71  energy = 0.5*zk3 + pot3
72  c1 = energy*(am(1)**2 + am(2)**2 + am(3)**2)
73  a(1) = m3 + m1
74  a(2) = 3.0d0*m3 + 2.0d0*m1
75  a(3) = 3.0d0*m3 + m1
76 * A(4) = - A(3)
77 * Note bug fix by Douglas Heggie 1/11/2000.
78  a(4) = - (3.0d0*m2 + m1)
79  a(5) = - (3.0d0*m2 + 2.0d0*m1)
80  a(6) = - (m2 + m1)
81 *
82 * Solve quintic equation for the central collinear configuration.
83  icount = 0
84  s = 1.0d0
85 * 50 F1 = A(1)
86 * FP1 = F1*5.0D0
87 * DO 60 I = 2,5
88 * F1 = F1*S + A(I)
89 * FP1 = FP1*S + (6-I)*A(I)
90 * 60 CONTINUE
91 *
92 * Replace by iteration of f(s)/s**2 = 0 for faster convergence (DCH).
93  50 f1 = ((a(1)*s + a(2))*s + a(3))*s + a(4) + a(5)/s + a(6)/s**2
94  fp1 = (3.0*a(1)*s + 2.0*a(2))*s + a(3) - (2.0*a(6)/s + a(5))/s**2
95 *
96 * F1 = F1*S + A(6)
97  dx = -f1/fp1
98  IF (abs(dx).GT.err) THEN
99  s = s + dx
100  icount = icount + 1
101  IF (icount.LT.20) go to 50
102  WRITE (6,70) s
103  70 FORMAT (' WARNING! ZARE NO CONVERGENCE S =',f8.4)
104  END IF
105 *
106 * Compute critical value CCRIT of C1.
107  y = 1.0d0 + s
108  fx = m1*m3 + m2*(m3/y + m1/s)
109  gx = m1*m3 + m2*(m3*y**2 + m1*s**2)
110  ccrit = -0.5d0*fx**2*gx/(m1 + m2 + m3)
111 *
112 * Form stability parameter SP.
113  sp = c1/ccrit
114 *
115  RETURN
116 *
117  END