Nbody6
 All Files Functions Variables
cstab5.f
Go to the documentation of this file.
1  SUBROUTINE cstab5(ITERM)
2 *
3 *
4 * Five-body chain stability test.
5 * -------------------------------
6 *
7  include 'commonc.h'
8  include 'common2.h'
9  REAL*8 m,mb,mb1,mb2,r2(nmx,nmx),xcm(3),vcm(3),xx(3,3),vv(3,3),
10  & xcm2(3),vcm2(3),xb(3),vb(3),m1,m2,m3
11  INTEGER ij(nmx),isort(nmx)
12 *
13 *
14 * Exclude unlikely stability configuration (ratio 1/4 - 4).
15  ir = 0
16  DO 1 i = 1,n-2
17  ratio = rinv(i+1)/rinv(i)
18  IF (ratio.GT.0.25.AND.ratio.LT.4.0) THEN
19  ir = ir + 1
20  END IF
21  1 CONTINUE
22 *
23 * Exit on at least one unfavourable distance ratio.
24  IF (ir.GT.0) THEN
25 * WRITE (6,2) I, (1.0/RINV(K),K=1,N-1)
26 * 2 FORMAT (' CSTAB5 TEST I R ',I4,1P,5E9.1)
27  go to 40
28  END IF
29 *
30 * Ensure a close binary at beginning or end of chain.
31 * RB = 0.04*RSUM
32 * IF (1.0/RINV(1).GT.RB.AND.1.0/RINV(N-1).GT.RB) GO TO 40
33 *
34 * Check stability of three innermost bodies and one perturber.
35  CALL cstab3(iterm)
36  IF (iterm.LT.0) go to 40
37 *
38  CALL hpsort(n-1,rinv,isort)
39  ib2 = isort(2)
40  ib3 = isort(3)
41  ratio = rinv(ib3)/rinv(ib2)
42  IF (ratio.GT.0.25.AND.ratio.LT.4.0) go to 40
43 *
44 * Sort particle separations (I1 & I2 form closest pair).
45  CALL r2sort(ij,r2)
46 *
47 * Save indices with I1 & I2 and I3 & I4 as inner & outer binary.
48  i1 = ij(1)
49  i2 = ij(2)
50  i3 = ij(3)
51  i4 = ij(4)
52  i5 = 0
53  IF (n.EQ.2) THEN
54  i3 = i2
55  i4 = i1
56  ELSE IF (n.EQ.3) THEN
57  i4 = i1
58 * Note N = 4 is already defined correctly (but redundant here).
59  ELSE IF (n.GT.4) THEN
60 * Determine indices of second closest pair (avoid pair I1-I2).
61  rx1 = 1.0
62  rx0 = r2(i1,i2)
63  DO 5 j1 = 1,n
64  IF (j1.EQ.i1.OR.j1.EQ.i2) go to 5
65  DO 4 j2 = j1+1,n
66  IF (j2.EQ.i1.OR.j2.EQ.i2) go to 4
67  IF (r2(j1,j2).LT.rx1.AND.r2(j1,j2).GT.rx0) THEN
68  rx1 = r2(j1,j2)
69  i3 = j1
70  i4 = j2
71  END IF
72  4 CONTINUE
73  5 CONTINUE
74 * Identify remaining single particle(s) by exclusion.
75  DO 8 i = 1,n
76  IF (i.EQ.i1.OR.i.EQ.i2.OR.i.EQ.i3.OR.i.EQ.i4) go to 8
77  IF (i5.EQ.0) THEN
78  i5 = i
79  ELSE
80  i6 = i
81  END IF
82  8 CONTINUE
83  END IF
84 *
85  mb = m(i1) + m(i2)
86  mb2 = m(i3) + m(i4)
87  m5 = m(i5)
88  mb1 = mb + mb2 + m5
89 *
90 * Form orbital parameters with c.m. of I3 & I4 as third body.
91  vrel2 = 0.0d0
92  vrel3 = 0.0d0
93  rdot = 0.0d0
94  d12 = 0.0
95  d15 = 0.0
96  d25 = 0.0
97  v12 = 0.0
98  v15 = 0.0
99  v25 = 0.0
100  DO 10 k = 1,3
101  j1 = 3*(i1-1) + k
102  j2 = 3*(i2-1) + k
103  j3 = 3*(i3-1) + k
104  j4 = 3*(i4-1) + k
105  j5 = 3*(i5-1) + k
106  vrel2 = vrel2 + (v(j1) - v(j2))**2
107  rdot = rdot + (x(j1) - x(j2))*(v(j1) - v(j2))
108  xcm(k) = (m(i1)*x(j1) + m(i2)*x(j2))/mb
109  vcm(k) = (m(i1)*v(j1) + m(i2)*v(j2))/mb
110  xcm2(k) = (m(i3)*x(j3) + m(i4)*x(j4))/mb2
111  vcm2(k) = (m(i3)*v(j3) + m(i4)*v(j4))/mb2
112  vrel3 = vrel3 + (v(j3) - v(j4))**2
113  d12 = d12 + (xcm(k) - xcm2(k))**2
114  d15 = d15 + (xcm(k) - x(j5))**2
115  d25 = d25 + (xcm2(k) - x(j5))**2
116  v12 = v12 + (vcm(k) - vcm2(k))**2
117  v15 = v15 + (vcm(k) - v(j5))**2
118  v25 = v25 + (vcm2(k) - v(j5))**2
119  10 CONTINUE
120 *
121 * Select inner binary and form relevant quantities.
122  rb1 = 0.0
123  vb1 = 0.0
124  rd1 = 0.0
125  IF (d12.LT.min(d15,d25)) THEN
126  semi = 2.0/sqrt(d12) - v12/(mb + mb2)
127 * Q0 = M5/(MB + MB2)
128  m1 = mb
129  m2 = mb2
130  m3 = m5
131  DO 11 k = 1,3
132  j5 = 3*(i5-1) + k
133  xb(k) = (mb*xcm(k) + mb2*xcm2(k))/(mb + mb2)
134  vb(k) = (mb*vcm(k) + mb2*vcm2(k))/(mb + mb2)
135  rd1 = rd1 + (xb(k) - x(j5))*(vb(k) - v(j5))
136  rb1 = rb1 + (xb(k) - x(j5))**2
137  vb1 = vb1 + (vb(k) - v(j5))**2
138  xx(k,1) = xcm(k)
139  vv(k,1) = vcm(k)
140  xx(k,2) = xcm2(k)
141  vv(k,2) = vcm2(k)
142  xx(k,3) = x(j5)
143  vv(k,3) = v(j5)
144  11 CONTINUE
145  ELSE IF (d15.LT.min(d12,d25)) THEN
146  semi = 2.0/sqrt(d15) - v15/(mb + m5)
147 * Q0 = MB2/(MB + M5)
148  m1 = mb
149  m2 = m5
150  m3 = mb2
151  DO 12 k = 1,3
152  j5 = 3*(i5-1) + k
153  xb(k) = (mb*xcm(k) + m5*x(j5))/(mb + m5)
154  vb(k) = (mb*vcm(k) + m5*v(j5))/(mb + m5)
155  rd1 = rd1 + (xb(k) - xcm2(k))*(vb(k) - vcm2(k))
156  rb1 = rb1 + (xb(k) - xcm2(k))**2
157  vb1 = vb1 + (vb(k) - vcm2(k))**2
158  xx(k,1) = xcm(k)
159  vv(k,1) = vcm(k)
160  xx(k,2) = x(j5)
161  vv(k,2) = v(j5)
162  xx(k,3) = xcm2(k)
163  vv(k,3) = vcm2(k)
164  12 CONTINUE
165  ELSE
166  semi = 2.0/sqrt(d25) - v25/(mb2 + m5)
167 * Q0 = MB/(MB2 + M5)
168  m1 = mb2
169  m2 = m5
170  m3 = mb
171  DO 13 k = 1,3
172  j5 = 3*(i5-1) + k
173  xb(k) = (mb2*xcm2(k) + m5*x(j5))/(mb2 + m5)
174  vb(k) = (mb2*vcm2(k) + m5*v(j5))/(mb2 + m5)
175  rd1 = rd1 + (xb(k) - xcm(k))*(vb(k) - vcm(k))
176  rb1 = rb1 + (xb(k) - xcm(k))**2
177  vb1 = vb1 + (vb(k) - vcm(k))**2
178  xx(k,1) = xcm2(k)
179  vv(k,1) = vcm2(k)
180  xx(k,2) = x(j5)
181  vv(k,2) = v(j5)
182  xx(k,3) = xcm(k)
183  vv(k,3) = vcm(k)
184  13 CONTINUE
185  END IF
186 *
187 * Specify inner & outer semi-major axis.
188  semi = 1.0/semi
189  rb1 = sqrt(rb1)
190  semi1 = 2.0/rb1 - vb1/mb1
191  semi1 = 1.0/semi1
192 *
193 * Obtain orbital elements for inner and outer motion.
194  rb0 = sqrt(r2(i1,i2))
195  semi0 = 2.0/rb0 - vrel2/mb
196  semi0 = 1.0/semi0
197  ecc0 = sqrt((1.0d0 - rb0/semi0)**2 + rdot**2/(semi0*mb))
198  ecc1 = sqrt((1.0d0 - rb1/semi1)**2 + rd1**2/(semi1*mb1))
199  rb2 = sqrt(r2(i3,i4))
200  semi2 = 2.0/rb2 - vrel3/mb2
201  semi2 = 1.0/semi2
202 *
203 * Obtain the inclination.
204  CALL inclin(xx,vv,xb,vb,alpha)
205 *
206 * Employ the general stability criterion.
207  pcrit = stability(m1,m2,m3,ecc0,ecc1,alpha)*semi
208 *
209 * Modify correction factor by widest binary (cf. routine IMPACT).
210  IF (semi2.GT.0.0) pcrit = (1.0 + 0.1*semi2/semi)*pcrit
211 *
212 * Check hierarchical stability condition for bound close pair.
213  pmin = semi1*(1.0d0 - ecc1)
214  iterm = 0
215  IF (pmin.GT.pcrit.AND.semi.GT.0.0.AND.semi1.GT.0.0.AND.
216  & rb0.GT.semi0.AND.semi2.GT.0.0) THEN
217  iterm = -1
218  WRITE (6,20) ecc0, ecc1, semi, semi1, pmin, pcrit, semi2,
219  & 180.0*alpha/3.14
220  20 FORMAT (' CSTAB5 E0 =',f6.3,' E1 =',f6.3,' A =',1p,e8.1,
221  & ' A1 =',e8.1,' PM =',e9.2,' PC =',e9.2,
222  & ' A2 =',e8.1,' IN =',0p,f6.1)
223  END IF
224 *
225  40 RETURN
226 *
227  END