Nbody6
 All Files Functions Variables
impact.f
Go to the documentation of this file.
1  SUBROUTINE impact(I)
2 *
3 *
4 * Multiple collision or merger search.
5 * ------------------------------------
6 *
7  include 'common6.h'
8  common/binary/ cm(4,mmax),xrel(3,mmax),vrel(3,mmax),
9  & hm(mmax),um(4,mmax),umdot(4,mmax),tmdis(mmax),
10  & namem(mmax),nameg(mmax),kstarm(mmax),iflagm(mmax)
11  common/clump/ bodys(ncmax,5),t0s(5),ts(5),steps(5),rmaxs(5),
12  & names(ncmax,5),isys(5)
13  CHARACTER*8 which1
14  REAL*8 xx(3,3),vv(3,3)
15  INTEGER listq(100)
16  SAVE listq,qcheck
17  DATA izare,listq(1),qcheck /0,0,0.0d0/
18 *
19 *
20 * Set index of KS pair & both components of c.m. body #I.
21  ipair = i - n
22  i1 = 2*ipair - 1
23  i2 = i1 + 1
24  nttry = nttry + 1
25  jcomp = ifirst
26  ks2 = 0
27  rmax2 = 1.0
28  ttot = time + toff
29  ri2 = (x(1,i) - rdens(1))**2 + (x(2,i) - rdens(2))**2 +
30  & (x(3,i) - rdens(3))**2
31 *
32 * Search c.m. neighbours if binary has at most two perturbers.
33  j1 = i1
34  IF (list(1,j1).LE.2) j1 = i
35  nnb2 = list(1,j1) + 1
36 *
37 * Find the dominant body (JCOMP) and nearest perturber (JMAX).
38  rcrit2 = 1.0d+04*rmin2
39  iter = 0
40  5 pert1 = 0.0
41  pert2 = 0.0
42  np = 0
43  DO 10 l = 2,nnb2
44  j = list(l,j1)
45  rij2 = (x(1,i) - x(1,j))**2 + (x(2,i) - x(2,j))**2 +
46  & (x(3,i) - x(3,j))**2
47  IF (j1.EQ.i.AND.rij2.GT.rcrit2) go to 10
48  np = np + 1
49  jlist(np) = j
50  pert = body(j)/(rij2*sqrt(rij2))
51  IF (pert.GT.pert2) THEN
52  IF (pert.GT.pert1) THEN
53  rjmin2 = rij2
54  jmax = jcomp
55  jcomp = j
56  pert2 = pert1
57  pert1 = pert
58  ELSE
59  jmax = j
60  pert2 = pert
61  rmax2 = rij2
62  END IF
63  END IF
64  10 CONTINUE
65 *
66 * Perform an iteration to ensure at least two perturbers (< 10 times).
67  IF (np.LT.2) THEN
68  rcrit2 = 4.0*rcrit2
69  iter = iter + 1
70  IF (iter.LT.10) go to 5
71  END IF
72 *
73  rdot = (x(1,i) - x(1,jcomp))*(xdot(1,i) - xdot(1,jcomp)) +
74  & (x(2,i) - x(2,jcomp))*(xdot(2,i) - xdot(2,jcomp)) +
75  & (x(3,i) - x(3,jcomp))*(xdot(3,i) - xdot(3,jcomp))
76 *
77 * Specify larger perturbation for optional chain regularization.
78  IF ((kz(30).GT.0.OR.kz(30).EQ.-1).AND.nch.EQ.0) THEN
79  gstar = 100.0*gmin
80  kchain = 1
81  ELSE
82  gstar = gmin
83  kchain = 0
84 * Specify indicator -1 for allowing TRIPLE & QUAD but not CHAIN.
85  IF (kz(30).EQ.-2) kchain = -1
86  END IF
87 *
88 * Only accept inward motion or small secondary perturbation.
89  pert3 = 2.0*r(ipair)**3*pert2/body(i)
90  IF (rdot.GT.0.0.OR.pert3.GT.100.0*gstar) go to 100
91 *
92 * Include impact parameter test to distinguish different cases.
93  a2 = (xdot(1,i) - xdot(1,jcomp))**2 +
94  & (xdot(2,i) - xdot(2,jcomp))**2 +
95  & (xdot(3,i) - xdot(3,jcomp))**2
96  rij = sqrt(rjmin2)
97  a3 = 2.0/rij - a2/(body(i) + body(jcomp))
98  semi1 = 1.0/a3
99  a4 = rdot**2/(semi1*(body(i) + body(jcomp)))
100  ecc1 = sqrt((1.0d0 - rij/semi1)**2 + a4)
101  pmin = semi1*(1.0d0 - ecc1)
102 *
103 * Set semi-major axis, eccentricity & apocentre of inner binary.
104  semi = -0.5d0*body(i)/h(ipair)
105  a0 = semi
106  ecc2 = (1.0d0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(body(i)*semi)
107  ecc = sqrt(ecc2)
108  apo = abs(semi)*(1.0 + ecc)
109 *
110 * Quit on hyperbolic orbit with large impact parameter.
111  IF (ecc1.GT.1.0.AND.pmin.GT.50.0*semi) go to 100
112 *
113 * Include rectification for non-circular binaries with KSTAR = 10 & 12.
114  IF (kz(27).EQ.2.AND.kstar(i).GE.10.AND.kstar(i).LE.18) THEN
115  IF (ecc.GT.0.1.AND.mod(kstar(i),2).EQ.0) THEN
116  rm = max(radius(i1),radius(i2))
117  icirc = -1
118  CALL induce(ipair,emax,emin,icirc,tc,angle,tg,edav)
119  WRITE (6,15) name(i1), name(i2), kstar(i1), kstar(i2),
120  & kstar(i), list(1,i1), ecc, semi, rm, pmin,
121  & gamma(ipair), tc, 360.0*angle/twopi, emax
122  15 FORMAT (' NON-CIRCULAR NM K* NP E A R* PM G TC IN EX ',
123  & 2i7,4i4,f7.3,1p,5e9.1,0p,2f7.2)
124 * Circularize the orbit instantaneously for short TC.
125  IF (tc.LT.-100.0) THEN
126 * Set temporary unperturbed orbit (avoids new KSPOLY in DEFORM).
127  np = list(1,i1)
128  list(1,i1) = 0
129  dt1 = step(i1)
130  time0 = time
131  CALL ksrect(ipair)
132  qp = semi*(1.0 - ecc)
133  err = abs(qp - r(ipair))/qp
134 * Deform orbit to circular eccentricity after determining apocentre.
135  IF (r(ipair).GT.semi.OR.err.GT.1.0d-04) THEN
136 * Reduce eccentric anomaly by pi for inward motion.
137  IF (tdot2(ipair).LT.0.0d0) THEN
138  CALL ksapo(ipair)
139  END IF
140  END IF
141 * Predict to pericentre and transform by pi to exact apocentre.
142  CALL ksperi(ipair)
143  CALL ksapo(ipair)
144  eccm = 0.002
145  CALL deform(ipair,ecc,eccm)
146  list(1,i1) = np
147 * Resolv X & XDOT and initialize KS polynomial at apocentre time.
148  IF (np.EQ.0) dt1 = 0.0
149  time = time0 - dt1
150  CALL resolv(ipair,1)
151  CALL kspoly(ipair,1)
152  ELSE
153  kstar(i) = 0
154  END IF
155  END IF
156  END IF
157 *
158 * Form binding energy of inner & outer binary.
159  eb = body(i1)*body(i2)*h(ipair)/body(i)
160  IF(abs(eb).LT.1.0d-10) eb = -1.0d-10
161  eb1 = -0.5*body(jcomp)*body(i)/semi1
162 *
163 * Obtain the perturbing force on body #I & JCOMP (note the skip).
164  CALL fpert(i,jcomp,np,pert)
165 *
166 * Choose maximum of dominant scalar & total vectorial perturbation.
167  pert = pert*rjmin2/(body(i) + body(jcomp))
168  pert4 = 2.0*rjmin2*rij*pert2/(body(i) + body(jcomp))
169  pertm = max(pert4,pert)
170 *
171 * Use combined semi-major axis for binary-binary collision.
172  IF (jcomp.GT.n) THEN
173  jpair = jcomp - n
174  semi2 = -0.5d0*body(jcomp)/h(jpair)
175  j1 = 2*jpair - 1
176  eb2 = -0.5*body(j1)*body(j1+1)/semi2
177 * Define SEMI0 as smallest binary in case IPAIR denotes widest pair.
178  semi0 = min(abs(semi),abs(semi2))
179  semix = max(semi,semi2)
180  apo = apo + max(abs(semi2),r(jpair))
181  semi = semi + semi2
182 * Do not allow negative or soft cross section.
183  IF (1.0/semi.LT.0.5/rmin) go to 100
184 * Consider merger for PMIN > SEMI and large semi-major axis ratio.
185  IF (pmin.GT.semi.AND.semi2.GT.20.0*semi0) go to 30
186  END IF
187 *
188 * Check separation in case of chain regularization.
189  IF (kchain.GT.0) THEN
190 * Form effective gravitational radius (combine triple & quad).
191  ebt = eb + eb1
192  zmm = body(i1)*body(i2) + body(i)*body(jcomp)
193 * Set length of chain for decision-making (also used at termination).
194  rsum = r(ipair) + rij
195  ri = r(ipair)
196  IF (jcomp.GT.n) THEN
197  ebt = ebt + eb2
198  zmm = zmm + body(j1)*body(j1+1)
199  rsum = rsum + r(jpair)
200  ri = max(r(jpair),ri)
201  END IF
202  rgrav = zmm/abs(ebt)
203 * Employ initial size as upper limit in case of weakly bound system.
204  rgrav = min(rgrav,rmin)
205 * Save initial energy in binaries for routine SETSYS.
206  ebch0 = ebt - eb1
207 * Use RIJ instead of RSUM in 3*RGRAV test (increases initial RIJ).
208  IF (rij.GT.max(3.0*rgrav,rmin).OR.rsum.GT.2.0*rmin) go to 30
209  gi = 2.0*body(jcomp)*(ri/rij)**3/body(i)
210 * Enforce KS orbit using MERGE for high eccentricity if PMIN > 10*RI.
211  IF (ecc1.GT.0.99.AND.pmin.GT.10.0*ri.AND.
212  & pertm.LT.gmax) go to 40
213  IF (kz(27).GT.0.AND.jcomp.GT.n) THEN
214  IF (semi0.LT.semi2) j1 = i1
215  rt = 4.0*max(radius(j1),radius(j1+1))
216 * Do not allow large distance ratio for nearly synchronous binary.
217  IF (semi0.GT.rt.AND.ri.GT.25.0*semi0) go to 30
218  IF (min(semi0,semi2).LT.0.05*rij) THEN
219  IF (max(semi0,semi2).LT.0.1*rij) go to 30
220  END IF
221  END IF
222  END IF
223 *
224 * Include special case of strong interraction and large ECC1.
225  IF (ecc1.GT.0.9.AND.gamma(ipair).GT.0.01) THEN
226  IF (apo.LT.0.01*rmin.AND.pmin.LT.2.5*apo) go to 16
227  END IF
228 *
229 * Adopt triple, quad or chain regularization for strong interactions.
230 * IF ((APO.GT.0.01*RMIN.OR.JCOMP.GT.N).AND.PMIN.GT.APO) GO TO 30
231  IF ((apo.GT.0.01*rmin.OR.jcomp.GT.n).AND.pmin.GT.1.5*apo) go to 30
232 * IF (APO.LE.0.01*RMIN.AND.PMIN.GT.2.0*APO) GO TO 30
233  IF ((rij.GT.rmin.AND.semi1.GT.0.0).OR.rij.GT.2.0*rmin) go to 100
234  IF (pertm.GT.100.0*gstar) go to 30
235  16 IF (jcomp.GT.n.AND.pmin.GT.0.1*rmin) THEN
236  IF (pmin.GT.a0 + semi2) go to 30
237  END IF
238  IF (jcomp.GT.n.AND.pmin.GT.4.0*semix.AND.
239  & (ecc1.GT.0.9.AND.ecc1.LT.1.0)) go to 30
240 *
241 * Check almost stable triples (factor 1.2 is experimental).
242  IF (jcomp.LE.n.AND.pmin.GT.2.5*semi) THEN
243  CALL histab(ipair,jcomp,pmin,rstab)
244  ra = semi1*(1.0 + ecc1)
245  IF (semi1.LT.0.0) ra = rij
246  gi = pert*(ra/rij)**3
247 * Use estimated apocentre perturbation for decision-making.
248  IF (pmin.GT.1.2*rstab) THEN
249  IF (gi.LT.0.05) go to 30
250 * Choose chain for critical case of highly eccentric outer orbit.
251  IF (ecc1.LT.0.95) go to 100
252  ELSE IF (pmin.GT.0.7*rstab) THEN
253 * Treat marginally stable triple according to external perturbation.
254  IF (gi.LT.0.05) go to 30
255  IF (gi.LT.1.0.OR.ecc1.LT.0.9) go to 100
256  END IF
257  IF (pmin.GT.0.6*rstab.AND.pmin.LT.0.9*rstab) go to 100
258 * Delay for large distance ratio outside 0.5*RMIN.
259  IF (rij.GT.max(10.0*apo,0.5*rmin)) go to 100
260  IF (pmin.GT.2.5*apo) go to 40
261  END IF
262 *
263  IF (jcomp.GT.n) THEN
264  IF (rij.GT.10.0*apo) go to 100
265  END IF
266 * Skip chain if merged binary or chain c.m. (defined by NAME <= 0).
267  IF (name(i).LE.0.OR.name(jcomp).LE.0) go to 100
268 *
269 * Compare with existing subsystem of same type (if any).
270  IF (nsub.GT.0.AND.kchain.LE.0) THEN
271  perim = r(ipair) + rij
272  IF (jcomp.GT.n) perim = perim + r(jpair)
273  igo = 0
274  CALL permit(perim,igo)
275  IF (igo.GT.0) go to 100
276  END IF
277 *
278 * Skip all multiple regs on zero option (mergers done by #15 > 0).
279  IF (kz(30).EQ.0) go to 100
280 * Allow CHAIN only (#30 = -1) or TRIPLE & QUAD only (#30 < -1).
281  IF (kz(30).EQ.-2.AND.kchain.EQ.0) go to 100
282 *
283  which1 = ' TRIPLE '
284  IF (jcomp.GT.n) which1 = ' QUAD '
285  IF (kchain.GT.0) which1 = ' CHAIN '
286 *
287  IF (h(ipair).GT.0.0) THEN
288  WRITE (6,18) i, jcomp, ecc, ecc1, semi1, rij, gamma(ipair)
289  18 FORMAT (' HYP CHAIN I J E E1 A1 RIJ G ',
290  & 2i6,2f7.3,1p,3e9.1)
291  END IF
292 *
293  IF (kz(15).GT.1.OR.kz(30).GT.1) THEN
294  WRITE (6,20) which1, ipair, ttot, h(ipair), r(ipair),
295  & body(i), body(jcomp), pert4, rij, pmin,
296  & eb1/eb, list(1,i1)
297  20 FORMAT (/,' NEW',a8,i4,' T =',f8.2,' H =',f6.0,
298  & ' R =',1p,e8.1,' M =',2e8.1,' G4 =',1p,e8.1,
299  & ' R1 =',e8.1,' P =',e8.1,' E1 =',0p,f6.3,
300  & ' NP =',i2)
301  CALL flush(3)
302  END IF
303 *
304 * Include any close single or c.m. perturber (cf. routine SETSYS).
305  IF (jmax.NE.jcomp.AND.sqrt(rmax2).LT.min(2.0d0*rsum,rmin).AND.
306  & name(jmax).GT.0) THEN
307  IF (jcomp.GT.n.AND.jmax.GT.n) THEN
308  jcmax = 0
309  ELSE
310  WRITE (6,21) name(jcomp), name(jmax), rsum, sqrt(rmax2)
311  21 FORMAT (' B+2 CHAIN NAM RSUM RMX ',2i7,1p,2e10.2)
312  CALL xvpred(jmax,-1)
313  jcmax = jmax
314  END IF
315  ELSE
316  jcmax = 0
317  END IF
318 *
319 * Save global index of intruder for TRIPLE or CHAIN.
320  jclose = jcomp
321 *
322 * Check B-B interaction for switch of IPAIR & JPAIR or inert binary.
323  IF (kchain.GT.0.AND.jcomp.GT.n) THEN
324  k1 = 2*jpair - 1
325  WRITE (6,22) name(i1), name(i2), name(k1), name(k1+1),
326  & kstar(i), kstar(jcomp), ecc, ecc1, a0, semi2,
327  & rij, semi1, pmin
328  22 FORMAT (' CHAIN B-B NAME K* E0 E1 A0 A2 RIJ A1 PM ',
329  & 4i7,2i4,2f7.3,1p,5e10.2)
330  rt = 4.0*max(radius(i1),radius(i2))
331  IF (semi0.LT.4.0*rt.AND.list(1,j1).EQ.0.OR.
332  & min(semi0,semi2).LT.0.01*rij) THEN
333 * Ensure that widest binary comes first (more similar to triple).
334  IF (semi0.LT.semi2) THEN
335  kpair = jpair
336  jpair = ipair
337  ipair = kpair
338  jclose = n + jpair
339  END IF
340 * Check reduction of c.m. index (JPAIR becomes JPAIR - 1 if > IPAIR).
341  IF (jpair.GT.ipair) jclose = jclose - 1
342 * Include extra condition for inert binary approximation (9/3/12).
343  IF (kz(26).LT.2.AND.rij.GT.250.0*semi0) THEN
344 * Replace unperturbed near-synchronous binary by inert body in CHAIN.
345  jcomp = 0
346  WRITE (6,25) semi0, rij, r(jpair), gamma(jpair)
347  25 FORMAT (' INERT BINARY A RIJ R G ',1p,4e10.2)
348 * Note compact binary may end up as KS if another component escapes.
349  ELSE
350  jclose = 0
351  END IF
352  ELSE
353  jclose = 0
354  END IF
355  END IF
356 *
357 * Set phase indicator for calling TRIPLE or QUAD from MAIN.
358  iphase = 4
359  kspair = ipair
360 *
361 * Include the case of two interacting KS pairs.
362  IF (jcomp.GT.n) THEN
363  iphase = 5
364 * Switch pair indices and rename JCOMP if JPAIR has smaller size.
365  IF (step(j1).LT.step(i1).AND.list(1,i1).GT.0) THEN
366  kspair = jpair
367  jcomp = i
368  ks2 = ipair
369  ELSE
370  ks2 = jpair
371  END IF
372  IF (kz(27).LE.0.AND.jpair.GT.ipair) THEN
373  IF (jclose.GT.0) jclose = jclose - 1
374  END IF
375 * Terminate smallest pair first and reduce second index if higher.
376 * CALL KSTERM
377  IF (ks2.GT.kspair) ks2 = ks2 - 1
378  END IF
379 *
380 * See whether chain regularization indicator should be switched on.
381  IF (kchain.GT.0) THEN
382  iphase = 8
383  END IF
384 *
385 * Save KS indices and delay initialization until end of block step.
386  CALL delay(kchain,ks2)
387 *
388 * Prepare procedure for chain between hierarchy and single body (9/99).
389  IF (name(i).LT.0.AND.name(i).GE.-nzero.AND.jcomp.LE.n) THEN
390 * Indentify merged ghost particle JG.
391  CALL findj(i1,jg,im)
392  WRITE (6,28) name(i), name(jcomp), name(jg),ecc1, pmin, rij
393  28 FORMAT (' HI CHAIN NAM E1 PM RIJ ',i7,2i6,f7.3,1p,2e10.2)
394  jj = jcomp
395 * Terminate the merger in the usual way.
396  kspair = ipair
397  iphase = 7
398  CALL reset
399  zmu = body(2*npairs-1)*body(2*npairs)/body(ntot)
400  ebch0 = ebch0 + zmu*h(npairs)
401 * Specify chain indicator and define the two single particles.
402  iphase = 8
403  jcmax = jg
404  jclose = jj
405  kspair = npairs
406 * Set relevant variables in DELAY before terminating inner binary.
407  CALL delay(kchain,ks2)
408  CALL delay(iphase,-1)
409 * Initialize new chain of the 4 members JMAX, JCLOSE & KS components.
410  isub = 0
411  CALL chain(isub)
412 * Note that IPHASE = -1 now and INTGRT goes back to the beginning.
413  ELSE IF (name(i).LT.-nzero.OR.name(jcomp).LT.0.OR.
414  & (name(i).LT.0.AND.jcomp.GT.n)) THEN
415 * Continue until KS termination on MERGE2 or merger with JCOMP > N.
416  iphase = 0
417  END IF
418 *
419  go to 100
420 *
421 * Begin check for merger of stable hierarchical configuration.
422  30 ra = semi1*(1.0 + ecc1)
423  IF (semi1.LT.0.0) ra = rij
424 *
425 * Identify formation of wide quadruple before merger is accepted.
426  IF (jcomp.GT.n.AND.ecc1.LT.1.0.AND.semi1.LT.0.1*rscale) THEN
427  nnb = listq(1) - 1
428  k = 0
429 * See whether current list contains first inner/outer component.
430  nam1 = name(2*jpair-1)
431  DO 32 l = 2,nnb+2
432  IF (nam1.EQ.listq(l)) k = k + 1
433  32 CONTINUE
434 * Generate diagnostics of first five outer orbits every half period.
435  IF (k.LE.5.AND.time.GT.qcheck.AND.kz(37).GT.0) THEN
436  zmb = body(i) + body(jcomp)
437  ri = sqrt(ri2)
438  tk = semi1*sqrt(semi1/zmb)
439  qcheck = time + min(0.5*twopi*tk,0.1*tcr)
440  tk = days*tk
441 * Check the stability criterion.
442  pcr = stability(body(i1),body(i2),body(jcomp),ecc,ecc1,
443  & 0.0d0)*semix
444  WRITE (89,33) ttot, name(2*ipair-1), nam1, k, ri,
445  & ecc1, eb, eb2, eb1, tk, pmin, pcr
446  33 FORMAT (' QUAD T NAM LQ RI E1 EB EB2 EB1 P1 PM PC ',
447  & f8.1,2i6,i4,f6.2,f8.4,1p,3e12.3,3e9.1)
448  CALL flush(89)
449 * Remove two oldest members if list is too big.
450  IF (nnb.GT.96) THEN
451  DO 34 k = 2,nnb
452  listq(k) = listq(k+2)
453  34 CONTINUE
454  nnb = nnb - 2
455  END IF
456 * Add current names (inner & outer) at end and update membership.
457  listq(nnb+3) = name(2*ipair-1)
458  listq(nnb+4) = name(2*jpair-1)
459  listq(1) = nnb + 3
460  END IF
461  END IF
462 *
463 * Allow temporary merger of inner part of extremely eccentric orbit.
464  rfac = 10.0*rmin
465  IF (ecc1.GT.0.99.AND.ra.GT.rfac) THEN
466  IF (rij.LT.0.1*semi1) rfac = ra
467  END IF
468 *
469 * Increase apocentre tolerance to local scale factor for EB1 < EBS.
470  ebs = 0.25*ebh/sqrt(1.0 + sqrt(ri2)/rscale)
471  IF (eb1.LT.ebs) THEN
472  h2 = (rc**2 + ri2)/float(nc+10)**0.66667
473  rh = 6.0*sqrt(h2/cmsep2)
474  rfac = max(rfac,rh)
475 * Extend maximum apocentre for massive systems (less perturbers).
476  IF (body(i) + body(jcomp).GT.10.0*bodym) rfac = 2.0*rfac
477  END IF
478 
479 * Skip merger for hyperbolic & soft binding energy or large apocentre.
480 * ZF = 1.0
481 * IF (BODY(I)*SMU.LT.0.4) ZF = 1.5
482  IF (min(body(i1),body(i2)).LT.0.05*bodym) THEN
483 * Use orbital velocity condition instead of binding energy for planets.
484  IF (semi1.GT.2.0*rmin) go to 100
485  ELSE IF (eb.GT.ebh.OR.eb1.GT.ebs.OR.ra.GT.rfac) THEN
486  go to 100
487  END IF
488 *
489 * Check tidal capture option (synchronous or evolving binary orbit).
490  IF (kz(27).GT.0) THEN
491 * Skip merger if outer component would suffer tidal dissipation.
492 *** IF (SEMI1*(1.0 - ECC1).LT.4.0*RADIUS(JCOMP)) GO TO 100
493 * Do not allow merger if Roche overflow or mass loss during next orbit.
494  tk = twopi*semi1*sqrt(semi1/(body(i) + body(jcomp)))
495  tm = min(tev(i1),tev(i2),tev(jcomp),tev(i))
496  IF (kz(34).GT.0.AND.tm - time.LT.stepx) THEN
497  rt = 5.0*max(radius(i1),radius(i2))
498  IF (a0.LT.rt.OR.kstar(i).GT.0) THEN
499  CALL trflow(ipair,dtr)
500  IF ((mod(kstar(i),2).EQ.1.AND.dtr.LT.stepx).OR.
501  & dtr.LT.tk) go to 100
502  END IF
503  IF (jcomp.GT.n.AND.kstar(jcomp).GT.0) THEN
504  CALL trflow(jpair,dtr)
505  IF(mod(kstar(jcomp),2).EQ.1.OR.dtr.LT.stepx) goto 100
506  END IF
507  END IF
508 *
509 * Ensure SLEEP for circularizing binary with TCIRC > 1.0E+06.
510  qperi = a0*(1.0 - ecc)
511  rm = max(radius(i1),radius(i2))
512  IF (kstar(i).EQ.-2.AND.qperi.GT.10.0*rm) THEN
513  icirc = -1
514  CALL induce(ipair,emax,emin,icirc,tc,angle,tg,edav)
515  IF (tc.GT.1.0d+06) THEN
516  WRITE (6,35) name(i1), kstar(i1), kstar(i2), ecc,
517  & emax, qperi/rm, edav, a0, pmin, tc
518  35 FORMAT (' IMPACT SLEEP NM K* E EX QP/R ED A PM TC',
519  & i7,2i4,2f8.4,f6.1,1p,4e10.2)
520  kstar(i) = 0
521  nslp = nslp + 1
522  ii = -i
523  CALL spiral(ii)
524  END IF
525  END IF
526 *
527 * Delay merger for recently updated standard binary and short TCIRC.
528  dt = min(tev(i1),tev(i2)) - time
529  IF (kstar(i).EQ.0.AND.name(i).GT.0.AND.dt.LT.tk) THEN
530  icirc = -1
531  CALL tcirc(qperi,ecc,i1,i2,icirc,tc)
532 * WRITE (6,36) NAME(I1), ECC, TTOT, RADIUS(I1)*SU, QPERI, TC
533 * 36 FORMAT (' TCIRC NAM E T R* QP TC ',
534 * & I6,F7.3,F8.3,F7.1,1P,2E10.2)
535 * Beware possible termination by routine HMDOT using QPERI < 3*RADIUS.
536  IF (tc.LT.2000.0.AND.ecc.GT.0.002) go to 100
537  END IF
538  IF (kz(19).GE.3) THEN
539  IF (min(tev(i1),tev(i2)).LT.time + tk) go to 100
540  END IF
541 * Skip chaotic binary (KSTAR = -2 is treated below).
542  IF (kstar(i).EQ.-1.OR.kstar(jcomp).EQ.-1) go to 100
543  END IF
544 *
545 * Skip merger if an outer binary is fairly perturbed or not hard.
546  40 IF (jcomp.GT.n) THEN
547  IF (gamma(jpair).GT.1.0e-03.OR.eb2.GT.ebh) go to 100
548  END IF
549 *
550 * Estimate relative perturbation at apocentre from actual value.
551  gi = pert*(semi1*(1.0 + ecc1)/rij)**3
552  IF (pert.GT.gmax.OR.gi.GT.0.02) go to 100
553 *
554 * Switch to direct integration for planetary systems if GI > 1D-04.
555  IF (min(body(i1),body(i2)).LT.0.05*bodym) THEN
556  IF (gi.GT.1.0d-04) go to 100
557  END IF
558 *
559 * Ensure the inner semi-major axis is used for subsequent tests.
560  semi = -0.5*body(i)/h(ipair)
561  IF (nmerge.EQ.mmax - 1) THEN
562  IF (nwarn.LT.1000) THEN
563  nwarn = nwarn + 1
564  WRITE (6,41) nmerge
565  41 FORMAT (5x,'WARNING! MERGER LIMIT NMERGE =',i4)
566  END IF
567  go to 100
568  END IF
569 *
570 * Do not allow merger in the inner region of perturbed eccentric orbit.
571  IF (rij.LT.semi1.AND.list(1,i1).GT.0) THEN
572 * Note: moved down from label 30 with 0.1*SEMI1 replacing 2*PMIN.
573  IF (ecc1.GT.0.95.AND.rij.LT.0.1*semi1) THEN
574  go to 100
575  END IF
576  END IF
577 *
578 * -----------------------------------------------------------------------
579 * Form coefficients for stability test (Valtonen, Vistas Ast 32, 1988).
580 * AM = (2.65 + ECC)*(1.0 + BODY(JCOMP)/BODY(I))**0.3333
581 * FM = (2.0*BODY(JCOMP) - BODY(I))/(3.0*BODY(I))
582 * Note that routine NEWTEV in MERGE2 replaces suppressed part.
583 * IF (KZ(19).GE.3) THEN
584 * TM = MIN(TEV(I1),TEV(I2),TEV(JCOMP),TEV(I))
585 * IF (MIN(NAME(I),NAME(JCOMP)).LT.0.AND.TM-TIME.LT.0.2) THEN
586 * GO TO 100
587 * END IF
588 * END IF
589 *
590 * Expand natural logarithm for small arguments.
591 * IF (ABS(FM).LT.0.67) THEN
592 * BM = FM*(1.0 - (0.5 - ONE3*FM)*FM)
593 * ELSE
594 * BM = LOG(1.0D0 + FM)
595 * END IF
596 *
597 * Adopt mass dependent criterion of Harrington (A.J. 82, 753) & Bailyn.
598 * PCRIT = AM*(1.0 + 0.7*BM)*SEMI
599 * -----------------------------------------------------------------------
600 *
601 * Form hierarchical stability ratio (Eggleton & Kiseleva 1995).
602 * QL = BODY(I)/BODY(JCOMP)
603 * Q1 = MAX(BODY(I2)/BODY(I1),BODY(I1)/BODY(I2))
604 * Q3 = QL**0.33333
605 * Q13 = Q1**0.33333
606 * AR = 1.0 + 3.7/Q3 - 2.2/(1.0 + Q3) + 1.4/Q13*(Q3 - 1.0)/(Q3 + 1.0)
607 * EK = AR*SEMI*(1.0D0 + ECC)
608 *
609 * Choose the most dominant triple in case of two binaries.
610  IF (jcomp.GT.n) THEN
611 * Adopt 10% fudge factor with linear dependence on smallest ratio.
612  yfac = 1.0 + 0.1*min(semi2/semi,semi/semi2)
613  ELSE
614  yfac = 1.0
615  END IF
616 *
617 * Prepare inclination evaluation for triple or widest inner binary.
618  IF (jcomp.GT.n) THEN
619 * Ensure widest inner binary (swap is OK for termination or ZARE).
620  IF (semi.LT.semi2) THEN
621  ecc2 = (1.0 - r(jpair)/semi2)**2 +
622  & tdot2(jpair)**2/(body(jcomp)*semi2)
623  ecc = sqrt(ecc2)
624  kpair = ipair
625  ipair = jpair
626  jpair = kpair
627  i1 = 2*ipair - 1
628  i2 = i1 + 1
629  jj = i
630  i = jcomp
631  jcomp = jj
632  semiz = semi2
633  semi2 = semi
634  semi = semiz
635  END IF
636  END IF
637 *
638 * Resolve binary (even perturbed KS not always done).
639  IF (ecc1.LT.1.0) THEN
640  CALL resolv(ipair,1)
641  END IF
642 *
643 * Copy coordinates and velocities to local variables.
644  DO 42 k = 1,3
645  xx(k,1) = x(k,i1)
646  xx(k,2) = x(k,i2)
647  xx(k,3) = x(k,jcomp)
648  vv(k,1) = xdot(k,i1)
649  vv(k,2) = xdot(k,i2)
650  vv(k,3) = xdot(k,jcomp)
651  42 CONTINUE
652 *
653 * Determine the inclination (in radians).
654  CALL inclin(xx,vv,x(1,i),xdot(1,i),angle)
655 *
656 * Employ the basic stability criterion for fast check (ECC1 < 1).
657 * IF (ECC1.LT.1.0) THEN
658 * Q = BODY(JCOMP)/BODY(I)
659 * XFAC = (1.0 + Q)*(1.0 + ECC1)/SQRT(1.0 - ECC1)
660 * PCRIT = 2.8*XFAC**0.4*SEMI*(1.0 - 0.6*ANGLE/TWOPI)
661 * IF (PCRIT.GT.2.0*PMIN) GO TO 100
662 * END IF
663 *
664 * Evaluate the general stability function (Mardling 2008).
665  IF (ecc1.LT.1.0.AND.yfac.LT.1.02) THEN
666  bj = body(jcomp)
667  eout = ecc1
668 * Increase tolerance near sensitive stability boundary (RM 10/2008).
669  IF (eout.GT.0.8) THEN
670  de = 0.5*(1.0 - eout)
671  de = min(de,0.01d0)
672 * Evaluate outer eccentricity derivative due to dominant perturber.
673  IF (jmax.NE.jcomp.AND.pert.GT.gmin) THEN
674  CALL edot(i,jcomp,jmax,semi1,ecc1,eccdot)
675 * Include a small effect of positive derivative over 10 orbits.
676  IF (eccdot.GT.0.0) THEN
677  tk1 = twopi*semi1*
678  & sqrt(semi1/(body(i) + body(jcomp)))
679  de = de - 10.0*min(eccdot*tk1,0.001d0)
680 * WRITE (6,443) ECC1, DE, ECCDOT, TK1, PERT
681 * 443 FORMAT (' EDOT!! E1 DE ED TK G ',
682 * & 2F9.5,1P,3E9.1)
683  END IF
684  END IF
685 * Allow extra tolerance after 1000 tries (suppressed 9/3/12).
686 * IF (NMTRY.GE.1000) DE = MIN(1.0D0 - EOUT,0.02D0)
687  eout = eout - de
688  pmin = semi1*(1.0 - eout)
689  END IF
690  nst = nstab(semi,semi1,ecc,eout,angle,body(i1),body(i2),bj)
691  IF (nst.EQ.0) THEN
692  pcrit = 0.98*pmin*(1.0 - pert)
693  pcr = stability(body(i1),body(i2),body(jcomp),ecc,eout,
694  & angle)
695  pcr = pcr*semi
696 * Specify reduced peri if old criterion < PCRIT/2 (avoids switching).
697  IF (pcr.LT.0.5*pcrit) THEN
698  pcrit = 0.75*pcrit
699  END IF
700  IF (pcrit.LT.yfac*pcr.AND.pert.LT.0.02.AND.
701  & nmtry.LT.10) THEN
702  alph = 360.0*angle/twopi
703  fail = pmin*(1-pert) - yfac*pcr
704  WRITE (6,43) ttot, alph, ecc, ecc1, pmin, fail, pert
705  43 FORMAT (' NEWSTAB T INC EI EO PMIN FAIL PERT ',
706  & f7.1,f7.2,2f8.4,1p,3e10.2)
707  END IF
708 * WRITE (57,444) BODY(I1),(X(K,I1),K=1,3),(XDOT(K,I1),K=1,3)
709 * WRITE (57,444) BODY(I2),(X(K,I2),K=1,3),(XDOT(K,I2),K=1,3)
710 * WRITE (57,444)BODY(JCOMP),(X(K,JCOMP),K=1,3),(XDOT(K,JCOMP),K=1,3)
711 * 444 FORMAT (' ',1P,E14.6,1P,3E18.10,3E14.6)
712  IF (nmtry.GE.1000) THEN
713 * WRITE (6,44) TTOT, NAME(I1), ECC1, EOUT, PCRIT, PMIN
714 * 44 FORMAT (' MARGINAL STAB T NM E1 EOUT PCR PMIN ',
715 * & F7.1,I7,2F8.4,1P,2E10.2)
716  nmtry = 0
717  END IF
718  ELSE
719  nmtry = nmtry + 1
720  pcrit = 1.01*pmin
721  END IF
722  ELSE
723  pcr = stability(body(i1),body(i2),body(jcomp),ecc,ecc1,angle)
724  pcrit = pcr*semi
725  END IF
726 *
727 * Check whether the main perturber dominates the outer component.
728  IF (jmax.NE.jcomp) THEN
729  rij2 = (x(1,jmax) - x(1,jcomp))**2 +
730  & (x(2,jmax) - x(2,jcomp))**2 +
731  & (x(3,jmax) - x(3,jcomp))**2
732  fmax = (body(jmax) + body(jcomp))/rij2
733  IF (fmax.GT.(body(i) + body(jcomp))/rjmin2) go to 100
734  END IF
735 *
736 * Determine time-scale for stability (absolute or approximate).
737  pm1 = pmin*(1.0 - 2.0*pert)
738  CALL tstab(i,ecc1,semi1,pm1,yfac,iterm)
739  IF (iterm.GT.0) go to 100
740 *
741 * Check perturbed stability condition.
742  IF (pmin*(1.0 - pert).LT.yfac*pcrit) go to 100
743 *
744 * Extend active Roche case up to end of look-up time (bug fix 01/09).
745  IF (kstar(i).GE.11.AND.mod(kstar(i),2).NE.0) THEN
746  dt = tev(i) - time
747  IF (dt.LT.10.0*stepx) go to 100
748  tmdis(nmerge+1) = min(time + dt, tmdis(nmerge+1))
749  END IF
750 *
751 * Obtain growth time and inclination for significant eccentricity.
752  IF (semi1.GT.0.0.AND.ecc.GT.0.1) THEN
753  icirc = -1
754  CALL induce(ipair,emax,emin,icirc,tc,angle,tg,edav)
755 * Prevent termination for TC < 2000 in HMDOT but allow small EMAX & DE.
756  IF (kz(27).EQ.2.AND.tc.LT.2000.0) THEN
757  de = abs(emax - ecc)
758  dt = min(tev(i1),tev(i2),tev(i)) - time
759 * Enforce an update next block-step for small remaining times.
760  IF (dt.LT.0.1.AND.mod(kstar(i),2).EQ.0) THEN
761  tev(i1) = time + 0.1
762  tev(i2) = time + 0.1
763  WRITE (6,46) time+toff, name(i1), kstar(i), ecc, tc
764  46 FORMAT (' ENFORCED TEV UPDATE T NM K* E TC ',
765  & f9.2,i8,i4,f8.4,f7.1)
766  END IF
767 * Accept KSTAR = 10, coasting Roche or small eccentricity.
768  IF (mod(kstar(i),2).EQ.0) THEN
769  CALL trflow(ipair,dtr)
770  IF (dtr.LT.step(i)) THEN
771  tev(i) = time + dtr
772  END IF
773  ELSE
774  IF (emax.GT.0.8.OR.de.GT.0.2.OR.dt.LT.0.1) go to 100
775  END IF
776  END IF
777  END IF
778 *
779 * Check circularization and dissipation time (exclude Roche stages).
780  IF (kz(27).EQ.2.AND.kstar(i).LT.10.AND.ecc1.LT.1.0.AND.
781  & name(i).GT.0) THEN
782  ecc0 = ecc
783  CALL decide(ipair,semi,ecc,emax,emin,tc,tg,edav,iq)
784  IF (iq.GT.0.OR.iphase.LT.0) go to 100
785  tk1 = twopi*semi1*sqrt(semi1/(body(i) + body(jcomp)))
786 * IF (TMDIS(NMERGE+1) - TIME.LT.TK1) GO TO 100
787 * EK = EK*(1.0 + ECC)/(1.0 + ECC0)
788  pcrit = pcrit*(1.0 + ecc)/(1.0 + ecc0)
789  END IF
790 *
791 * Perform safety check on radii for case of no circularization.
792  IF (kz(27).EQ.0) THEN
793  IF (semi*(1.0 - ecc).LT.2.0*max(radius(i1),radius(i2))) THEN
794  IF (kz(19).EQ.0) THEN
795  go to 100
796  END IF
797  END IF
798  END IF
799 *
800 * Include rare case of circularizing outer binary.
801  IF (kz(27).EQ.2.AND.jcomp.GT.n.AND.kstar(jcomp).EQ.-2) THEN
802  ecc2 = (1.0 - r(jpair)/semi2)**2 +
803  & tdot2(jpair)**2/(body(jcomp)*semi2)
804  eccj = sqrt(ecc2)
805 * See whether to reduce the look-up time TMDIS (no skip here).
806  CALL decide(jpair,semi2,eccj,emaxj,emin,tc,tg,edav,iq)
807  IF (iphase.LT.0) go to 100
808  END IF
809 *
810 * Check Zare exchange stability criterion and create diagnostics.
811  IF (semi1.GT.0.0) THEN
812  CALL zare(i1,i2,sp)
813  q = body(jcomp)/body(i)
814 * Note inclination is determined by routine INCLIN for ECC < 0.1.
815  IF (sp.LT.1.0.AND.angle.LT.0.17) THEN
816  izare = izare + 1
817  IF (izare.LT.200) THEN
818  WRITE (7,48) ttot, q, ecc, ecc1, semi, pmin, pcrit,
819  & yfac, sp
820  WRITE (7,47) i,jcomp,n,i1,i2,rij,semi1
821  47 FORMAT (' I JCOMP N I1 I2 RIJ A1 ',5i6,1p,2e10.2)
822  CALL flush(7)
823  WRITE (6,48) ttot, q, ecc, ecc1, semi, pmin, pcrit,
824  & yfac, sp
825  48 FORMAT (' ZARE TEST T Q E E1 A PM PCR YF SP ',
826  & f8.2,f5.1,2f7.3,1p,3e9.1,0p,2f6.2)
827  END IF
828  END IF
829  WRITE (73,49) ttot, q, ecc, ecc1, semi, pmin, pcrit,
830  & tg, sp, 360.0*angle/twopi, kstar(i)
831  49 FORMAT (' STAB T Q E E1 A PM PCR TG SP IN K* ',
832  & f8.2,f5.1,2f7.3,1p,4e9.1,0p,f6.2,f7.1,i4)
833  CALL flush(73)
834 * IF (KSTAR(I1).GE.10) TEV(I1) = 1.0E+10
835 * IF (KSTAR(I2).GE.10) TEV(I2) = 1.0E+10
836  END IF
837 *
838 * Specify the final critical pericentre using the fudge factor.
839  pcrit = yfac*pcrit
840 *
841  IF (name(i).GT.n.AND.name(jcomp).GT.n.AND.ecc1.GT.1.0) go to 100
842 * Skip if #JCOMP is a chain c.m. but allow bound double hierarchy.
843  IF (name(jcomp).EQ.0) go to 100
844  IF (ecc1.GT.1.0.AND.min(name(i),name(jcomp)).LT.0) go to 100
845  DO 55 isub = 1,nsub
846  IF (name(jcomp).EQ.names(1,isub)) go to 100
847  55 CONTINUE
848 * Do not allow the formation of a SEPTUPLET.
849 * IF ((NAME(I).LT.-2*NZERO.AND.JCOMP.GT.N).OR.
850 * & NAME(JCOMP).LT.-2*NZERO) GO TO 100
851 *
852 * Include diagnostics for double hierarchy or optional standard case.
853  IF (name(i).LT.0.OR.name(jcomp).LT.0) THEN
854  IF (kz(15).GT.1) THEN
855  which1 = ' MERGE2 '
856  WRITE (6,20) which1, ipair, ttot, h(ipair), r(ipair),
857  & body(i), body(jcomp), pert4, rij, pmin,
858  & eb1/eb, list(1,i1)
859  END IF
860 * Note rare case of two hierarchies merging and identify ghost names.
861  IF (name(i).LT.0.AND.name(jcomp).LT.0) THEN
862  CALL findj(i1,ji,im)
863  j1 = 2*jpair - 1
864  CALL findj(j1,jj,jm)
865  WRITE (6,60) name(i1), name(ji), name(i1+1), name(j1),
866  & name(jj), name(j1+1), ecc, ecc1, semi,
867  & semi1, pmin, pcrit
868  60 FORMAT (' HI MERGE NAM E E1 A A1 PM PC ',
869  & 6i7,2f7.3,1p,4e10.2)
870  END IF
871  ELSE IF (kz(15).GT.1) THEN
872  which1 = ' MERGER '
873  IF (jcomp.GT.n) which1 = ' QUAD '
874  WRITE (6,20) which1, ipair, ttot, h(ipair), r(ipair),
875  & body(i), body(jcomp), pert4, rij, pmin,
876  & eb1/eb, list(1,i1)
877  END IF
878 *
879 * Check for diagnostic output of quadruples.
880  IF (semi1.GT.0.0.AND.jcomp.GT.n.AND.kz(37).GT.0) THEN
881  zmb = body(i) + body(jcomp)
882  tk = days*semi1*sqrt(semi1/zmb)
883  WRITE (89,65) ttot, name(2*ipair-1), name(2*jpair-1),
884  & listq(1), sqrt(ri2), ecc1, eb, eb2, eb1,
885  & tk, pmin, pcrit
886  65 FORMAT (' QUAD# T NAM LQ RI E1 EB EB2 EB1 P1 PM PC ',
887  & f8.1,2i6,i4,f6.2,f8.4,1p,3e12.3,3e9.1)
888  END IF
889 *
890 * Generate a diagnostic file of stable hierarchies (suppressed).
891  IF (ecc1.LT.-1.0) THEN
892  ri = sqrt(ri2)/rc
893  WRITE (80,70) tphys, ri, name(jcomp), ql, q1, ecc, ecc1,
894  & semi, semi1, pcrit/pmin, 360.0*angle/twopi,emax
895  70 FORMAT (f8.1,f5.1,i6,2f6.2,2f6.3,1p,2e10.2,0p,f5.2,f6.1,f6.3)
896  CALL flush(80)
897  END IF
898 *
899 * Copy pair index and set indicator for calling MERGE from MAIN.
900  kspair = ipair
901  iphase = 6
902 *
903 * Save KS indices and delay merger until end of block step.
904  CALL delay(ks2,ks2)
905 *
906  100 IF (iphase.NE.8) jcmax = 0
907 *
908  RETURN
909 *
910  END