Nbody6
 All Files Functions Variables
cmbody.f
Go to the documentation of this file.
1  SUBROUTINE cmbody(ENERGY,NSYS)
2 *
3 *
4 * Formation of c.m. body by collision.
5 * ------------------------------------
6 *
7  include 'common6.h'
8  parameter(nmx=10,nmx4=4*nmx)
9  common/close/ rij(4,4),rcoll4,qperi4,size4(4),ecoll4,ip(4)
10  common/ccoll2/ qk(nmx4),pk(nmx4),rik(nmx,nmx),SIZE(nmx),vstar1,
11  & ecoll1,rcoll,qperi,istar(nmx),icoll,isync,ndiss1
12  common/ebsave/ ebs
13  REAL*8 cm(6),a0(3),a2(3)
14  CHARACTER*8 which1
15  LOGICAL first
16  SAVE first
17  DATA first /.true./
18 *
19 *
20 * Distinguish between chain and triple or quad case (ICH > 0 or = 0).
21  IF (iphase.EQ.9) THEN
22  ich = 1
23  ELSE
24 * Activate collision indicator (otherwise done in CHTERM).
25  ich = 0
26  iphase = 9
27  END IF
28 *
29 * Specify global indices of subsystem (membership: NSYS = 2 - 5).
30  IF (nsys.EQ.2) THEN
31 *
32 * Define discrete time for prediction & new polynomials (T <= TBLOCK).
33  i = n + kspair
34  dt = min(tblock-tprev,dtmin)
35  IF (dt.GT.2.4e-11) THEN
36  time2 = time - tprev
37  CALL stepk(dt,dtn)
38  time = tprev + int((time2 + dt)/dtn)*dtn
39  time = min(tblock,time)
40  ELSE
41  time = min(t0(i) + step(i),tblock)
42  END IF
43  time0 = time
44 *
45 * Check for hierarchical configuration.
46  i1 = 2*kspair - 1
47  i2 = i1 + 1
48  jcl = 0
49  np1 = list(1,i1) + 1
50  semi = -0.5*body(i)/h(kspair)
51  rx = 100.0
52  DO 5 l = 2,np1
53  j = list(l,i1)
54  rij2 = 0.0
55  vij2 = 0.0
56  rdot = 0.0
57  a12 = 0.0
58  a22 = 0.0
59  a1a2 = 0.0
60  DO 2 k = 1,3
61  rij2 = rij2 + (x(k,i) - x(k,j))**2
62  vij2 = vij2 + (xdot(k,i) - xdot(k,j))**2
63  rdot = rdot + (x(k,i) - x(k,j))*(xdot(k,i) -xdot(k,j))
64  k1 = k + 1
65  IF (k1.GT.3) k1 = 1
66  k2 = k1 + 1
67  IF (k2.GT.3) k2 = 1
68  a0(k) = (x(k1,i1)-x(k1,i2))*(xdot(k2,i1)-xdot(k2,i2))
69  & - (x(k2,i1)-x(k2,i2))*(xdot(k1,i1)-xdot(k1,i2))
70  a2(k) = (x(k1,j) - x(k1,i))*(xdot(k2,j) - xdot(k2,i))
71  & - (x(k2,j) - x(k2,i))*(xdot(k1,j) - xdot(k1,i))
72  a12 = a12 + a0(k)**2
73  a22 = a22 + a2(k)**2
74  a1a2 = a1a2 + a0(k)*a2(k)
75  2 CONTINUE
76  rip = sqrt(rij2)
77  a1 = 2.0/rip - vij2/(body(i) + body(j))
78  a1 = 1.0/a1
79  ecc2 = (1.0 - rip/a1)**2 +
80  & rdot**2/(a1*(body(i) + body(j)))
81  IF (1.0/a1.GT.0.5/rmin) THEN
82  ecc1 = sqrt(ecc2)
83  rp1 = a1*(1.0 - ecc1)
84  ecc = 1.0 - r(kspair)/semi
85  ra = semi*(1.0 + ecc)
86  sr = rp1/ra
87  ga = 2.0*body(j)*(ra/rp1)**3/body(i)
88 * Determine inclination (8 bins of 22.5 degrees).
89  fac = a1a2/sqrt(a12*a22)
90  angle = 360.0*acos(fac)/twopi
91  WRITE (6,4) kspair, name(j), h(kspair), ecc, semi,
92  & a1, rp1, ga, ecc1, sr, angle
93  4 FORMAT (' HIERARCHY: KS NMJ H E A0 A1 RP GA E1 SR ',
94  & 'IN',2i6,f7.0,f9.5,1p,4e9.1,0p,f6.2,f6.1,f7.1)
95  ELSE
96  IF (rip.LT.rx) THEN
97  rx = rip
98  jx = j
99  eccx = sqrt(ecc2)
100  semix = a1
101  rdx = rdot/rx
102  END IF
103  END IF
104 *
105 * Select closest single body inside 0.5*RMIN as KS component.
106  IF (rip.LT.0.5*rmin.AND.j.LE.n) THEN
107  IF (jcl.GT.0) THEN
108  IF (rip.GT.rip0) go to 5
109  jcl = j
110  rip0 = rip
111  ELSE
112  jcl = j
113  rip0 = rip
114  END IF
115  END IF
116  5 CONTINUE
117 *
118  IF (rx.LT.20.0*rmin) THEN
119  WRITE (6,200) name(jx), eccx, rx, semix*(1.0 - eccx),
120  & rdx, semi
121  200 FORMAT (' PERTURBER: NM E1 RX PM RD A0 ',
122  & i6,f7.3,1p,4e9.1)
123  END IF
124 *
125 * Search for evidence of recent regularization.
126  nam1 = name(2*kspair-1)
127  nam2 = name(2*kspair)
128  nnb = listd(1)
129  DO 7 k = 2,nnb+1
130  IF (listd(k).EQ.nam1.OR.listd(k).EQ.nam2) THEN
131  WRITE (6,6) nam1, nam2, listd(k), k
132  6 FORMAT (' KS REMNANT: NAM LISTD K ',3i6,i4)
133  END IF
134  7 CONTINUE
135 *
136 * Predict body #JCL to current time in case of no collision.
137  IF (jcl.GT.0) CALL xvpred(jcl,-1)
138 *
139 * Ensure orbit is at pericentre (perturbed hyperbolic case is OK).
140  semi = -0.5*body(n+kspair)/h(kspair)
141  IF (r(kspair).GT.semi.AND.semi.GT.0.0) THEN
142  CALL ksapo(kspair)
143  CALL ksperi(kspair)
144 * Restore quantized time to avoid small STEP (KSTERM needs T0 = TIME).
145  time = time0
146  t0(i1) = time
147  END IF
148 *
149 * Save collision distance and VINF before any common envelope stage.
150  rcoll = r(kspair)
151  vinf = 0.0
152  IF (h(kspair).GT.0.0) vinf = sqrt(2.0*h(kspair))*vstar
153  ecc = 1.0 - r(kspair)/semi
154 *
155 * Include special procedure for common envelope stage with mass loss.
156  IF (kz(19).GE.3) THEN
157  k1 = kstar(i1)
158  k2 = kstar(i2)
159 * PreMS collision scheme: -1, -1 ==> -1; -1, 0 ==> -1; -1, x ==> x.
160  IF (k1.LT.0.OR.k2.LT.0) THEN
161  icase = -1
162  IF (max(k1,k2).EQ.0) icase = -1
163  IF (max(k1,k2).GT.0) icase = max(k1,k2)
164  ELSE
165  icase = ktype(k1,k2)
166  END IF
167  IF(icase.GT.100)THEN
168  iqcoll = 4
169  CALL expel(i1,i2,icase)
170  IF (icase.LT.0) go to 100
171 * Treat collision as before in case of CE without coalescence.
172  icomp = i1
173  END IF
174  END IF
175 *
176 * Update body #JCL to current time for new KS with combined c.m.
177  IF (jcl.GT.0) THEN
178 * CALL XVPRED(JCL,-1)
179  t0(jcl) = time
180  CALL dtchck(time,step(jcl),dtk(40))
181  DO 8 k = 1,3
182  x0dot(k,jcl) = xdot(k,jcl)
183  x0(k,jcl) = x(k,jcl)
184  8 CONTINUE
185  END IF
186 *
187 * Check diagnostics of degenerate binary (skip case of velocity kick).
188  IF(kz(8).GT.3.AND.max(kstar(i1),kstar(i2)).GE.10)THEN
189  IF(kstar(i1).LE.12)THEN
190  CALL degen(kspair,kspair,5)
191  END IF
192  END IF
193 *
194 * Check optional diagnostics for final stage of binary evolution.
195  IF (kz(8).GT.3) THEN
196  CALL binev(kspair)
197  END IF
198 *
199 * Save binding energy (BODY(I2) = 0 is OK).
200  eb = body(i1)*body(i2)*h(kspair)/body(i)
201  which1 = ' BINARY '
202  IF (h(kspair).GT.0.0) THEN
203  which1 = ' HYPERB '
204  nhyp = nhyp + 1
205  END IF
206 *
207 * Terminate KS pair and set relevant indices for collision treatment.
208  kstari = kstar(i)
209  t0(i1) = time
210  semi = -0.5*body(i)/h(kspair)
211  tk = days*semi*sqrt(semi/body(i))
212  CALL dtchck(time,step(i1),dtk(40))
213  CALL ksterm
214  i1 = 2*npairs + 1
215  i2 = i1 + 1
216  i3 = 0
217  icomp = i1
218  dmin2 = min(dmin2,rcoll)
219  ELSE
220 * Note JLIST(1->NCH) contains global indices (JLIST(4)=0 for NCH=3).
221  i1 = jlist(1)
222  i2 = jlist(2)
223  i3 = jlist(3)
224  i4 = jlist(4)
225  IF (nsys.GT.4) i5 = jlist(5)
226  iqcoll = 5
227  kstari = 0
228  vinf = 0.0
229  ecc = 1.0 + 2.0*ebs*dminc/(body(i1)*body(i2))
230  ecc = max(ecc,0.001d0)
231  IF (ebs.GT.0) THEN
232  hi = ebs*(body(i1) + body(i2))/(body(i1)*body(i2))
233  vinf = sqrt(2.0*hi)*vstar
234  END IF
235 *
236 * Set new quantized time (note: suppress if problems in CHTERM).
237  time = tblock
238 *
239 * Include special treatment for common envelope stage inside chain.
240  IF (ich.GT.0.AND.kz(19).GE.3) THEN
241  zm1 = body(i1)*zmbar
242  zm2 = body(i2)*zmbar
243  r1 = radius(i1)*su
244  r2 = radius(i2)*su
245  WRITE (86,9) tphys, name(i1), name(i2), kstar(i1),
246  & kstar(i2), zm1, zm2, r1, r2, dminc*su, ecc
247  9 FORMAT (' COLL: TPH NAM K* M R* QP E ',
248  & f8.1,2i7,2i4,2f6.2,3f7.1,f9.5)
249 *
250  k1 = kstar(i1)
251  k2 = kstar(i2)
252  icase = ktype(k1,k2)
253  IF(icase.GT.100)THEN
254  iqcoll = 6
255  CALL expel2(i1,i2,icase)
256 * Decide between chain restart, coalescence or collision.
257  IF (icase.GT.0) THEN
258 * Adopt negative membership and reverse NSYS to denote chain restart.
259  nch = -nsys
260  nsys = -nsys
261  go to 100
262  END IF
263 * Check for coalescence (one or two bodies of zero mass).
264  IF (body(i1).EQ.0.0d0.OR.body(i2).EQ.0.0d0) THEN
265 * Include rare case of two massless stars (defined by KSTAR = 15).
266  IF (body(i1) + body(i2).EQ.0.0d0) THEN
267  ri = sqrt(x(1,i1)**2 + x(2,i1)**2 +x(3,i1)**2)
268  vi = sqrt(xdot(1,i1)**2 + xdot(2,i1)**2 +
269  & xdot(3,i1)**2)
270  dtmax = dtk(1)
271  CALL dtchck(time,dtmax,dtk(40))
272  WRITE (6,10) kstar(i1), kstar(i2), dtmax
273  10 FORMAT (' MASSLESS CM K* DTX ',2i4,1p,e10.2)
274  go to 32
275  END IF
276 * Reduce the membership (< 0 for SETSYS) and remove ghost from chain.
277  nch = -(nsys - 1)
278  jlist(1) = i1
279  IF (body(i1).EQ.0.0d0) jlist(1) = i2
280  DO 11 l = 2,nsys-1
281  jlist(l) = jlist(l+1)
282  11 CONTINUE
283  go to 100
284  END IF
285 * Treat collision as before in case of CE without coalescence.
286  icomp = i1
287  jlist(1) = i1
288  jlist(2) = i2
289  END IF
290  END IF
291  END IF
292 *
293 * Obtain mass loss and evolution epoch of composite star.
294  dm = 0.0d0
295  IF (kz(19).GE.3) THEN
296  CALL mix(i1,i2,dm)
297  icomp = i1
298 * Note possible switching of I1 and I2 (cf. JLIST).
299  END IF
300 *
301 * Define global c.m. coordinates & velocities from body #I1 & I2.
302  zm = body(i1) + body(i2)
303  DO 12 k = 1,3
304  cm(k) = (body(i1)*x(k,i1) + body(i2)*x(k,i2))/zm
305  cm(k+3) = (body(i1)*xdot(k,i1) + body(i2)*xdot(k,i2))/zm
306  12 CONTINUE
307 *
308 * Set T0 = TIME for correct potential energy correction in FCORR.
309  IF (ich.GT.0) THEN
310  DO 14 l = 1,nch
311  j = jlist(l)
312  t0(j) = time
313  CALL dtchck(time,step(j),dtk(40))
314  14 CONTINUE
315  END IF
316 *
317 * Ensure NAME of the heaviest body is new progenitor.
318  IF (body(i2).GT.body(i1)) THEN
319  nam1 = name(i1)
320  name(i1) = name(i2)
321  name(i2) = nam1
322  END IF
323 *
324 * Evaluate potential energy with respect to colliding bodies.
325  IF (nsys.EQ.2) THEN
326 * Copy perturber list to JPERT.
327  nnb = list(1,i1)
328  DO 15 l = 1,nnb
329  jpert(l) = list(l+1,i1)
330  15 CONTINUE
331  jlist(1) = i1
332  jlist(2) = i2
333 * Replace second old KS component temporarily by arbitrary body.
334  jpert(1) = n
335  IF (i2.EQ.n) jpert(1) = n + 1
336  CALL nbpot(2,nnb,pot1)
337  ELSE
338 * Obtain differential effect on #I1 & #I2 due to other members.
339  DO 16 l = 3,nsys
340  jpert(l-2) = jlist(l)
341  16 CONTINUE
342  np = nsys - 2
343  CALL nbpot(2,np,pot1)
344 * Compensate EGRAV for any chain (or TRIPLE/QUAD) mass loss.
345  IF (dm.GT.0.0d0) THEN
346  eg0 = egrav
347  DO 18 l = 1,nsys
348  j = jlist(l)
349  IF (j.EQ.i1.OR.j.EQ.i2) go to 18
350  rij2 = 0.0
351  DO 17 k = 1,3
352  rij2 = rij2 + (cm(k) - x(k,j))**2
353  17 CONTINUE
354  egrav = egrav - dm*body(j)/sqrt(rij2)
355  18 CONTINUE
356  WRITE (6,19) nsys, egrav, egrav - eg0
357  19 FORMAT (' CHAIN/TRIPLE MASS LOSS NSYS EGRAV DEGR ',
358  & i4,1p,2e10.2)
359  END IF
360  END IF
361 *
362 * Create new body from c.m. and initialize zero mass ghost in #I2.
363  zm1 = body(i1)*zmbar
364  zm2 = body(i2)*zmbar
365  body(i1) = zm
366  body(i2) = 0.0d0
367  spin(i1) = (spin(i1) + spin(i2))*(1.0 - dm/zm)
368  t0(i2) = tadj + dtadj
369  IF (kz(23).EQ.0.OR.rtide.GT.1000.0*rscale) t0(i2) = 1.0d+10
370  dtmax = dtk(1)
371  CALL dtchck(time,dtmax,dtk(40))
372  step(i2) = dtmax
373  ri = sqrt((x(1,i2) - rdens(1))**2 + (x(2,i2) - rdens(2))**2
374  & + (x(3,i2) - rdens(3))**2)
375  vi = sqrt(xdot(1,i2)**2 + xdot(2,i2)**2 + xdot(3,i2)**2)
376  name1 = name(i1)
377  name2 = name(i2)
378 *
379  DO 20 k = 1,3
380  x(k,i1) = cm(k)
381  x0(k,i1) = cm(k)
382  xdot(k,i1) = cm(k+3)
383  x0dot(k,i1) = cm(k+3)
384 * Ensure that ghost will escape next output (far from fast escapers).
385  x0(k,i2) = 1000.0*rscale*(x(k,i2) - rdens(k))/ri
386  x(k,i2) = x0(k,i2)
387  x0dot(k,i2) = sqrt(0.004*zmass/rscale)*xdot(k,i2)/vi
388  xdot(k,i2) = x0dot(k,i2)
389  f(k,i2) = 0.0d0
390  fdot(k,i2) = 0.0d0
391  d2(k,i2) = 0.0d0
392  d3(k,i2) = 0.0d0
393  20 CONTINUE
394 *
395 * Set appropriate parameters for coalescing GR binaries.
396  IF (kstar(i1).GE.13.AND.kz(28).GT.0) THEN
397  tev(i1) = 1.0d+10
398  END IF
399  IF (kstar(i2).GE.13.AND.kz(28).GT.0) THEN
400  tev(i2) = 1.0d+10
401  kstar(i2) = 0
402  END IF
403 *
404 * Refresh index of new dominant body in case of switch in routine MIX.
405  jlist(1) = i1
406 * Obtain potential energy w.r.t. new c.m. and apply tidal correction.
407  IF (nsys.EQ.2) THEN
408  CALL nbpot(1,nnb,pot2)
409  ELSE
410  CALL nbpot(1,np,pot2)
411  END IF
412  dp = pot2 - pot1
413  ecoll = ecoll + dp
414 *
415 * Remove the ghost particle from perturber lists containing #I1.
416  IF (nsys.EQ.2) THEN
417  jpert(1) = i2
418  jlist(1) = i2
419  CALL nbrem(i1,1,nnb)
420 * Remove ghost from list of I1 (use NTOT as dummy here).
421  jpert(1) = i1
422  CALL nbrem(ntot,1,1)
423  ELSE
424 * Determine index and set neighbour membership of original c.m.
425  jlist(2) = i2
426  nnb = 0
427  DO 21 l = 1,nsys
428  j = jlist(l)
429  IF (list(1,j).GT.nnb) THEN
430  nnb = list(1,j)
431  icm = j
432  END IF
433  21 CONTINUE
434 *
435 * Update neighbour lists of current c.m. and remove ghost I2.
436  DO 22 l = 1,nnb
437  jpert(l) = list(l+1,icm)
438  22 CONTINUE
439  CALL nbrest(icm,nsys,nnb)
440  jlist(1) = i2
441  jpert(1) = icm
442  CALL nbrem(ntot,1,1)
443  END IF
444 *
445 * Determine new neighbour list for first CHAIN/TRIPLE member #I1.
446  IF (nsys.GT.2) THEN
447  rsi = rscale*(10.0/float(n - npairs))**0.3333
448  CALL nblist(i1,rsi)
449  END IF
450 *
451 * Include correction procedure in case of mass loss (cf routine MIX).
452  IF (kz(19).GE.3.AND.dm.GT.0.0d0) THEN
453 *
454  nnb = list(1,i1)
455  ilist(1) = nnb
456  DO 24 l = 2,nnb+1
457  ilist(l) = list(l,i1)
458  24 CONTINUE
459 *
460 * Reduce mass of composite body and update total mass (check SN mass).
461  body(i1) = zm - dm
462  body(i1) = max(body(i1),0.0d0)
463  zmass = zmass - dm
464  zm = zm - dm
465 *
466 * Perform total force & energy corrections (new polynomial set later).
467  kw = kstar(i1)
468  CALL fcorr(i1,dm,kw)
469 *
470 * Initialize new polynomials of neighbours & #I for DM > 0.1 DMSUN.
471  IF (dm*zmbar.GT.0.1) THEN
472 *
473 * Include body #I at the end (counting from location #2; not KS case).
474  nnb2 = nnb + 2
475  ilist(nnb2) = i1
476  IF (nsys.EQ.2) nnb2 = nnb2 - 1
477 *
478 * Obtain new F & FDOT and time-steps.
479  DO 30 l = 2,nnb2
480  j = ilist(l)
481  DO 25 k = 1,3
482  x0dot(k,j) = xdot(k,j)
483  25 CONTINUE
484  CALL fpoly1(j,j,0)
485  CALL fpoly2(j,j,0)
486  30 CONTINUE
487  END IF
488 * Note danger of multiple collisions on same block-step (early times).
489 * TPREV = TIME - STEPX
490  END IF
491 *
492 * Check creation of ghost(s) after collision of two white dwarfs.
493  32 IF (kstar(i1).EQ.15) THEN
494  t0(i1) = tadj + dtadj
495  step(i1) = dtmax
496  DO 35 k = 1,3
497  x0(k,i1) = 1000.0*rscale*x(k,i1)/ri
498  x(k,i1) = x0(k,i1)
499  x0dot(k,i1) = sqrt(0.004*zmass/rscale)*xdot(k,i2)/vi
500  xdot(k,i1) = x0dot(k,i1)
501  f(k,i1) = 0.0d0
502  fdot(k,i1) = 0.0d0
503  d2(k,i1) = 0.0d0
504  d3(k,i1) = 0.0d0
505  35 CONTINUE
506 * Include case of two ghost stars.
507  IF (body(i2).EQ.0.0d0.AND.i1.NE.i2) THEN
508  i1 = i2
509  go to 32
510  END IF
511 * Perform initialization of single particle in three-body chain.
512  IF (nsys.EQ.3) THEN
513  icomp = i3
514 * Make quick exit from routine CHAIN on zero membership.
515  nsys = 0
516  nch = 0
517  ech = 0.0
518  go to 40
519  END IF
520 * Include treatment for larger chain here when needed (rare case).
521  IF (nsys.GT.3) THEN
522  WRITE (6,38) name(i1)
523  38 FORMAT (' DANGER! ZERO MASS IN CHAIN NAME',i7)
524  stop
525  END IF
526  END IF
527 *
528 * Decide appropriate path for each case.
529  IF (nsys.EQ.2) go to 40
530  IF (nsys.EQ.3) go to 45
531 *
532 * Switch KS components if body #I3 & I4 is closer than #I1 & I3.
533  IF (jlist(7).LT.0) THEN
534  i4 = i1
535  i1s = i1
536  i1 = jlist(4)
537  jlist(4) = i1s
538  END IF
539  icomp = i4
540 *
541 * Check KS case for new regularization with close hierarchical body.
542  40 IF (nsys.EQ.2) THEN
543  IF (jcl.GT.0) THEN
544  icomp = i1
545  jcomp = jcl
546  CALL ksreg
547  go to 80
548  END IF
549  END IF
550 *
551 * Initialize force polynomial for new single, third or fourth body.
552  CALL fpoly1(icomp,icomp,0)
553  CALL fpoly2(icomp,icomp,0)
554  IF (nsys.EQ.5) THEN
555  CALL fpoly1(i5,i5,0)
556  CALL fpoly2(i5,i5,0)
557  WRITE (6,42) name(i5), step(i5)
558  42 FORMAT (' 5-CHAIN! NM DT ',i7,1p,e10.2)
559  END IF
560  IF (nsys.EQ.0) go to 95
561  IF (nsys.LE.2) go to 80
562 *
563 * Add kinetic energy from last body and check DMIN in TRIPLE or QUAD.
564  45 IF (ich.GT.0) THEN
565  which1 = ' CHAIN '
566 * Specify new membership (< 0 for SETSYS) and remove ghost from chain.
567  nch = -(nsys - 1)
568  jlist(1) = i1
569  DO 50 l = 2,nsys-1
570  jlist(l) = jlist(l+1)
571  50 CONTINUE
572 * Copy well defined binding energy and skip explicit evaluation.
573  eb = ebs
574  chcoll = chcoll + eb
575  go to 80
576  ELSE IF (nsys.EQ.3) THEN
577  i = i3
578  which1 = ' TRIPLE '
579  dmin3 = min(dmin3,rcoll4)
580  rcoll = rcoll4
581  eb = ebs
582  ELSE
583  i = i4
584  which1 = ' QUAD '
585  dmin4 = min(dmin4,rcoll4)
586  rcoll = rcoll4
587  eb = ebs
588  END IF
589 *
590 * Set global components for new KS regularization (ICOMP < JCOMP).
591  icomp = min(i1,i3)
592  jcomp = max(i1,i3)
593 *
594 * Initialize new KS pair.
595  CALL ksreg
596 *
597 * Update energy loss & collision counters.
598  80 ecoll = ecoll + eb
599  e(10) = e(10) + eb + dp
600  egrav = egrav + eb
601 *
602 * Open the second coalescence unit #26 first time.
603  IF (first.AND.(iqcoll.EQ.3.OR.kstari.GE.10)) THEN
604  OPEN (unit=26,status='NEW',form='FORMATTED',file='COAL2')
605  first = .false.
606 *
607 * Print cluster scaling parameters at start of the run.
608  WRITE (26,82) rbar, bodym*zmbar, body1*zmbar, tscale,
609  & nbin0, nzero
610  82 FORMAT (/,4x,'MODEL: RBAR =',f5.1,' <M> =',f6.2,
611  & ' M1 =',f6.1,' TSCALE =',f6.2,
612  & ' NB =',i4,' N0 =',i6,//)
613 *
614  WRITE (26,84)
615  84 FORMAT (' TIME NAME NAME K1 K2 IQ M1 M2',
616  & ' DM R1 R2 r/Rc R ECC P',/)
617  END IF
618 *
619 * Distinguish case of contact binary (i.e. coalescence).
620  IF (iqcoll.EQ.3.OR.kstari.GE.10) THEN
621  npop(8) = npop(8) + 1
622  ncoal = ncoal + 1
623  WRITE (6,85) iqcoll, name1, name2, zm*smu, rcoll, eb, dp, ecc
624  85 FORMAT (/,' BINARY COAL IQCOLL =',i3,' NAME =',2i6,
625  & ' M =',f6.2,' RCOLL =',1p,e8.1,' EB =',e9.1,
626  & ' DP =',e9.1,' E =',0p,f8.4)
627 *
628  r1 = radius(i1)*su
629  r2 = radius(i2)*su
630  WRITE (26,86) ttot, name1, name2, kstar(i1), kstar(i2),
631  & iqcoll, zm1, zm2, dm*zmbar, r1, r2, ri/rc,
632  & rcoll*su, ecc, tk
633  86 FORMAT (1x,f7.1,2i6,3i4,3f5.1,2f7.2,f6.2,f7.2,f9.5,1p,e9.1)
634  CALL flush(26)
635  go to 95
636  END IF
637 *
638  npop(8) = npop(8) + 1
639  ncoll = ncoll + 1
640 *
641  WRITE (6,90) which1, nsys, name1, name2, zm*smu, rcoll, eb,
642  & vinf, ecc, dp
643  90 FORMAT (/,a8,'COLLISION NSYS =',i3,' NAME =',2i6,
644  & ' M =',f6.2,' RCOLL =',1p,e8.1,' EB =',e9.1,
645  & ' VINF =',0p,f5.1,' ECC =',f9.5,' DP =',1p,e9.1)
646 *
647 * Specify IPHASE < 0 for new sorting.
648  95 iphase = -1
649 *
650 * Reduce NSUB for chain (temporary increase by CHINIT before CHTERM).
651  100 IF (ich.GT.0) THEN
652  nsub = max(nsub - 1,0)
653  END IF
654 * Skip NSUB reduction for continuation of CHAIN (bug fix 26/8/06).
655  ttot = time + toff
656 *
657  RETURN
658 *
659  END