Nbody6
 All Files Functions Variables
chaos2.f
Go to the documentation of this file.
1  SUBROUTINE chaos2(I1,I2,ECC,HI,IS,BODYI,ZMU,RSTAR,SEMI1,ECC1,DH,
2  & idis,kstari)
3 *
4 *
5 * Chain chaotic treatment of tidal capture.
6 * -----------------------------------------
7 *
8 * Theory of Rosemary Mardling, Ap. J. XX, YYY, 1995.
9 * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
10 *
11  include 'common6.h'
12  parameter(nmx=10,nmx2=2*nmx,nmx3=3*nmx,nmx4=4*nmx,
13  & nmx8=8*nmx,nmxm=nmx*(nmx-1)/2)
14  REAL*8 m,mass,mc,mij,mkk,rstar(2)
15  common/chain1/ xch(nmx3),vch(nmx3),m(nmx),
16  & zz(nmx3),wc(nmx3),mc(nmx),
17  & xi(nmx3),pi(nmx3),mass,rinv(nmxm),rsum,mkk(nmx),
18  & mij(nmx,nmx),tkk(nmx),tk1(nmx),iname(nmx),nn
19 * Note renaming of NAMEC to NAMEX in COMMON/CHREG/ to avoid conflict.
20  common/chreg/ timec,tmax,rmaxc,cm(10),namex(6),nstep1,kz27,kz30
21  common/ccoll2/ qk(nmx4),pk(nmx4),rij(nmx,nmx),SIZE(nmx),vstar1,
22  & ecoll1,rcoll,qperi,istar(nmx),icoll,isync,ndiss1
23  common/modes/ eb0(ntmax),zj0(ntmax),ecrit(ntmax),ar(ntmax),
24  & br(ntmax),eosc(4,ntmax),edec(ntmax),tosc(ntmax),
25  & rp(ntmax),es(ntmax),cx(2,ntmax),iosc(ntmax),
26  & namec(ntmax)
27  REAL*8 de2(2),de3(2),ww(6),w(4),alf(4),tl(2),at(2),tdyn(2),wg(2),
28  & qg(2),wscale(2),qscale(2),eosc0(2),ql(2),td(2)
29  REAL*8 ran2
30  INTEGER is(2),kg(2)
31  CHARACTER*8 which1
32  SAVE kicks,ndec
33  DATA ww /2.119,3.113,8.175,3.742,4.953,9.413/
34  DATA eccm2 /0.00000399/
35 *
36 *
37 * Skip treatment for circular or hyperbolic stage (ISYNC = 1 or -1).
38  IF (isync.NE.0) go to 100
39 *
40  IF (ecc.LT.0.0021) THEN
41  kstari = 10
42  isync = 1
43  go to 100
44  END IF
45 *
46 * See whether binary c.m. can be identified from the chaos table.
47  ic = 0
48  nam1 = nzero + namex(i1)
49  nam2 = nzero + namex(i2)
50  DO 1 k = 1,nchaos
51  IF (namec(k).EQ.nam1.OR.namec(k).EQ.nam2) ic = k
52  1 CONTINUE
53 *
54 * Skip existing spiral (long t_circ for both original and new case).
55  IF (ic.GT.0) THEN
56  k = 0
57  IF (namex(i1) + namex(i2).EQ.ksave(2)) k = 1
58  IF (namex(i1) + namex(i2).EQ.ksave(4)) k = 3
59  IF (k.GT.0) THEN
60  kstari = ksave(k)
61  IF (kstari.EQ.-2) THEN
62  isync = 1
63  go to 100
64  END IF
65  END IF
66  END IF
67 *
68 * Increase counter for new chaos and initialize index & variables.
69  IF (ic.EQ.0) THEN
70  nchaos = nchaos + 1
71  ic = nchaos
72  namec(ic) = nam1
73  ndec = 0
74  iosc(ic) = 0
75  kstari = -1
76  tosc(ic) = time + timec
77  DO 5 k = 1,4
78  eosc(k,ic) = 0.0d0
79  5 CONTINUE
80  IF (nchaos.GT.ntmax) THEN
81  WRITE (6,6) name(i1), nchaos, ecc
82  6 FORMAT (' FATAL ERROR! CHAOS2 NM NCH E ',
83  & i6,i4,f8.4)
84  stop
85  END IF
86  END IF
87 *
88 * Define oscillation period (dimensionless time) and damping constants.
89  DO 10 k = 1,2
90  j = k + 2
91  ik = i1
92  IF (k.EQ.2) ik = i2
93  tdyn(k) = sqrt(rstar(k)**3/m(ik))
94 * Specify polytropic index for each star (n = 2, 3 or 3/2).
95  IF (is(k).EQ.3.OR.is(k).EQ.5) THEN
96 * Set fudged KS index = 0 to identify second particle in routine GIANT.
97  CALL giant2(k,ik,wg,qg,wscale,qscale,xn,qd)
98  w(k) = wg(1)
99  w(j) = wg(2)
100  ql(k) = qd
101  kg(k) = 1
102  ELSE
103  ql(k) = 1.0d+04
104  kg(k) = 0
105  ip = 3
106  IF (is(k).GE.3) ip = 2
107  IF (is(k).EQ.4.OR.is(k).EQ.6) ip = 3
108  IF (is(k).EQ.0) ip = 1
109  w(k) = ww(ip)
110  w(j) = ww(ip+3)
111  END IF
112  alf(k) = 2.0*tdyn(k)/sqrt(w(k))
113  alf(j) = 3.0*tdyn(k)/sqrt(w(j))
114  tl(k) = twopi*tdyn(k)/sqrt(w(k))
115  10 CONTINUE
116 *
117 * Save initial eccentricity, binding energy & J0.
118  semi = -0.5*bodyi/hi
119  cj = zmu*sqrt(bodyi)
120  IF (iosc(ic).EQ.0) THEN
121  ecc0 = ecc
122  eb0(ic) = zmu*hi
123  zj0(ic) = cj*sqrt(qperi*(1.0 + ecc0))
124  edec(ic) = 0.0
125  iosc(ic) = 1
126  tosc(ic) = time + timec
127  rp(ic) = qperi
128  es(ic) = ecc0
129  kicks = 0
130 * Initialize chaos boundary parameters (ECRIT, AR & BR).
131  CALL chaos0(qperi,ecc,eb0(ic),zj0(ic),m(i1),m(i2),
132  & rstar(1),rstar(2),w,ecrit(ic),ar(ic),br(ic),idis)
133 *
134 * Check whether chaotic stage should be replaced by spiralling.
135  IF (idis.EQ.-1) THEN
136  iosc(ic) = 2
137  WRITE (6,12) is(1), is(2), m(i1), m(i2), rstar(1),
138  & rstar(2), qperi, semi, ecc
139  12 FORMAT (' CHAIN SPIRAL K* M1 M2 R* QP A E ',
140  & 2i4,1p,6e10.2,0p,f8.4)
141 * Activate spiral indicator and save time, pericentre & eccentricity.
142  kstari = -2
143  isync = 1
144  tosc(ic) = time + timec
145  rp(ic) = qperi
146  es(ic) = ecc0
147  namc = 0
148  DO 14 j = 1, nchaos
149  IF (namec(j).EQ.nam1.OR.namec(j).EQ.nam2) THEN
150  namc = namec(j)
151  END IF
152  14 CONTINUE
153  ic = 0
154 * Obtain diagnostic output after strong interaction.
155  CALL recoil(3)
156  go to 90
157  END IF
158 *
159  ncha = ncha + 1
160  which1 = ' CHAOS '
161  WRITE (6,15) which1, namex(i1), namex(i2), is(1), is(2),
162  & m(i1)*zmbar, m(i2)*zmbar, rstar(1), rstar(2),
163  & qperi, semi, ecc
164  15 FORMAT (' CHAIN',a8,' NAM K* M1 M2 R* QP A E ',
165  & 2i6,2i4,2f5.1,1p,4e10.2,0p,f8.4)
166 *
167 * Obtain diagnostic output after strong interaction.
168  CALL recoil(3)
169 * Exit on positive indicator (collision is too difficult here).
170  IF (idis.GT.1) go to 90
171  END IF
172 *
173 * Obtain energy dissipation from separate modes.
174  CALL tides2(qperi,m(i1),m(i2),rstar(1),rstar(2),is,ecc,kg,
175  & wscale,qscale,de2,de3)
176 *
177 * Evaluate time-scale for Kochanek-type damping (linear & non-linear).
178  eosc0(1) = eosc(1,ic)**2 + eosc(2,ic)**2
179  eosc0(2) = eosc(3,ic)**2 + eosc(4,ic)**2
180  eosc0(1) = sqrt(eosc0(1))
181  eosc0(2) = sqrt(eosc0(2))
182 * Adopt non-linear dissipation time scale of Kumar & Goodman 1995.
183  qnl1 = ql(1)/max(sqrt(eosc0(1)),0.00001d0)
184  qnl2 = ql(2)/max(sqrt(eosc0(2)),0.00001d0)
185  semi = abs(semi)
186  tk = twopi*semi*sqrt(semi/bodyi)
187  td(1) = (1.0/ql(1) + 1.0/qnl1)*tk
188  td(2) = (1.0/ql(2) + 1.0/qnl2)*tk
189 *
190 * Include check on wide interaction near the boundary.
191  rm = max(rstar(1),rstar(2))
192  IF (ndec.GT.10000.AND.qperi.GT.4.0*rm) THEN
193  DO 18 k = 1,2
194  de2(k) = 1000.0*de2(k)
195  de3(k) = 1000.0*de3(k)
196  18 CONTINUE
197  END IF
198 *
199 * Sum old and new oscillation energies for all modes.
200  e20 = 0.0
201  e30 = 0.0
202  e2t = 0.0
203  e3t = 0.0
204  zjosc = 0.0
205  DO 20 k = 1,2
206  j = k + 2
207  at(k) = -td(k)/tl(k)
208  e20 = e20 + eosc(k,ic)
209  e30 = e30 + eosc(j,ic)
210  eosc(k,ic) = eosc(k,ic)*exp(at(k))
211  eosc(j,ic) = eosc(j,ic)*exp(at(k))
212  edec(ic) = edec(ic) - (eosc(k,ic) + eosc(j,ic))
213  IF (iosc(ic).EQ.-1) THEN
214  delta = 0.5*twopi
215  ELSE
216  delta = twopi*ran2(idum1)
217  END IF
218  eosc(k,ic) = eosc(k,ic) +
219  & 2.0*sqrt(eosc(k,ic)*de2(k))*cos(delta) + de2(k)
220  IF (iosc(ic).EQ.-1) THEN
221  IF (k.EQ.2) iosc(ic) = 1
222  ELSE
223  delta = twopi*ran2(idum1)
224  END IF
225  eosc(j,ic) = eosc(j,ic) +
226  & 2.0*sqrt(eosc(j,ic)*de3(k))*cos(delta) + de3(k)
227 * Ensure that oscillation energies are not negative.
228  eosc(k,ic) = max(eosc(k,ic),0.0d0)
229  eosc(j,ic) = max(eosc(j,ic),0.0d0)
230  e2t = e2t + eosc(k,ic)
231  e3t = e3t + eosc(j,ic)
232  zjosc = zjosc + alf(k)*eosc(k,ic) + alf(j)*eosc(j,ic)
233  20 CONTINUE
234 *
235 * Specify change in oscillation energy and sum decayed energy.
236 * DET = (E2T - E20) + (E3T - E30)
237  edec(ic) = edec(ic) + e20 + e30
238 *
239 * Set new binding energy & semi-major axis.
240  hnew = (eb0(ic) - edec(ic) - (e2t + e3t))/zmu
241  dh = hnew - hi
242 * Note that the net change in EGRAV is taken care of at termination.
243  semi1 = -0.5*bodyi/hnew
244 *
245 * Calculate the new eccentricity and specify pericentre.
246  efac = (zj0(ic) - zjosc)/cj
247  ecc2 = 1.0 - efac**2/semi1
248  ecc2 = max(ecc2,eccm2)
249  ecc1 = sqrt(ecc2)
250  peri1 = semi1*(1.0d0 - ecc1)
251 *
252 * Switch off chaos indicator on transition to hyperbolic orbit.
253  IF (hnew.GT.0.0) THEN
254  kstari = 0
255  isync = -1
256 * Reduce index if current case is last (otherwise updated later).
257  IF (ic.EQ.nchaos) THEN
258  nchaos = nchaos - 1
259  END IF
260  ic = 0
261  WRITE (6,25) namex(i1), namex(i2), ndec, kicks, ecc, ecc1,
262  & semi1, qperi
263  25 FORMAT (' TERMINATED CHAOS NAM NDEC KICK E E1 A1 QP ',
264  & 2i6,2i4,2f8.4,1p,2e10.2)
265  ELSE IF (ecc.GT.1.0) THEN
266  vinf = sqrt(2.0*hi)*vstar1
267  p = days*semi1*sqrt(semi1/bodyi)
268  WRITE (6,26) time+toff, namex(i1), namex(i2), ecc, ecc1,
269  & vinf, semi1, p
270  26 FORMAT (' CHAIN CAPTURE T NAM E E1 VINF A1 P ',
271  & f10.3,2i6,f9.5,f8.4,f5.1,1p,2e9.1)
272  END IF
273 *
274 * Check energy or eccentricity criterion for chaotic case.
275  IF (iosc(ic).EQ.1.OR.iosc(ic).EQ.-1) THEN
276  IF (edec(ic).GT.-(ecrit(ic) - eb0(ic))) THEN
277  iosc(ic) = 2
278  WRITE (6,30) namex(i1), namex(i2), ndec, kicks, ecc1,
279  & semi1, ecrit(ic), edec(ic)
280  30 FORMAT (' END CHAOS NAM NDEC KICKS E A ECRIT EDEC ',
281  & 2i6,2i5,f8.4,1p,3e10.2)
282 * Activate spiral indicator and save time, pericentre & eccentricity.
283  kstari = -2
284  isync = 1
285  tosc(ic) = time + timec
286  rp(ic) = peri1
287  es(ic) = ecc1
288  ic = 0
289  ELSE
290  eps = (2.0*ecrit(ic) - eb0(ic) + edec(ic))/(zmu*bodyi)
291  eccm = (1.0 + eps*br(ic))/(1.0 - eps*ar(ic))
292  IF (ecc1.LT.eccm) THEN
293  iosc(ic) = -1
294  kicks = kicks + 1
295  IF (kicks.LT.3) THEN
296  WRITE (6,32) ndec, eps, eccm, ecc1
297  32 FORMAT (' NEW KICK NDEC EPS ECCM E ',
298  & i5,1p,e10.2,0p,2f8.4)
299  END IF
300  END IF
301  END IF
302  END IF
303 *
304 * Reduce chaos index on disruption if current case is last.
305 * IF (IDIS.GT.0.AND.IC.EQ.NCHAOS) THEN
306 * NCHAOS = NCHAOS - 1
307 * END IF
308  ndec = ndec + 1
309 *
310 * Skip energy change for small ECC1 and reduce NCHAOS if new SPIRAL.
311  IF (ecc1.LT.0.0021) THEN
312  kstari = 10
313  IF (ic.EQ.0) nchaos = nchaos - 1
314  ic = 0
315  END IF
316 *
317 * Update dissipation indicator in case of change.
318  90 k = 0
319  IF (namex(i1) + namex(i2).EQ.ksave(2)) k = 1
320  IF (namex(i1) + namex(i2).EQ.ksave(4)) k = 3
321  IF (k.GT.0) THEN
322  ksave(k) = kstari
323  ELSE IF (kstari.LT.0) THEN
324 * Include possibility of exchange with KSTAR < 0.
325  l = 1
326  IF (nn.GT.3) l = 3
327  ksave(l) = kstari
328  ksave(l+1) = namex(i1) + namex(i2)
329  END IF
330 *
331  100 RETURN
332 *
333  END