Nbody6
 All Files Functions Variables
chaos.f
Go to the documentation of this file.
1  SUBROUTINE chaos(IPAIR,I1,I2,QPERI,ECC,IS,ZMU,RKS,SEMI1,ECC1,IDIS)
2 *
3 *
4 * Chaotic tidal interactions.
5 * ---------------------------
6 *
7 * Theory of Rosemary Mardling, Ap. J. XX, YYY, 1995.
8 * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
9 *
10  include 'common6.h'
11  common/modes/ eb0(ntmax),zj0(ntmax),ecrit(ntmax),ar(ntmax),
12  & br(ntmax),eosc(4,ntmax),edec(ntmax),tosc(ntmax),
13  & rp(ntmax),es(ntmax),cm(2,ntmax),iosc(ntmax),
14  & namec(ntmax)
15  REAL*8 de2(2),de3(2),ww(6),w(4),alf(4),tl(2),at(2),tdyn(2),
16  & wg(2),qg(2),wscale(2),qscale(2),a0(3),a2(3),eosc0(2),
17  & ql(2),td(2)
18  REAL*8 ran2
19  INTEGER is(2),kg(2)
20  CHARACTER*8 which1
21  SAVE kicks,ndec,eosc0,zjosc
22  DATA ww /2.119,3.113,8.175,3.742,4.953,9.413/
23  DATA eccm2 /0.00000399/
24 *
25 *
26 * Define c.m. index and search current names for chaos index.
27  idis = 0
28  i = n + ipair
29  ic = 0
30  DO 1 k = 1,nchaos
31  IF (namec(k).EQ.name(i)) ic = k
32  1 CONTINUE
33 *
34 * See whether new case could be an old SPIRAL exchanged in chain.
35  IF (ic.GT.0) THEN
36 * Remove tables and re-initialize on large epoch or standard binary.
37  IF (time - tosc(ic).GT.1.0.OR.kstar(i).EQ.0) THEN
38  ii = -i
39  CALL spiral(ii)
40  ic = 0
41  END IF
42  END IF
43 *
44 * Increase counter for new chaos and initialize index & variables.
45  IF (ic.EQ.0) THEN
46  nchaos = nchaos + 1
47  ic = nchaos
48  namec(ic) = name(i)
49  kicks = 0
50  ndec = 0
51  iosc(ic) = 0
52  kstar(i) = -1
53  tosc(ic) = time
54 * Ensure next location contains zero core mass (avoids confusion).
55  cm(1,nchaos+1) = 0.0
56  cm(2,nchaos+1) = 0.0
57  DO 5 k = 1,4
58  eosc(k,ic) = 0.0d0
59  5 CONTINUE
60  IF (kz(8).GT.3.AND.h(ipair).LT.0.0) THEN
61  CALL binev(ipair)
62  END IF
63  IF (nchaos.GT.ntmax) THEN
64  WRITE (6,6) name(i1), nchaos, ecc, qperi
65  6 FORMAT (' FATAL ERROR! CHAOS NM NCH E QP ',
66  & i6,i4,f8.4,1p,e9.1)
67  WRITE (6,7) (namec(k),k=1,nchaos)
68  7 FORMAT (' NAMEC ',12i6,(/,12i6))
69  stop
70  END IF
71  END IF
72 *
73 * Check termination for rare case of wide chaos (minimum 1000 calls).
74  qp = qperi/max(radius(i1),radius(i2))
75  IF (qp.GT.7.0.AND.ndec.GT.1000) THEN
76 * Activate spiral indicator and save time, pericentre & eccentricity.
77  kstar(i) = -2
78  tosc(ic) = time
79  rp(ic) = qperi
80  es(ic) = ecc
81  iosc(ic) = 2
82  nsp = nsp + 1
83  semi = -0.5*body(i)/h(ipair)
84  WRITE (6,8) name(i1), name(i2), ndec, list(1,i1), qp, ecc,
85  & semi
86  8 FORMAT (' WIDE CHAOS NAM NDEC NP QP/S E A ',
87  & 3i6,i4,f5.1,f8.4,1p,e10.2)
88  ndec = 0
89  kicks = 0
90  go to 34
91  END IF
92 *
93 * Define oscillation period (dimensionless time) and damping constants.
94  xn = 0.0
95  qd = 0.0
96  DO 10 k = 1,2
97  j = k + 2
98  ik = i1 + k - 1
99  tdyn(k) = radius(ik)*sqrt(radius(ik)/body(ik))
100 * Specify polytropic index for each star (n = 2, 3 or 3/2).
101  IF (kstar(ik).EQ.3.OR.kstar(ik).EQ.5.OR.
102  & kstar(ik).EQ.6.OR.kstar(ik).EQ.9) THEN
103  CALL giant(ipair,ik,wg,qg,wscale,qscale,xn,qd)
104  w(k) = wg(1)
105  w(j) = wg(2)
106  ql(k) = qd
107  kg(k) = 1
108  ELSE
109  ql(k) = 1.0d+04
110  kg(k) = 0
111  ip = 3
112  IF (kstar(ik).GE.3) ip = 2
113  IF (kstar(ik).EQ.4.OR.kstar(ik).EQ.7) ip = 3
114  IF (kstar(ik).EQ.8) ip = 3
115  IF (kstar(ik).EQ.0) ip = 1
116  w(k) = ww(ip)
117  w(j) = ww(ip+3)
118  END IF
119  alf(k) = 2.0*tdyn(k)/sqrt(w(k))
120  alf(j) = 3.0*tdyn(k)/sqrt(w(j))
121  tl(k) = twopi*tdyn(k)/sqrt(w(k))
122  10 CONTINUE
123 *
124 * Save initial eccentricity, binding energy & J0 the first time.
125  cj = zmu*sqrt(body(i))
126  IF (iosc(ic).EQ.0.OR.list(1,i1).GT.0) THEN
127  semi = -0.5*body(i)/h(ipair)
128  ecc0 = ecc
129  eb0(ic) = zmu*h(ipair)
130  zj0(ic) = cj*sqrt(qperi*(1.0 + ecc0))
131  edec(ic) = 0.0
132  iosc(ic) = 1
133 * Reset diagnostic indicator for CHAOS0 (output only first time).
134  IF (ndec.GT.0) idis = kstar(i)
135 *
136 * Initialize chaos boundary parameters (ECRIT, AR & BR).
137  CALL chaos0(qperi,ecc,eb0(ic),zj0(ic),body(i1),body(i2),
138  & radius(i1),radius(i2),w,ecrit(ic),ar(ic),br(ic),idis)
139 *
140 * Begin spiralling stage if chaos boundary has been crossed.
141  IF (idis.EQ.-1) THEN
142 * Include safety check on skipping rare case of hyperbolic orbit.
143  IF (ecc.GT.1.0) THEN
144  nchaos = nchaos - 1
145  kstar(i) = 0
146  go to 80
147  END IF
148 * Restrict tidal circularization to short time-scales (< 100 Myr).
149  icirc = -1
150  CALL tcirc(qperi,ecc,i1,i2,icirc,tc)
151  IF (tc.GT.100.0) THEN
152  nchaos = nchaos - 1
153  kstar(i) = 0
154  go to 80
155  END IF
156 * Skip initialization for difficult hierarchy and small EMAX.
157  IF (list(1,i1).EQ.1) THEN
158  icirc = -1
159  jcomp = list(2,i1)
160  CALL induce(ipair,emax,emin,icirc,tc,angle,tg,edav)
161  IF (emax.LT.0.9) THEN
162  nchaos = nchaos - 1
163  kstar(i) = 0
164  go to 80
165  END IF
166  END IF
167  iosc(ic) = 2
168  WRITE (6,12) ttot, name(i1), name(i2), kstar(i1),
169  & kstar(i2), list(1,i1), body(i1)*zmbar,
170  & body(i2)*zmbar, radius(i1)*su,
171  & radius(i2)*su, qperi, semi, ecc, xn, qd
172  12 FORMAT (' NEW SPIRAL T NM K* NP M1 M2 R* QP A E n Q ',
173  & f9.2,2i6,3i3,2f5.1,2f6.1,
174  & 1p,2e10.2,0p,f7.3,f5.1,f7.1)
175 * Activate spiral indicator and save time, pericentre & eccentricity.
176  kstar(i) = -2
177  tosc(ic) = time
178  rp(ic) = qperi
179  es(ic) = ecc0
180  nsp = nsp + 1
181 * Rectify the KS solution on transition to standard circularization.
182  CALL ksrect(ipair)
183 * Initialize perturbed KS (small STEP after integration in KSPERI).
184  IF (list(1,i1).GT.0) THEN
185  imod = 1
186  CALL kspoly(ipair,imod)
187  kslow(ipair) = 1
188  END IF
189  go to 34
190  END IF
191 *
192 * Reduce chaos index on disruption if current case is last.
193  IF (idis.GT.0.AND.ic.EQ.nchaos) THEN
194  nchaos = nchaos - 1
195  go to 80
196  END IF
197 *
198 * Print NEW CHAOS/CAPTURE the first time.
199  IF (ndec.EQ.0) THEN
200  ncha = ncha + 1
201  which1 = ' CHAOS '
202  IF (semi.LT.0.0) which1 = ' CAPTURE'
203  ga = gamma(ipair)*(rmin/rks)**3
204  WRITE (6,15) which1, ttot, name(i1), name(i2), kstar(i1),
205  & kstar(i2), list(1,i1), body(i1)*zmbar,
206  & body(i2)*zmbar, radius(i1)*su,
207  & radius(i2)*su, qperi, semi, ecc, ga, xn, qd
208  15 FORMAT (' NEW',a8,' T NM K* NP M1 M2 R* QP A E G n Q ',
209  & f9.2,2i6,3i3,2f5.1,2f6.1,1p,2e10.2,
210  & 0p,f9.5,f7.3,f5.1,f7.1)
211  END IF
212  END IF
213 *
214 * Obtain energy dissipation from separate modes.
215  CALL tides2(qperi,body(i1),body(i2),radius(i1),radius(i2),is,ecc,
216  & kg,wscale,qscale,de2,de3)
217 *
218 * Employ an emergency procedure for zero energy change (TIDES2 bug).
219  IF (de2(1).EQ.0.0d0.OR.de2(2).EQ.0.0d0) THEN
220  eb = bodym*h(ipair)
221  de2(1) = -0.001*eb*ran2(idum)
222  de2(2) = -0.001*eb*ran2(idum)
223  END IF
224 *
225 * Evaluate time-scale for Kochanek-type damping (linear & non-linear).
226  eosc0(1) = eosc(1,ic)**2 + eosc(2,ic)**2
227  eosc0(2) = eosc(3,ic)**2 + eosc(4,ic)**2
228  eosc0(1) = sqrt(eosc0(1))
229  eosc0(2) = sqrt(eosc0(2))
230 * Adopt non-linear dissipation time scale of Kumar & Goodman 1995.
231  qnl1 = ql(1)/max(sqrt(eosc0(1)),0.00001d0)
232  qnl2 = ql(2)/max(sqrt(eosc0(2)),0.00001d0)
233  semi = -0.5*body(i)/h(ipair)
234  semi = abs(semi)
235  tk = twopi*semi*sqrt(semi/body(i))
236  td(1) = (1.0/ql(1) + 1.0/qnl1)*tk
237  td(2) = (1.0/ql(2) + 1.0/qnl2)*tk
238 *
239  rm = max(radius(i1),radius(i2))
240  IF (ndec.GT.10000.AND.qperi.GT.4.0*rm) THEN
241  DO 18 k = 1,2
242  de2(k) = 1000.0*de2(k)
243  de3(k) = 1000.0*de3(k)
244  18 CONTINUE
245  END IF
246 * Sum old and new oscillation energies for all modes.
247  e20 = 0.0
248  e30 = 0.0
249  e2t = 0.0
250  e3t = 0.0
251  zjosc = 0.0
252  DO 20 k = 1,2
253  j = k + 2
254  at(k) = -td(k)/tl(k)
255  e20 = e20 + eosc(k,ic)
256  e30 = e30 + eosc(j,ic)
257  eosc(k,ic) = eosc(k,ic)*exp(at(k))
258  eosc(j,ic) = eosc(j,ic)*exp(at(k))
259  edec(ic) = edec(ic) - (eosc(k,ic) + eosc(j,ic))
260  IF (iosc(ic).EQ.-1) THEN
261  delta = 0.5*twopi
262  ELSE
263  delta = twopi*ran2(idum1)
264  END IF
265  eosc(k,ic) = eosc(k,ic) +
266  & 2.0*sqrt(eosc(k,ic)*de2(k))*cos(delta) + de2(k)
267  IF (iosc(ic).EQ.-1) THEN
268  IF (k.EQ.2) iosc(ic) = 1
269  ELSE
270  delta = twopi*ran2(idum1)
271  END IF
272  eosc(j,ic) = eosc(j,ic) +
273  & 2.0*sqrt(eosc(j,ic)*de3(k))*cos(delta) + de3(k)
274 * Ensure that oscillation energies are not negative.
275  eosc(k,ic) = max(eosc(k,ic),0.0d0)
276  eosc(j,ic) = max(eosc(j,ic),0.0d0)
277  e2t = e2t + eosc(k,ic)
278  e3t = e3t + eosc(j,ic)
279  zjosc = zjosc + alf(k)*eosc(k,ic) + alf(j)*eosc(j,ic)
280  20 CONTINUE
281 *
282 * Specify change in oscillation energy and sum decayed energy.
283 * DET = (E2T - E20) + (E3T - E30)
284  edec(ic) = edec(ic) + e20 + e30
285  tosc(ic) = time
286 *
287 * Set new binding energy & semi-major axis.
288  hi = h(ipair)
289  hnew = (eb0(ic) - edec(ic) - (e2t + e3t))/zmu
290  semi1 = -0.5*body(i)/hnew
291 *
292 * Calculate the new eccentricity.
293  efac = (zj0(ic) - zjosc)/cj
294  ecc2 = 1.0 - efac**2/semi1
295  ecc2 = max(ecc2,eccm2)
296  ecc1 = sqrt(ecc2)
297  peri1 = semi1*(1.0d0 - ecc1)
298 *
299 * Switch off chaos indicator on transition to hyperbolic orbit.
300  IF (hnew.GT.0.0) THEN
301  kstar(i) = 0
302 * Reduce index if current case is last (otherwise updated in KSTERM).
303  IF (ic.EQ.nchaos) THEN
304  nchaos = nchaos - 1
305  END IF
306  WRITE (6,25) ipair, ndec, kicks, ecc, ecc1, semi1, qperi,
307  & hnew - hi
308  25 FORMAT (' TERMINATED CHAOS IPAIR NDEC KICK E E1 A1 QP DH ',
309  & 3i4,2f9.5,1p,3e10.2)
310  go to 80
311  END IF
312 *
313 * Update total energy loss due to change in binding energy.
314  ndec = ndec + 1
315  h(ipair) = hnew
316  deb = zmu*(hi - h(ipair))
317  ecoll = ecoll + deb
318  egrav = egrav + deb
319  e(10) = e(10) + deb
320 *
321 * Check energy or eccentricity criterion for chaotic case.
322  IF (kz(27).EQ.2.AND.(iosc(ic).EQ.1.OR.iosc(ic).EQ.-1)) THEN
323  IF (edec(ic).GT.-(ecrit(ic) - eb0(ic))) THEN
324  iosc(ic) = 2
325  WRITE (6,30) ttot, name(i1), name(i2), ndec, kicks,
326  & ecc1, semi1, ecrit(ic), edec(ic)
327  30 FORMAT (' END CHAOS T NM NDEC KICKS E A ECRIT EDEC ',
328  & f9.2,4i6,f8.4,1p,3e10.2)
329 * Activate spiral indicator and save time, pericentre & eccentricity.
330  kstar(i) = -2
331  tosc(ic) = time
332  rp(ic) = peri1
333  es(ic) = ecc1
334  nsp = nsp + 1
335  ndec = 0
336  kicks = 0
337  ELSE
338  eps = (2.0*ecrit(ic) - eb0(ic) + edec(ic))/(zmu*body(i))
339  eccm = (1.0 + eps*br(ic))/(1.0 - eps*ar(ic))
340  IF (ecc1.LT.eccm) THEN
341  iosc(ic) = -1
342  kicks = kicks + 1
343  IF (kicks.LT.3) THEN
344  WRITE (6,32) ndec, eps, eccm, ecc1
345  32 FORMAT (' NEW CHAOS KICK NDEC EPS ECCM E ',
346  & i5,1p,e10.2,0p,2f8.4)
347  END IF
348  END IF
349  END IF
350  END IF
351 *
352 * Check for hierarchical configuration on first call.
353  34 IF (ndec.LE.1.AND.semi.GT.0.0.AND.semi.LT.2.0*rmin) THEN
354  np1 = list(1,i1) + 1
355  DO 39 l = 2,np1
356  j = list(l,i1)
357  rij2 = 0.0
358  vij2 = 0.0
359  rdot = 0.0
360  a12 = 0.0
361  a22 = 0.0
362  a1a2 = 0.0
363  DO 35 k = 1,3
364  rij2 = rij2 + (x(k,i) - x(k,j))**2
365  vij2 = vij2 + (xdot(k,i) - xdot(k,j))**2
366  rdot = rdot + (x(k,i) - x(k,j))*(xdot(k,i) -xdot(k,j))
367  k1 = k + 1
368  IF (k1.GT.3) k1 = 1
369  k2 = k1 + 1
370  IF (k2.GT.3) k2 = 1
371  a0(k) = (x(k1,i1)-x(k1,i2))*(xdot(k2,i1)-xdot(k2,i2))
372  & - (x(k2,i1)-x(k2,i2))*(xdot(k1,i1)-xdot(k1,i2))
373  a2(k) = (x(k1,j) - x(k1,i))*(xdot(k2,j) - xdot(k2,i))
374  & - (x(k2,j) - x(k2,i))*(xdot(k1,j) - xdot(k1,i))
375  a12 = a12 + a0(k)**2
376  a22 = a22 + a2(k)**2
377  a1a2 = a1a2 + a0(k)*a2(k)
378  35 CONTINUE
379  rip = sqrt(rij2)
380  a1 = 2.0/rip - vij2/(body(i) + body(j))
381  a1 = 1.0/a1
382 * Include impact parameter test for recoil check.
383  a4 = rdot**2/(a1*(body(i) + body(j)))
384  eccp = sqrt((1.0d0 - rip/a1)**2 + a4)
385  pmin = a1*(1.0d0 - eccp)
386  IF (pmin.LT.3.0*semi) THEN
387  tm = min(tev(i1),tev(i2)) - time
388  WRITE (6,36) ipair, name(j), eccp, pmin/semi,
389  & rdot/rip, semi, a1, rip, tm
390  36 FORMAT (' RECOIL: KS NMJ E1 PM/A RD A0 A1 RP TM ',
391  & i4,i6,f7.3,f5.1,f6.1,1p,4e10.2)
392  END IF
393 * Accept semi-major axis ratio below 25.
394  IF (1.0/a1.GT.0.04/semi.AND.idis.LE.0) THEN
395  ra = semi*(1.0 + ecc)
396  sr = pmin/ra
397  ga = 2.0*body(j)*(ra/pmin)**3/body(i)
398 * Determine inclination (8 bins of 22.5 degrees).
399  fac = a1a2/sqrt(a12*a22)
400  fac = acos(fac)
401  in = 1 + fac*360.0/(twopi*22.5)
402  WRITE (6,38) ipair, name(j), h(ipair), semi, a1,
403  & pmin, ga, eccp, sr, in
404  38 FORMAT (' HIERARCHY: KS NMJ H A0 A1 RP GA E1 SR ',
405  & 'IN',i5,i6,f7.0,1p,3e9.1,0p,2f6.2,f6.1,i3)
406  END IF
407  39 CONTINUE
408  END IF
409 *
410 * Check for terminated or escaped chaotic binaries first time.
411  IF (ndec.EQ.1.AND.nchaos.GT.1) THEN
412  j = 1
413 * See whether case #J indicates current or escaped/disrupted KS binary.
414  40 DO 45 jpair = 1,npairs
415  IF (namec(j).EQ.name(n+jpair)) THEN
416 * Update #J if KSTAR > 0, otherwise consider next member.
417  IF (kstar(n+jpair).GT.0) THEN
418  go to 50
419  ELSE
420  go to 70
421  END IF
422  END IF
423  45 CONTINUE
424 *
425 * Skip during mergers or multiple regularizations (KSTAR not visible).
426  IF (nmerge.GT.0.OR.nsub.GT.0) go to 70
427 *
428 * Update chaos variables for #IC and any disrupted or escaped binaries.
429  50 nchaos = nchaos - 1
430  DO 60 l = j,nchaos
431  l1 = l + 1
432  DO 55 k = 1,4
433  eosc(k,l) = eosc(k,l1)
434  55 CONTINUE
435  eb0(l) = eb0(l1)
436  zj0(l) = zj0(l1)
437  ecrit(l) = ecrit(l1)
438  ar(l) = ar(l1)
439  br(l) = br(l1)
440  edec(l) = edec(l1)
441  tosc(l) = tosc(l1)
442  rp(l) = rp(l1)
443  es(l) = es(l1)
444  cm(1,l) = cm(1,l1)
445  cm(2,l) = cm(2,l1)
446  iosc(l) = iosc(l1)
447  namec(l) = namec(l1)
448  60 CONTINUE
449 * Consider the same location again after each removal (J <= NCHAOS).
450  j = j - 1
451  70 j = j + 1
452  IF (j.LE.nchaos) go to 40
453  END IF
454 *
455 * Check optional binary diagnostics on transition from chaotic state.
456  80 IF (kz(8).GT.3.AND.kstar(i).NE.-1) THEN
457 * Skip output for unperturbed case (CALL KSTIDE from UNPERT).
458  IF (list(1,i1).GT.0.AND.h(ipair).LT.0.0) THEN
459  CALL binev(ipair)
460  END IF
461  END IF
462 *
463  RETURN
464 *
465  END