Nbody6
 All Files Functions Variables
coal.f
Go to the documentation of this file.
1  SUBROUTINE coal(IPAIR,KW1,KW2,MASS)
2 *
3 *
4 * Coalescence of Roche/CE binary.
5 * -------------------------------
6 *
7  include 'common6.h'
8  CHARACTER*8 which1
9  REAL*8 cm(6)
10  REAL*8 mass(2)
11  LOGICAL first
12  SAVE first
13  DATA first /.true./
14 *
15 *
16 * Distinguish between KS and chain regularization.
17  IF (ipair.GT.0) THEN
18 * Define discrete time for prediction & new polynomials (T <= TBLOCK).
19  i = n + ipair
20  dt = 0.1d0*step(i)
21  IF (dt.GT.2.4e-11) THEN
22  time2 = time - tprev
23  CALL stepk(dt,dtn)
24  time = tprev + int((time2 + dt)/dtn)*dtn
25  time = min(tblock,time)
26  ELSE
27  time = min(t0(i) + step(i),tblock)
28  END IF
29 *
30 * Set zero energy (EB correction done in routines EXPEL & EXPEL2).
31  eb = 0.d0
32  rcoll = r(ipair)
33  dmin2 = min(dmin2,rcoll)
34  vinf = 0.0
35 *
36 * Define indicator for different cases, including hyperbolic KS.
37  IF ((kstar(i).LE.10.AND.iqcoll.NE.0).OR.iqcoll.EQ.-2) THEN
38  which1 = ' BINARY '
39  iqcoll = 0
40  eb1 = body(2*ipair-1)*body(2*ipair)*h(ipair)/body(i)
41  ELSE
42  which1 = ' ROCHE '
43 * Save energy for correction (Roche COAL but not via CMBODY & EXPEL).
44  IF (iqcoll.EQ.0) THEN
45  eb = body(2*ipair-1)*body(2*ipair)*h(ipair)/body(i)
46  eb1 = eb
47  egrav = egrav + eb
48  END IF
49  iqcoll = 3
50  END IF
51  IF (h(ipair).GT.0.0) THEN
52  which1 = ' HYPERB '
53  nhyp = nhyp + 1
54  iqcoll = -1
55  vinf = sqrt(2.0*h(ipair))*vstar
56  eb1 = body(2*ipair-1)*body(2*ipair)*h(ipair)/body(i)
57  END IF
58 *
59 * Check optional diagnostics for binary evolution.
60  IF (kz(8).GT.3) THEN
61  CALL binev(ipair)
62  END IF
63 *
64 * Remove any circularized KS binary from the CHAOS/SYNCH table.
65  IF (kstar(i).GE.10.AND.nchaos.GT.0) THEN
66  ii = -i
67  CALL spiral(ii)
68  END IF
69 *
70 * Terminate KS pair and set relevant indices for collision treatment.
71  iphase = 9
72  kspair = ipair
73  t0(2*ipair-1) = time
74  CALL ksterm
75  i1 = 2*npairs + 1
76  i2 = i1 + 1
77  ELSE
78 * Copy dominant indices, two-body separation and binding energy.
79  i1 = jlist(1)
80  i2 = jlist(2)
81  rcoll = dminc
82  eb = 0.0
83  vinf = 0.d0
84  iqcoll = 5
85  which1 = ' CHAIN '
86 * Note that new chain TIME already quantized in routine CHTERM.
87  END IF
88 *
89 * Define global c.m. coordinates & velocities from body #I1 & I2.
90  icomp = i1
91  zm = body(i1) + body(i2)
92  DO 5 k = 1,3
93  cm(k) = (body(i1)*x(k,i1) + body(i2)*x(k,i2))/zm
94  cm(k+3) = (body(i1)*xdot(k,i1) + body(i2)*xdot(k,i2))/zm
95  5 CONTINUE
96 *
97 * Form central distance (scaled by RC), central distance and period.
98  ri2 = 0.d0
99  rij2 = 0.d0
100  vij2 = 0.d0
101  DO 10 k = 1,3
102  ri2 = ri2 + (x(k,i1) - rdens(k))**2
103  rij2 = rij2 + (x(k,i1) - x(k,i2))**2
104  vij2 = vij2 + (xdot(k,i1) - xdot(k,i2))**2
105  10 CONTINUE
106  ri = sqrt(ri2)
107  rij = sqrt(rij2)
108  semi = 2.d0/rij - vij2/zm
109  semi = 1.d0/semi
110  tk = days*semi*sqrt(abs(semi)/zm)
111  ecc = max(1.0 - rcoll/semi,0.001d0)
112  IF (iqcoll.EQ.5) THEN
113  eb1 = -0.5*body(i1)*body(i2)/semi
114  END IF
115 *
116 * Choose appropriate reference body containing #ICH neighbour list.
117  IF (ipair.GT.0) THEN
118  icm = i1
119  ELSE
120  icm = jlist(3)
121  IF (list(1,i1).GT.list(1,icm)) icm = i1
122  IF (list(1,i2).GT.list(1,icm)) icm = i2
123 * Include possible 4th chain member just in case.
124  IF (jlist(4).GT.0) THEN
125  i4 = jlist(4)
126  IF (list(1,i4).GT.list(1,icm)) icm = 4
127  END IF
128  END IF
129 *
130 * Form perturber list for corrections and add chain #I1 to NB lists.
131  nnb = list(1,icm)
132  l2 = 1
133  jmin = n + 1
134  rij2 = 1.0d+10
135 * Copy neighbour list of former c.m. (less any #I2) to JPERT.
136  DO 15 l = 1,nnb
137  jpert(l) = list(l+1,icm)
138  IF (jpert(l).EQ.i2) THEN !unlikely condition but no harm.
139  l2 = l
140  jpert(l) = n
141  IF (i2.EQ.n) jpert(l) = n - 1
142  ELSE
143 * Determine closest external member for possible KS.
144  jj = jpert(l)
145  ri2 = 0.d0
146  DO 12 k = 1,3
147  ri2 = ri2 + (cm(k) - x(k,jj))**2
148  12 CONTINUE
149  IF (ri2.LT.rij2) THEN
150  jmin = jj
151  rij2 = ri2
152  ENDIF
153  END IF
154  15 CONTINUE
155 *
156 * Restore at least chain collider mass in local neighbour lists.
157  IF (ipair.LE.0) THEN
158  nnb = nnb + 1
159  jpert(nnb) = icm
160 * Include the case of old c.m. coming from #I1 or #I2.
161  IF (icm.EQ.i1.OR.icm.EQ.i2) jpert(nnb) = jlist(3)
162  jlist(1) = i1
163  IF (icm.EQ.i1) jlist(1) = jlist(3)
164  CALL nbrest(icm,1,nnb)
165  END IF
166 *
167 * Evaluate potential energy with respect to colliding bodies.
168  jlist(1) = i1
169  jlist(2) = i2
170  CALL nbpot(2,nnb,pot1)
171 *
172 * Specify new mass from sum and initialize zero mass ghost in #I2.
173  zmnew = (mass(1) + mass(2))/zmbar
174  dm = zm - zmnew
175  IF (dm.LT.1.0d-10) dm = 0.d0
176 * Delay inclusion of any mass loss until after energy correction.
177  zm1 = body(i1)*zmbar
178  zm2 = body(i2)*zmbar
179  body(i1) = zm
180  body(i2) = 0.d0
181  name1 = name(i1)
182  name2 = name(i2)
183  IF(body0(i1).LT.body0(i2))THEN
184  body0(i1) = body0(i2)
185  epoch(i1) = epoch(i2)
186  tev(i1) = tev(i2)
187  spin(i1) = spin(i2)
188  radius(i1) = radius(i2)
189  radius(i2) = 0.0
190  name2 = name(i2)
191  name(i2) = name(i1)
192  name(i1) = name2
193  ENDIF
194  t0(i1) = time
195  t0(i2) = tadj + dtadj
196  IF (kz(23).EQ.0.OR.rtide.GT.1000.0*rscale) t0(i2) = 1.0d+10
197 * CALL DTCHCK(TIME,STEP(I2),DTK(40))
198  step(i2) = 1.0d+06
199 *
200 * Start new star from current time unless ROCHE case with TEV0 > TIME.
201  tev0(i1) = max(tev0(i1),tev0(i2))
202  IF(iqcoll.NE.3) tev(i1) = max(time,tev0(i1))
203  tev(i2) = 1.0d+10
204  vi = sqrt(xdot(1,i2)**2 + xdot(2,i2)**2 + xdot(3,i2)**2)
205 *
206 * Set T0 = TIME for any other chain members.
207  IF (ipair.LT.0) THEN
208  DO 18 l = 1,nch
209  j = jlist(l)
210  IF (j.NE.i1.AND.j.NE.i2) THEN
211  t0(j) = time
212  END IF
213  18 CONTINUE
214  END IF
215 *
216 * Check that a mass-less primary has type 15 for kick velocity.
217  IF(zmnew*zmbar.LT.0.001.AND.kw1.NE.15)THEN
218  WRITE(6,*)' ERROR COAL: mass1 = 0.0 and kw1 is not equal 15'
219  WRITE(6,*)' I KW ',i1,kw1
220  stop
221  END IF
222 *
223  DO 20 k = 1,3
224  x(k,i1) = cm(k)
225  x0(k,i1) = cm(k)
226  xdot(k,i1) = cm(k+3)
227  x0dot(k,i1) = cm(k+3)
228 * Ensure that ghost will escape next output (far from fast escapers).
229  x0(k,i2) = min(1.0d+04 + (x(k,i2)-rdens(k)),
230  & 1000.d0*rscale*(x(k,i2)-rdens(k))/ri)
231  x(k,i2) = x0(k,i2)
232  x0dot(k,i2) = sqrt(0.004d0*zmass/rscale)*xdot(k,i2)/vi
233  xdot(k,i2) = x0dot(k,i2)
234  f(k,i2) = 0.d0
235  fdot(k,i2) = 0.d0
236  d0(k,i2) = 0.0
237  d1(k,i2) = 0.0
238  d2(k,i2) = 0.d0
239  d3(k,i2) = 0.d0
240  d0r(k,i2) = 0.0
241  d1r(k,i2) = 0.0
242  d2r(k,i2) = 0.d0
243  d3r(k,i2) = 0.d0
244  20 CONTINUE
245 *
246 * Obtain potential energy w.r.t. new c.m. and apply tidal correction.
247  CALL nbpot(1,nnb,pot2)
248  dp = pot2 - pot1
249  ecoll = ecoll + dp
250 *
251  j2 = jpert(l2)
252  jpert(l2) = i2
253 *
254 * Remove the ghost particle #I2 from perturber lists containing #I1.
255  jlist(1) = i2
256  CALL nbrem(i1,1,nnb)
257  jlist(1) = i1
258  jpert(l2) = j2
259 *
260 * Include correction procedure in case of mass loss (cf routine MIX).
261  IF (kz(19).GE.3.AND.dm.GT.0.0) THEN
262 *
263 * Reduce mass of composite body and update total mass (check SN mass).
264  body(i1) = zmnew
265  body(i1) = max(body(i1),0.d0)
266  IF (abs(body(i1)).LT.1.0d-10) tev(i1) = 1.0d+10
267  zmass = zmass - dm
268 *
269 * Adopt ILIST procedure from NBODY4 with NNB available in FCORR.
270  ilist(1) = nnb
271  DO 22 l = 1,nnb
272  ilist(l+1) = jpert(l)
273  22 CONTINUE
274 *
275 * Delay velocity kick until routine MDOT on type 13/14/15 in ROCHE.
276 * [Note: NS,BH should not be created here unless possibly for
277 * AIC of WD or TZ object from COMENV.]
278  kw = kw1
279  IF (kw1.EQ.13.OR.kw1.EQ.14) THEN
280  IF(kstar(i1).GE.13.OR.kstar(i2).GE.13) kw = 0
281 * IF(KSTAR(I1).GE.10.OR.KSTAR(I2).GE.10) KW = 0
282  END IF
283  IF (kw1.GE.10.AND.kw1.LE.12) THEN
284  IF(kstar(i1).GE.10.OR.kstar(i2).GE.10) kw = 0
285  END IF
286 *
287 * Perform total force & energy corrections (new polynomial set later).
288  CALL fcorr(i1,dm,kw)
289 *
290 * Specify commensurate time-step (not needed for BODY(I1) = 0).
291  IF (body(i1).GT.0.0d0) THEN
292  CALL dtchck(time,step(i1),dtk(40))
293  END IF
294 *
295 * Set IPHASE = -3 to preserve ILIST.
296  iphase = -3
297 *
298 * Initialize new polynomials of neighbours & #I1 for DM > 0.1 DMSUN.
299  IF (dm*zmbar.GT.0.1) THEN
300 *
301 * Include body #I1 at the end (counting from location #2).
302  nnb2 = nnb + 2
303  ilist(nnb2) = i1
304 *
305 * Obtain new F & FDOT and time-steps.
306  DO 30 l = 2,nnb2
307  j = ilist(l)
308  IF (l.EQ.nnb2) THEN
309  j = i1
310  ELSE IF (t0(j).LT.time) THEN
311  CALL xvpred(j,-2)
312  CALL dtchck(time,step(j),dtk(40))
313  END IF
314  DO 25 k = 1,3
315  x0dot(k,j) = xdot(k,j)
316  25 CONTINUE
317 * Create ghost for rare case of massless first component.
318  IF (body(j).EQ.0.0d0) THEN
319  DO 26 k = 1,3
320  x0(k,i1) = min(1.0d+04 + (x(k,i1)-rdens(k)),
321  & 1000.d0*rscale*(x(k,i1)-rdens(k))/ri)
322  x(k,i1) = x0(k,i1)
323  x0dot(k,i1) = sqrt(0.004d0*zmass/rscale)*
324  & xdot(k,i1)/vi
325  xdot(k,i1) = x0dot(k,i1)
326  f(k,i1) = 0.d0
327  fdot(k,i1) = 0.d0
328  d2(k,i1) = 0.d0
329  d3(k,i1) = 0.d0
330  26 CONTINUE
331  t0(i1) = 1.0d+06
332  WRITE (6,28) name(i1), kw1
333  28 FORMAT (' MASSLESS PRIMARY! NAM KW ',i8,i4)
334  ELSE
335  CALL fpoly1(j,j,0)
336  CALL fpoly2(j,j,0)
337  END IF
338  30 CONTINUE
339  END IF
340 * TPREV = TIME - STEPX
341  END IF
342 *
343 * See whether closest neighbour forms a KS pair (skip chain).
344  IF (ipair.GT.0.AND.body(i1).GT.0.0d0) THEN
345  IF (jmin.LE.n.AND.rij2.LT.rmin2) THEN
346  DO 35 k = 1,3
347  x0dot(k,jmin) = xdot(k,jmin)
348  35 CONTINUE
349  icomp = min(i1,jmin)
350  jcomp = max(i1,jmin)
351  CALL ksreg
352  WRITE (6,36) name(icomp), name(jcomp), list(1,2*npairs-1),
353  & r(npairs), h(npairs), step(ntot)
354  36 FORMAT (' COAL KS NM NP R H DTCM ',2i6,i4,1p,3e10.2)
355  i2 = jmin
356 * Note that T0(I2) may not have a large value after #I2 is exchanged.
357  t0(i2) = tadj + dtadj
358  ELSE
359 * Initialize force polynomial for new single body.
360  CALL fpoly1(icomp,icomp,0)
361  CALL fpoly2(icomp,icomp,0)
362  END IF
363  END IF
364 *
365 * Update energy loss & collision counters (EB = 0 for CHAIN COAL).
366  ecoll = ecoll + eb
367  e(10) = e(10) + eb
368  npop(9) = npop(9) + 1
369  ncoal = ncoal + 1
370  eb = eb1
371 *
372 * Open unit #12 the first time.
373  IF (first) THEN
374  OPEN (unit=12,status='NEW',form='FORMATTED',file='COAL')
375  first = .false.
376 *
377 * Print cluster scaling parameters at start of the run.
378  IF (ncoal.EQ.1) THEN
379  WRITE (12,40) rbar, bodym*zmbar, body1*zmbar, tscale,
380  & nbin0, nzero
381  40 FORMAT (/,6x,'MODEL: RBAR =',f5.1,' <M> =',f6.2,
382  & ' M1 =',f6.1,' TSCALE =',f6.2,
383  & ' NB =',i5,' N0 =',i6,//)
384  WRITE (12,45)
385  45 FORMAT (' TIME NAME NAME K1 K2 IQ M1 M2',
386  & ' DM R1 R2 r/Rc R ECC P',/)
387  END IF
388  END IF
389 *
390  WRITE (12,50) ttot, name1, name2, kstar(i1), kstar(i2), iqcoll,
391  & zm1, zm2, dm*zmbar, radius(i1)*su, radius(i2)*su,
392  & ri/rc, rij*su, ecc, tk
393  50 FORMAT (1x,f7.1,2i6,3i4,3f5.1,2f7.2,f6.1,f7.2,f9.5,1p,e9.1)
394  CALL flush(12)
395 *
396  WRITE (6,55) which1, iqcoll, name1, name2, kstar(i1), kstar(i2),
397  & kw1, zmnew*zmbar, rcoll, eb, dp, dm*zmbar, vinf
398  55 FORMAT (/,a8,'COAL IQ =',i3,' NAME =',2i6,' K* =',3i3,
399  & ' M =',f6.2,' RCOLL =',1p,e8.1,' EB =',e9.1,
400  & ' DP =',e9.1,' DM =',0p,f6.2,' VINF =',f4.1)
401 *
402  kstar(i1) = kw1
403  kstar(i2) = 15
404 * Specify IPHASE < 0 for new sorting.
405  iphase = -1
406  iqcoll = 0
407 *
408  RETURN
409 *
410  END