Nbody6
 All Files Functions Variables
chstab.f
Go to the documentation of this file.
1  SUBROUTINE chstab(ITERM)
2 *
3 *
4 * Chain stability test.
5 * ---------------------
6 *
7  include 'commonc.h'
8  include 'common2.h'
9  common/chreg/ timec,tmax,rmaxc,cm(10),namec(6),nstep1,kz27,kz30
10  common/cpert/ rgrav,gpert,ipert,npert
11  REAL*8 m,mb,mb1,r2(nmx,nmx),xcm(3),vcm(3),xx(3,3),vv(3,3),
12  & a1(3),a2(3),xrel(3),vrel(3),ei(3),hi(3),ho(3)
13  INTEGER ij(nmx)
14 *
15 *
16 * Sort particle separations (I1 & I2 form closest pair).
17  CALL r2sort(ij,r2)
18  i1 = ij(1)
19  i2 = ij(2)
20  i3 = ij(3)
21  mb = m(i1) + m(i2)
22  mb1 = mb + m(i3)
23 *
24 * Form output diagnostics.
25  vrel2 = 0.0d0
26  vrel21 = 0.0d0
27  rdot = 0.0d0
28  rdot3 = 0.0d0
29  r3 = 0.0d0
30  DO 5 k = 1,3
31  j1 = 3*(i1-1) + k
32  j2 = 3*(i2-1) + k
33  j3 = 3*(i3-1) + k
34  vrel2 = vrel2 + (v(j1) - v(j2))**2
35  rdot = rdot + (x(j1) - x(j2))*(v(j1) - v(j2))
36  xcm(k) = (m(i1)*x(j1) + m(i2)*x(j2))/mb
37  vcm(k) = (m(i1)*v(j1) + m(i2)*v(j2))/mb
38  r3 = r3 + (x(j3) - xcm(k))**2
39  vrel21 = vrel21 + (v(j3) - vcm(k))**2
40  rdot3 = rdot3 + (x(j3) - xcm(k))*(v(j3) - vcm(k))
41  xx(k,1) = x(j1)
42  xx(k,2) = x(j2)
43  xx(k,3) = x(j3)
44  vv(k,1) = v(j1)
45  vv(k,2) = v(j2)
46  vv(k,3) = v(j3)
47  5 CONTINUE
48 *
49 * Evaluate orbital elements for inner and outer motion.
50  rb = sqrt(r2(i1,i2))
51  r3 = sqrt(r3)
52  semi = 2.0d0/rb - vrel2/mb
53  semi = 1.0/semi
54  ecc = sqrt((1.0d0 - rb/semi)**2 + rdot**2/(semi*mb))
55  semi1 = 2.0/r3 - vrel21/mb1
56  semi1 = 1.0/semi1
57  ecc1 = sqrt((1.0d0 - r3/semi1)**2 + rdot3**2/(semi1*mb1))
58  pmin = semi1*(1.0d0 - ecc1)
59 *
60 * Form hierarchical stability ratio (Eggleton & Kiseleva 1995).
61 * QL = MB/M(I3)
62 * Q1 = MAX(M(I2)/M(I1),M(I1)/M(I2))
63 * Q3 = QL**0.33333
64 * Q13 = Q1**0.33333
65 * AR = 1.0 + 3.7/Q3 - 2.2/(1.0 + Q3) + 1.4/Q13*(Q3 - 1.0)/(Q3 + 1.0)
66 * EK = AR*SEMI*(1.0D0 + ECC)
67 *
68 * Obtain the inclination.
69  CALL inclin(xx,vv,xcm,vcm,alpha)
70 *
71 * Replace the EK criterion by the MA 1999 stability expression.
72 * PCRIT = stability(M(I1),M(I2),M(I3),ECC,ECC1,ALPHA)*SEMI
73 * Add 1% for perturbation to avoid repeated switching.
74 * PCRIT = 1.01*PCRIT
75 *
76 * Evaluate the general stability function (Mardling, Cambody 2008).
77  IF (ecc1.LT.1.0.AND.ecc.LT.1.0) THEN
78  nst = nstab(semi,semi1,ecc,ecc1,alpha,m(i1),m(i2),m(i3))
79  IF (nst.EQ.0) THEN
80  pcrit = 0.99*pmin
81  ELSE
82  pcrit = 1.01*pmin
83  END IF
84  ELSE
85 * Set nominal failed value for SEMI1 > 0 test.
86  pcrit = 1.01*pmin
87  END IF
88 *
89 * Prepare evaluation of maximum eccentricity (see routine HIGROW).
90  DO 10 k = 1,3
91  xrel(k) = xx(k,1) - xx(k,2)
92  vrel(k) = vv(k,1) - vv(k,2)
93  10 CONTINUE
94  a12 = 0.0
95  a22 = 0.0
96  a1a2 = 0.0
97  ri2 = 0.0
98  vi2 = 0.0
99  rvi = 0.0
100  DO 12 k = 1,3
101  k1 = k + 1
102  IF (k1.GT.3) k1 = 1
103  k2 = k1 + 1
104  IF (k2.GT.3) k2 = 1
105  a1(k) = xrel(k1)*vrel(k2) - xrel(k2)*vrel(k1)
106  a2(k) = (xx(k1,3) - xcm(k1))*(vv(k2,3) - vcm(k2))
107  & - (xx(k2,3) - xcm(k2))*(vv(k1,3) - vcm(k1))
108  a12 = a12 + a1(k)**2
109  a22 = a22 + a2(k)**2
110  a1a2 = a1a2 + a1(k)*a2(k)
111  ri2 = ri2 + xrel(k)**2
112  vi2 = vi2 + vrel(k)**2
113  rvi = rvi + xrel(k)*vrel(k)
114  12 CONTINUE
115 *
116 * Construct the Runge-Lenz vector (Heggie & Rasio 1995, Eq.(5)).
117  ei2 = 0.0
118  DO 15 k = 1,3
119  ei(k) = (vi2*xrel(k) - rvi*vrel(k))/mb - xrel(k)/sqrt(ri2)
120  ei2 = ei2 + ei(k)**2
121  15 CONTINUE
122  ei2 = min(ei2,0.9999d0)
123 *
124 * Define unit vectors for inner eccentricity and angular momenta.
125  cosj = 0.0
126  sjsg = 0.0
127  DO 18 k = 1,3
128  ei(k) = ei(k)/sqrt(ei2)
129  hi(k) = a1(k)/sqrt(a12)
130  ho(k) = a2(k)/sqrt(a22)
131  cosj = cosj + hi(k)*ho(k)
132  sjsg = sjsg + ei(k)*ho(k)
133  18 CONTINUE
134 *
135 * Evaluate the expressions A & Z.
136  a = cosj*sqrt(1.0 - ei2)
137  z = (1.0 - ei2)*(2.0 - cosj**2) + 5.0*ei2*sjsg**2
138 *
139 * Obtain maximum inner eccentricity (Douglas Heggie, Sept. 1995).
140  z2 = z**2 + 25.0 + 16.0*a**4 - 10.0*z - 20.0*a**2 - 8.0*a**2*z
141  z2 = max(z2,0.0d0)
142  emax = (z + 1.0 - 4.0*a**2 + sqrt(z2))/6.0
143  emax = max(emax,0.0001d0)
144  emax = sqrt(emax)
145 *
146 * Check hierarchical stability condition for bound close pair (RB > a).
147  iterm = 0
148  alpha = 180.0*alpha/3.1415
149  IF (pmin.GT.pcrit.AND.semi.GT.0.0.AND.semi1.GT.0.0.AND.
150  & rb.GT.semi) THEN
151 * Delay termination until perturbation < 10^{-6} (includes CHAIN B-B).
152  IF (r3.GT.semi1.AND.gpert.LT.1.0d-06) THEN
153  iterm = -1
154  WRITE (6,20) namec(i1), namec(i2), namec(i3), ecc, emax,
155  & ecc1, semi, semi1, pmin, pcrit, alpha
156  20 FORMAT (' NEW HIARCH NM =',3i6,' E =',f6.3,' EX =',f7.4,
157  & ' E1 =',f6.3,' A =',1p,e8.1,
158  & ' A1 =',e8.1,' PM =',e9.2,
159  & ' PC =',e9.2,' IN =',0p,f7.1)
160  ri = sqrt(cm(1)**2 + cm(2)**2 + cm(3)**2)
161  q0 = m(i3)/mb
162  WRITE (81,30) timec, ri, namec(i3), q0, ecc, emax, ecc1,
163  & semi, semi1, pcrit/pmin, alpha
164  30 FORMAT (2f8.4,i6,f6.2,3f6.3,1p,2e10.2,0p,f5.2,f8.1)
165  CALL flush(81)
166  END IF
167 * Include termination test for wide triple system (exclude ECC1 > 0.9).
168  ELSE IF (pmin.GT.3.0*semi*(1.0 + ecc).AND.semi.GT.0.0.AND.
169  & ecc1.LT.0.9) THEN
170 * Wait for favourable configuration (R > SEMI, R3 > SEMI1 & RDOT3 > 0.
171  IF (rb.GT.semi.AND.r3.GT.semi1.AND.rdot3.GT.0.0) THEN
172  apo = semi*(1.0 + ecc)
173  WRITE (6,40) ecc, ecc1, alpha, rb, r3, pcrit, pmin, apo
174  40 FORMAT (' WIDE CHAIN E E1 IN RB R3 PC PM APO ',
175  & 2f7.3,f7.1,1p,5e10.2)
176  CALL flush(6)
177  iterm = -1
178  END IF
179 * Enforce termination of long-lived configuration using two limits.
180  ELSE IF (nstep1.GT.50000.AND.pmin.GT.0.9*pcrit.AND.
181  & r3.GT.semi1) THEN
182  iterm = -1
183  ELSE IF (nstep1.GT.500000.AND.r3.GT.semi1) THEN
184  iterm = -1
185  END IF
186 *
187  RETURN
188 *
189  END