Nbody6
 All Files Functions Variables
ksint.f
Go to the documentation of this file.
1  SUBROUTINE ksint(I1)
2 *
3 *
4 * Regularized integration.
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/modes/ eb0(ntmax),zj0(ntmax),ecrit(ntmax),ar(ntmax),
12  & br(ntmax),eosc(4,ntmax),edec(ntmax),tosc(ntmax),
13  & zp(ntmax),es(ntmax),cz(2,ntmax),iosc(ntmax),
14  & namec(ntmax)
15  common/slow0/ range,islow(10)
16  common/gamdot/ dgam
17  REAL*8 ui(4),uidot(4),xi(6),vi(6),fp(6),fd(6)
18  LOGICAL iq
19  SAVE iterm
20 *
21 *
22 * Set second component, pair index & c.m. index.
23  i2 = i1 + 1
24  ipair = kvec(i1)
25  i = n + ipair
26 *
27 * Define perturber membership & inverse c.m. mass.
28  nnb0 = list(1,i1)
29  bodyin = 1.0/body(i)
30 *
31 * Check for further unperturbed motion.
32  IF (nnb0.EQ.0.AND.h(ipair).LT.0.0) THEN
33 * Include possible eccentricity modulation of hierarchical binary.
34  IF (name(i).LT.0) THEN
35  im = 0
36  DO 1 k = 1,nmerge
37  IF (namem(k).EQ.name(i)) im = k
38  1 CONTINUE
39  IF (im.GT.0.AND.time.GT.tmdis(im)) THEN
40  semi = -0.5*body(i)/h(ipair)
41  ecc2 = (1.0 - r(ipair)/semi)**2 +
42  & tdot2(ipair)**2/(body(i)*semi)
43  rp = semi*(1.0 - sqrt(ecc2))*(1.0 - 2.0*gamma(ipair))
44  IF (rp.LT.r0(ipair)) THEN
45  WRITE (77,4) name(i1), sqrt(ecc2), rp, r0(ipair)
46  4 FORMAT (' TMDIS TERM NM E RP R0',i7,1p,3e10.2)
47  CALL flush(77)
48  instab = instab + 1
49  go to 90
50  ELSE IF (kz(27).EQ.2) THEN
51  CALL eccmod(i,iterm)
52 * Update time on termination to prevent looping.
53  IF (iterm.GT.0) THEN
54  t0(i1) = time
55  go to 90
56  END IF
57  END IF
58 * Check any inner and outer circularizing binary using NAMEM & NAMEG.
59  DO 3 k = 1,nchaos
60  IF (namec(k).EQ.nzero - namem(im)) THEN
61 * Update unperturbed binary if T - TOSC > 10 Myr (cf. IMPACT & DECIDE).
62  IF ((time - tosc(k))*tstar.GT.10.0) THEN
63  t0(i1) = time
64  go to 90
65  END IF
66  END IF
67  IF (namec(k).EQ.nameg(im)) THEN
68  IF ((time - tosc(k))*tstar.GT.10.0) THEN
69  t0(i1) = time
70  go to 90
71  END IF
72  END IF
73  3 CONTINUE
74  END IF
75  END IF
76  dt1 = step(i1)
77  CALL unpert(ipair)
78 *
79 * Check updating of unperturbed relativistic KS binary.
80  IF (kz(11).NE.0.AND.list(1,i1).EQ.0) THEN
81  CALL brake4(i1,i2,dt1)
82  IF (iphase.LT.0) go to 100
83  END IF
84 *
85 * Try re-initialize chain WD/BH system after dormant KS (#11 only).
86  IF (kz(11).NE.0.AND.nch.EQ.0.AND.list(1,i1).GT.0) THEN
87  IF (min(kstar(i1),kstar(i2)).GE.10) THEN
88  semi = -0.5*body(i)/h(ipair)
89  IF (semi.GT.10.0*rmin) go to 100
90  WRITE (6,222) time+toff, name(jclose), kstar(i1),
91  & kstar(i2), list(1,i1), gamma(ipair),
92  & semi, r(ipair)
93  222 FORMAT (' ACTIVATE CHAIN T NMJ K* NP G A R ',
94  & f9.1,i7,3i4,1p,3e10.2)
95  kspair = ipair
96 * Restore unperturbed motion from BRAKE4 (NP = 0 fixes some problem).
97  IF (gamma(ipair).LT.1.0d-10) THEN
98  jclose = 0
99  list(1,i1) = 0
100  END IF
101  ks2 = 0
102 * Include case of binary as dominant perturber.
103  IF (jclose.GT.n) THEN
104  ks2 = jclose - n
105  IF (ks2.GT.kspair) ks2 = ks2 - 1
106  jcomp = jclose
107  jp = jclose - n
108  WRITE (6,223) kspair, ks2, jclose, gamma(jp)
109  223 FORMAT (' BINARY PERT KSP KS2 JCLOSE GJP ',
110  & 2i4,i7,1p,e10.2)
111  ELSE
112 * Avoid JCOMP > N & JCLOSE < N for spurious CALL KSTERM in DELAY.
113  jcomp = 0
114  END IF
115  iphase = 8
116 * Save KS parameters until end of block-step (JCMAX=0: no extra pert).
117  CALL delay(1,ks2)
118  jcmax = 0
119  END IF
120  END IF
121  go to 100
122  END IF
123 *
124 * Perform KS prediction of U, UDOT & H.
125  CALL kspred(ipair,i1,i,bodyin,ui,uidot,xi,vi)
126 *
127 * Obtain the perturbing force & derivative.
128  CALL kspert(i1,nnb0,xi,vi,fp,fd)
129 *
130 * Save old radial velocity & relative perturbation and set new GAMMA.
131  rdot = tdot2(ipair)
132 * G0 = GAMMA(IPAIR)
133  gi = sqrt(fp(1)**2 + fp(2)**2 + fp(3)**2)*r(ipair)**2*bodyin
134  gamma(ipair) = gi
135 *
136 * Apply the Hermite corrector.
137  CALL kscorr(ipair,ui,uidot,fp,fd,td2,tdot4,tdot5,tdot6)
138 *
139 * Increase regularization time-step counter and update the time.
140  nstepu = nstepu + 1
141  t0(i1) = time
142 * Check for early return during termination (called from KSTERM).
143  IF (iphase.NE.0) go to 100
144 *
145 * Define useful scalars.
146  ri = r(ipair)
147  hi = h(ipair)
148 *
149 * Initialize termination indicator and check for large perturbation.
150  iq = .false.
151  IF (gi.GT.0.25) go to 2
152  IF (gi.LT.0.03) THEN
153  jcomp = 0
154  go to 20
155  END IF
156  CALL flyby(i,iterm)
157  IF (iterm.EQ.0.AND.kstar(i).LT.0) THEN
158  iq = .true.
159  END IF
160  IF (iterm.EQ.1) THEN
161 * Delay chain regularization search until near end of block-step.
162  IF (time + step(i1).GT.tblock) THEN
163  CALL impact(i)
164  IF (iphase.GT.0) go to 100
165  END IF
166  ELSE IF (iterm.EQ.2) THEN
167  iq = .true.
168  go to 20
169  ELSE
170  go to 20
171  END IF
172 *
173 * Find the dominant body for large perturbations.
174  2 s = 4.0*step(i)
175  fmax = body(i)/ri**2
176 * Initialize JCOMP for prediction and optional diagnostics in KSTERM.
177  jcomp = 0
178  DO 10 l = 2,nnb0+1
179  j = list(l,i1)
180 * Only search bodies within twice the c.m. time-step.
181  IF (step(j).GT.s) go to 10
182 * Compare strong perturber and either component with current pair.
183  DO 5 k = i1,i2
184  rij2 = (x(1,j) - x(1,k))**2 + (x(2,j) - x(2,k))**2 +
185  & (x(3,j) - x(3,k))**2
186  IF (body(j) + body(k).GT.rij2*fmax) jcomp = j
187  5 CONTINUE
188  10 CONTINUE
189 *
190 * Set termination if strong perturber <= N forms dominant pair.
191  IF (jcomp.GT.0.OR.gi.GT.1.0) THEN
192 * Check optional binary search.
193 * IF (KZ(4).GT.0) THEN
194 * DGAM = GI - G0
195 * K = KZ(4)
196 * CALL EVOLVE(IPAIR,K)
197 * END IF
198  IF (jcomp.LE.n.OR.gi.GT.2.0) iq = .true.
199  END IF
200 *
201 * Check termination of hyperbolic encounter (R > R0 or R > 2*RMIN).
202  20 IF (hi.GT.0.0d0.AND.name(i).GT.0) THEN
203  IF ((ri.GT.r0(ipair).AND.gi.GT.gmax).OR.ri.GT.2.0*rmin.OR.
204  & (gi.GT.0.5.AND.td2.GT.0.0)) THEN
205 * Skip termination delay in case of velocity kick (cf. routine KSTERM).
206  IF (hi.LT.100.0.OR.gi.GT.0.1.OR.ri.GT.5.0*rmin) THEN
207  iq = .true.
208  END IF
209  END IF
210  END IF
211 *
212 * Choose basic regularized step using binding energy or lower limit.
213  IF (abs(hi).GT.eclose) THEN
214  w1 = 0.5/abs(hi)
215  ELSE
216  w1 = r0(ipair)*bodyin
217  w2 = 0.5/abs(hi)
218  w1 = min(w1,w2)
219  IF (ri.GT.r0(ipair)) w1 = w1*r0(ipair)/ri
220 * Maximum square step for soft binaries & weak hyperbolic pairs.
221  IF (hi.LT.0.0d0) THEN
222 * Include case of hard binary with massive components or merger.
223  w2 = -0.5/hi
224  IF (name(i).LT.0) THEN
225  w1 = ri*bodyin
226  END IF
227  w1 = min(w2,w1)
228  END IF
229  END IF
230 *
231 * Include perturbation factor in predicted step.
232  IF (gi.LT.0.0005) THEN
233 * Use second-order expansion of cube root for small perturbations.
234  w3 = 333.3*gi
235  w2 = sqrt(w1)/(1.0 + w3*(1.0 - w3))
236  ELSE
237  w3 = 1.0 + 1000.0*gi
238  w2 = sqrt(w1)/w3**0.3333
239  END IF
240 *
241 * Form new regularized step (include inertial factor).
242  dtu = etau*w2
243  dtu = min(1.2*dtau(ipair),dtu)
244 *
245 * Include convergence criterion DH = H'*DTU + H''*DTU**2/2 = 0.001*|H|.
246  IF (gi.GT.1.0d-04) THEN
247  dh = 1.0e-03*max(abs(h(ipair)),eclose)
248  xf = 2.0*dh/abs(hdot2(ipair))
249  yf = hdot(ipair)/hdot2(ipair)
250  dtu1 = sqrt(xf + yf**2) - abs(yf)
251  dtu = min(dtu1,dtu)
252  END IF
253 *
254 * Check pericentre step reduction for perturbed spiral.
255  IF (kstar(i).EQ.-2.AND.td2.LT.0.0) THEN
256  semi = -0.5*body(i)/hi
257  IF (ri.LT.semi) THEN
258 * Predict radial velocity for step DTU (note scaled coefficients).
259  rd = ((one3*tdot5*dtu + tdot4)*dtu +
260  & tdot3(ipair))*dtu + 2.0*td2
261 * Adopt pericentre step with 1% safety factor (small TDOT4 is OK).
262  IF (rd.GT.0.0) THEN
263  a2 = 0.5*tdot3(ipair)/tdot4
264  dtu1 = sqrt(a2**2 - 2.0*td2/tdot4) - a2
265  dtu1 = 1.01*max(dtu1,1.0d-10)
266  dtu = min(dtu1,dtu)
267  END IF
268  END IF
269  END IF
270 *
271 * Reset reference energy and generate new Stumpff coefficients.
272  h0(ipair) = h(ipair)
273  30 z = -0.5d0*h(ipair)*dtu**2
274  CALL stumpf(ipair,z)
275  z5 = sf(6,ipair)
276  z6 = sf(7,ipair)
277  dt12 = one12*dtu*z6
278 *
279 * Convert to physical time units modified by Stumpff coefficients.
280  step(i1) = (((((tdot6*dt12 + tdot5*z5)*0.2*dtu + 0.5d0*tdot4)*dtu
281  & + tdot3(ipair))*one6*dtu + td2)*dtu + ri)*dtu
282 *
283 * Ensure that regularized step is smaller than the c.m. step.
284  IF (step(i1).GT.step(i).AND.hi.LT.0) THEN
285  dtu = 0.5d0*dtu
286  go to 30
287  END IF
288  dtau(ipair) = dtu
289 *
290 * See whether the KS slow-down procedure is activated.
291  imod = kslow(ipair)
292  IF (imod.GT.1) THEN
293  zmod = float(islow(imod))
294  step(i1) = zmod*step(i1)
295  END IF
296 *
297 * Check diagnostics print option.
298  IF (kz(10).GE.3) THEN
299  WRITE (6,40) ipair, time+toff, h(ipair), ri, dtau(ipair),
300  & gi, step(i1), list(1,i1), imod
301  40 FORMAT (3x,'KS MOTION',i6,2f10.4,1p,4e10.2,2i4)
302  END IF
303 *
304 * Employ special termination criterion in merger case.
305  IF (name(i).LT.0) THEN
306 * Terminate if apocentre perturbation > 0.25 (R > SEMI) or GI > 0.25.
307  IF (hi.LT.0.0) THEN
308  semi = -0.5*body(i)/hi
309 * ECC2 = (1.0 - RI/SEMI)**2 + TDOT2(IPAIR)**2/(BODY(I)*SEMI)
310 * A0 = SEMI*(1.0 + SQRT(ECC2))/RI
311 * Replace eccentricity calculation with typical value.
312  a0 = 1.5*semi/ri
313  ga = gi*a0*a0*a0
314  IF (ga.GT.0.25.AND.ri.GT.semi) iq = .true.
315  IF (ri.GT.20*rmin.AND.nnb0.GT.0.8*list(1,i)) iq = .true.
316  IF (gi.GT.0.1.AND.ri.GT.rmin) iq = .true.
317  IF (gi.GT.0.25) iq = .true.
318 * Include extra condition for planet case.
319  IF (min(body(i1),body(i2)).LT.0.05*bodym) THEN
320  IF (gi.GT.2.0d-04) iq = .true.
321  END IF
322  ELSE
323  IF (td2.GT.0.0.AND.(gi.GT.gmax.OR.ri.GT.rmin)) iq = .true.
324  END IF
325  IF (.NOT.iq) go to 60
326  END IF
327 *
328 * Delay termination until end of block for large perturbations.
329  IF (iq) THEN
330  dtr = tblock - time
331 * WRITE (6,45) IPAIR, TTOT, GI, RI, DTR, STEP(I1)
332 * 45 FORMAT (' TERM TEST KS T G R DTR DT ',
333 * & I4,F10.4,F7.3,1P,E10.2,2E9.1)
334  IF (dtr.LT.step(i1)) go to 90
335  END IF
336 *
337 * Check standard termination criterion (suppress on IQ = .true.).
338  IF (ri.GT.r0(ipair).AND.ri.GT.2.0*rmin.AND..NOT.iq) THEN
339 * Include termination for rare tidal capture starting at pericentre.
340  IF (kstar(i).LT.0.AND.ri.GT.5.0*rmin) go to 90
341 * Impose a limit using size of neighbour sphere (100*R > 0.80*RS).
342  IF (ri.GT.8.0d-03*rs(i).AND.gi.GT.1.0d-04) go to 90
343 * See whether termination can be delayed for sufficient perturbers.
344  IF (nnb0.LT.0.80*list(1,i).AND.gi.LT.0.1) go to 60
345 * Check updating of R0 for newly hardened binary orbit.
346  IF (hi.LT.-eclose) THEN
347  semi = -0.5*body(i)/hi
348  r0(ipair) = max(rmin,2.0d0*semi)
349  r0(ipair) = min(r0(ipair),5.0*rmin)
350  go to 70
351  END IF
352 *
353 * IF (KZ(4).GT.0.AND.GI.GT.GPRINT(1)) THEN
354 * DGAM = GI - G0
355 * K = KZ(4)
356 * DO 55 L = 2,K
357 * IF (GI.LT.GPRINT(L)) THEN
358 * CALL EVOLVE(IPAIR,L-1)
359 * GO TO 90
360 * END IF
361 * 55 CONTINUE
362 * CALL EVOLVE(IPAIR,K)
363 * END IF
364  go to 90
365  END IF
366 *
367 * End integration cycle for hyperbolic motion.
368  60 IF (hi.GE.0.0d0) THEN
369  IF (rdot*td2.LT.0.0d0) THEN
370 * Determine pericentre for hyperbolic two-body motion.
371  semi = -0.5d0*body(i)/hi
372  ecc2 = (1.0 - ri/semi)**2 + tdot2(ipair)**2/(body(i)*semi)
373  qperi = semi*(1.0d0 - sqrt(ecc2))
374  dmin2 = min(dmin2,qperi)
375 *
376 * Check optional tidal interaction or stellar collision.
377  IF (kz(19).GE.3.AND.kstar(i).LT.10) THEN
378  vinf = sqrt(2.0*hi)*vstar
379  ks1 = kstar(i1)
380  ks2 = kstar(i2)
381  rx = max(radius(i1),radius(i2))
382 * Determine maximum periastron factor for capture (VINF in km/sec).
383  IF (kz(27).LE.2) THEN
384  rfac = rpmax2(radius(i1),radius(i2),body(i1),
385  & body(i2),ks1,ks2,vinf)
386  rcap = rfac*rx
387  ELSE
388  dv = sqrt(2.0*hi)
389 * Note that Quinlan & Shapiro function returns actual distance.
390  rcap = rpmax(body(i1),body(i2),vstar,dv,qperi)
391  END IF
392  IF (qperi.LT.5.0*rx) THEN
393  WRITE (54,54) ttot, name(i1), name(i2), ks1,
394  & ks2, vinf, rcap*su, rx*su, qperi*su
395  54 FORMAT (' CLOSE T NAM K* VINF RCAP RX QP ',
396  & f7.1,2i6,2i4,f6.2,3f6.1)
397  END IF
398  IF (qperi.LT.rcap) THEN
399  j1 = i1
400  IF (radius(i2).GT.radius(i1)) j1 = i2
401  fac = 0.5*body(i)/body(j1)
402 * Adopt collision criterion of Kochanek (Ap.J. 385, 604, 1992).
403  IF (kz(27).LE.2) THEN
404  rcoll = 1.7*fac**0.3333*radius(j1)
405  ELSE
406  rcoll = 6.0*body(i)/clight**2
407  END IF
408  WRITE (55,58) ttot, ipair, name(i1), name(i2),
409  & ks1, ks2, kstar(i), vinf
410  58 FORMAT (' RPMAX: T KS NAM K* VINF ',
411  & f7.1,i5,2i6,3i4,f7.2)
412  WRITE (55,59) sqrt(ecc2), hi, r(ipair), semi,
413  & qperi, body(i1), body(i2),
414  & body(i)*zmbar
415  59 FORMAT (' E H R A QP BODY MT ',
416  & f9.5,1p,6e10.2,0p,f6.1)
417  ri2 = 0.0
418  vi2 = 0.0
419  DO 61 k = 1,3
420  ri2 = ri2 + (x(k,i) - cmr(k))**2
421  vi2 = vi2 + xdot(k,i)**2
422  61 CONTINUE
423  WRITE (55,62) sqrt(ri2)/rc, sqrt(vi2)*vstar,
424  & rhod, radius(i1)*su, radius(i2)*su,
425  & rcap, radius(j1)/qperi, rcoll/qperi
426  62 FORMAT (' r/RC V* <C> R* RCAP R1/QP RCOLL/QP ',
427  & 2f5.1,3f6.1,3f5.1)
428  CALL flush(55)
429  IF (qperi.LT.rcoll) THEN
430 * Obtain KS variables at pericentre before merging into one body.
431  CALL ksperi(ipair)
432  kspair = ipair
433  iqcoll = -2
434  CALL cmbody(qperi,2)
435  ELSE IF (kstar(i).GE.0.AND.kz(27).GT.0) THEN
436  CALL kstide(ipair,qperi)
437  END IF
438  END IF
439 * Check options for artificial collisions.
440  ELSE IF (kz(27).EQ.-1.AND.kz(13).LT.0) THEN
441  rfac = 2.0
442  IF (qperi.LT.rfac*max(radius(i1),radius(i2))) THEN
443  j1 = i1
444  IF (radius(i2).GT.radius(i1)) j1 = i2
445  fac = 0.5*body(i)/body(j1)
446 * Adopt collision criterion of Kochanek (Ap.J. 385, 604, 1992).
447  rcoll = 1.7*fac**0.3333*radius(j1)
448  IF (qperi.LT.rcoll) THEN
449  CALL touch(ipair,i1,i2,rcoll)
450  END IF
451  END IF
452  END IF
453  END IF
454  go to 100
455  END IF
456 *
457 * Check perturbation threshold (H < 0 & GAMMA > GMAX).
458 * IF (KZ(4).EQ.0.OR.G0.LT.GMAX) GO TO 70
459 *
460 * K = KZ(4)
461 * DO 65 L = 1,K
462 * IF ((G0 - GPRINT(L))*(GI - GPRINT(L)).LE.0.0) THEN
463 * IF (L.EQ.1) THEN
464 * W1 = -0.5*BODY(I)/HI
465 * W2 = W1*BODYIN
466 * TK = TWOPI*W1*SQRT(W2)
467 * END IF
468 *
469 * Estimate smallest permitted output interval at new level.
470 * DTCRIT = TK*ORBITS(L)
471 * IF (TIME - TLASTB(L).LT.DTCRIT) GO TO 70
472 * DGAM = GI - G0
473 * CALL EVOLVE(IPAIR,L)
474 * GO TO 70
475 * END IF
476 * 65 CONTINUE
477 *
478 * Check for partial reflection during approach (NB! only IMOD = 1).
479 * 70 IF (GI.LT.GMIN.AND.TD2.LT.0.0D0) THEN
480 * Skip apocentre position itself.
481 * IF (RDOT.LT.0.0D0.AND.IMOD.EQ.1) THEN
482 * IF (KZ(25).GT.0) CALL FREEZE(IPAIR)
483 * GO TO 100
484 * END IF
485 * END IF
486 *
487 * Determine new perturbers for binary at apocentre turning point.
488  70 IF (rdot*td2.GE.0.0d0) go to 100
489  semi = -0.5d0*body(i)/hi
490 *
491 * Check minimum two-body separation just after pericentre.
492  IF (rdot.LT.0.0d0) THEN
493 * Obtain pericentre by Mikkola's algorithm (GAMMA < 0.001).
494  IF (gi.LT.0.001) THEN
495  CALL peri(ui,uidot,ri,body(i1),body(i2),qperi)
496  ELSE
497  qperi = ri
498  END IF
499  dmin2 = min(dmin2,qperi)
500 *
501 * Check optional tidal interaction or stellar collision (skip merger).
502  IF (kz(19).GE.3.AND.kstar(i).LE.10.AND.name(i).GT.0) THEN
503  rfac = 5.0
504  IF (kz(27).LE.2) THEN
505  IF (kz(27).EQ.1) rfac = 4.0
506  rx = rfac*max(radius(i1),radius(i2))
507  ELSE
508  rx = rpmin(body(i1),body(i2),vstar,hi,qperi)
509  END IF
510  IF (qperi.LT.rx) THEN
511  j1 = i1
512  IF (radius(i2).GT.radius(i1)) j1 = i2
513  fac = 0.5*body(i)/body(j1)
514 * Adopt collision criterion of Kochanek (Ap.J. 385, 604, 1992).
515  IF (kz(27).LE.2) THEN
516  rcoll = 1.7*fac**0.3333*radius(j1)
517  ELSE
518  rcoll = 6.0*body(i)/clight**2
519  END IF
520  IF (qperi.LT.rcoll) THEN
521 * Obtain KS variables at pericentre before merging into one body.
522  CALL ksperi(ipair)
523  kspair = ipair
524  iqcoll = -2
525  CALL cmbody(qperi,2)
526  ELSE IF (kstar(i).GE.0) THEN
527 * Distinguish between sequential, standard and GR circularization.
528  IF (kz(27).EQ.1) THEN
529  icirc = 1
530  tc = 0.0
531  ELSE IF (kz(27).EQ.2.AND.kstar(i).LT.10) THEN
532  ecc2 = (1.0 - ri/semi)**2 +
533  & tdot2(ipair)**2/(body(i)*semi)
534  ecc = sqrt(ecc2)
535  icirc = 0
536  CALL tcirc(qperi,ecc,i1,i2,icirc,tc)
537  ELSE
538  icirc = 1
539  tc = 0.0
540  END IF
541  IF (kstar(i).GE.10) icirc = 0
542 * Skip tidal effects for circularization time above 100 Myr (07/08).
543  IF (icirc.GT.0.AND.kz(27).GT.0.AND.
544  & tc.LT.100.0) THEN
545  CALL kstide(ipair,qperi)
546  END IF
547  END IF
548  END IF
549 * Check for perturbed spiral or chaos case (skip collision).
550  IF (kstar(i).EQ.-2.AND.iphase.EQ.0) THEN
551  CALL spiral(ipair)
552  ELSE IF (kstar(i).EQ.-1.AND.iphase.EQ.0) THEN
553  CALL kstide(ipair,qperi)
554  END IF
555 * Check options for artificial collisions.
556  ELSE IF (kz(27).EQ.-1.AND.kz(13).LT.0) THEN
557  rfac = 2.0
558  IF (qperi.LT.rfac*max(radius(i1),radius(i2))) THEN
559  j1 = i1
560  IF (radius(i2).GT.radius(i1)) j1 = i2
561  fac = 0.5*body(i)/body(j1)
562 * Adopt collision criterion of Kochanek (Ap.J. 385, 604, 1992).
563  rcoll = 1.7*fac**0.3333*radius(j1)
564  IF (qperi.LT.rcoll) THEN
565  CALL touch(ipair,i1,i2,rcoll)
566  END IF
567  END IF
568  END IF
569  go to 100
570  END IF
571 *
572 * Save maximum separation of persistent binary.
573  rmax = max(rmax,ri)
574 *
575 * Check binary reference radius or merger stability criterion.
576  IF (name(i).GT.0) THEN
577 * Update termination length scale in case of initial soft binary.
578  eb = body(i1)*body(i2)*hi*bodyin
579  IF (eb.LT.ebh) r0(ipair) = max(rmin,2.0*semi)
580  ELSE
581  ecc2 = (1.0 - ri/semi)**2 + tdot2(ipair)**2/(body(i)*semi)
582  ecc = sqrt(ecc2)
583  rp = semi*(1.0 - ecc)*(1.0 - 2.0*gi)
584 * Find merger index.
585  im = 0
586  DO 72 k = 1,nmerge
587  IF (namem(k).EQ.name(i)) im = k
588  72 CONTINUE
589 * Exclude inner planets from the general stability test.
590  IF (min(body(i1),body(i2)).LT.0.05*bodym) THEN
591  IF (rp.LT.r0(ipair)) go to 90
592  END IF
593 * Include optional Kozai diagnostics (large EMAX) & circularization.
594  IF (kz(42).GT.0) THEN
595  CALL kozai(ipair,im,ecc,semi,iterm)
596 * Check merger termination and perform circularization or collision.
597  IF (iterm.GT.0) THEN
598  iphase = 7
599  kspair = ipair
600  CALL reset
601  IF (iterm.GT.1) CALL circ
602  go to 100
603  END IF
604  END IF
605 * Assess the stability inside critical pericentre (safety factor 1.04).
606  IF (rp.LT.1.04*r0(ipair)) THEN
607 * Note: assessment needs to use same eccentricity as for acceptance.
608  CALL assess(ipair,im,ecc,semi,iterm)
609  IF (iterm.GT.0) THEN
610  instab = instab + 1
611  go to 90
612  END IF
613  END IF
614 * Check possible eccentricity modulation or t_circ update.
615  IF (im.GT.0.AND.(time.GT.tmdis(im).OR.
616  & tmdis(im).GT.1.0d+06)) THEN
617  IF (kz(27).EQ.2) THEN
618  CALL eccmod(i,iterm)
619  IF (iterm.GT.0) THEN
620 * WRITE (6,76) RP, R0(IPAIR)
621 * 76 FORMAT (' ECCMOD TERM RP R0 ',1P,2E10.2)
622  go to 90
623  END IF
624 * Consider both inner and possible outer circularizing binary.
625  DO 78 k = 1,nchaos
626  IF (namec(k).EQ.nzero - namem(im).AND.
627  & kstarm(im).EQ.-2) THEN
628 * Update unperturbed binary if T - TOSC > 10 Myr (cf. IMPACT & DECIDE).
629  IF ((time - tosc(k))*tstar.GT.10.0) go to 90
630  END IF
631  IF (namec(k).EQ.nameg(im).AND.
632  & kstarm(im).EQ.-2) THEN
633  IF ((time - tosc(k))*tstar.GT.10.0) go to 90
634  END IF
635 * Note: perturbed binary is treated if pericentre before next IMPACT.
636  78 CONTINUE
637  END IF
638  END IF
639  END IF
640 *
641 * Produce diagnostics for any circularizing perturbed binary.
642  IF (kstar(i).EQ.-2.AND.gi.GT.0.01) THEN
643  ecc = ri/semi - 1.0
644  qps = semi*(1.0 - ecc)/max(radius(i1),radius(i2))
645  zm = body(i1)/body(i2)
646  WRITE (21,80) ttot, name(i1), name(i2), kstar(i1), kstar(i2),
647  & list(1,i1), qps, zm, gi, ecc, semi
648  80 FORMAT (' PERT SPIRAL T NAM K* NP QP/S M1/M2 G E A ',
649  & f11.4,2i6,3i3,2f5.1,2f8.4,1p,e10.2)
650  CALL flush(21)
651  END IF
652 *
653 * See whether KS slow-down procedure should be (re)-checked (no Chaos).
654  IF (kz(26).GT.0.AND.kstar(i).GE.0) THEN
655  kmod = range*gmin/max(gi,1.0d-10)
656  IF (kmod.GT.1.OR.imod.GT.1) THEN
657  CALL ksmod(ipair,kmod)
658  IF (kmod.LT.0) go to 100
659  go to 82
660  END IF
661  END IF
662 *
663 * Set approximate value of next period with perturbation included.
664  tk = twopi*semi*sqrt(semi*bodyin)*(1.0 + gi)
665  IF (imod.GT.1) THEN
666  tk = zmod*tk
667  END IF
668 *
669 * Use old perturber list if next apocentre is before the c.m. step.
670  IF (time + tk.LT.t0(i) + step(i)) THEN
671  go to 100
672  END IF
673 *
674 * Select new perturbers (J = N adopted for unperturbed Chaos).
675  82 CALL kslist(ipair)
676 *
677 * Check rectification of chaotic spiral at start of unperturbed motion.
678  IF (kstar(i).EQ.-2.AND.list(1,i1).EQ.0) THEN
679  dmr = 0.d0
680  CALL chrect(ipair,dmr)
681  IF (iphase.LT.0) go to 100
682  ELSE
683  CALL ksrect(ipair)
684  END IF
685 *
686 * See whether a massive BH subsystem can be selected.
687  IF (kz(11).NE.0.AND.nch.EQ.0.AND.body(i).GT.10.0*bodym.AND.
688  & semi.LT.rmin.AND.list(1,i1).LE.5.AND.name(i).GT.0) THEN
689 *
690  IF (semi.GT.0.1*rmin) go to 88
691 * Check optional BH condition (prevents mass-loss complications).
692  IF (kz(11).LE.-2) THEN
693  IF (kstar(i1).NE.14.OR.kstar(i2).NE.14) go to 88
694  END IF
695 *
696  rcr2 = rmin2*body(i)/(10.0*bodym)
697  nnb1 = list(1,i1) + 1
698  rx2 = 1000.0
699  jclose = 0
700  DO 85 l = 2,nnb1
701  j = list(l,i1)
702  rij2 = (x(1,i) - x(1,j))**2 + (x(2,i) - x(2,j))**2 +
703  & (x(3,i) - x(3,j))**2
704  rd = (x(1,i)-x(1,j))*(xdot(1,i)-xdot(1,j)) +
705  & (x(2,i)-x(2,j))*(xdot(2,i)-xdot(2,j)) +
706  & (x(3,i)-x(3,j))*(xdot(3,i)-xdot(3,j))
707  IF (rij2.LT.rcr2.AND.rd.LT.0.0.OR.
708  & (rij2.LT.4.0*semi**2)) THEN
709  IF (rij2.LT.rx2) THEN
710  rx2 = rij2
711  rrd = rd
712  jclose = j
713  END IF
714  END IF
715  85 CONTINUE
716  IF (jclose.GT.0) THEN
717  IF (name(jclose).LE.0) go to 88
718  rx = sqrt(rx2)
719 * Limit energy of triple system (< 50*EBH) using radial velocity.
720  zmu = body(i)*body(jclose)/(body(i) + body(jclose))
721  ebt = eb + zmu*(0.5*(rrd/rx)**2-(body(i)+body(jclose)/rx))
722  IF (ebt.GT.50.0*ebh) go to 88
723  WRITE (6,86) ttot, name(jclose), list(1,i1), step(i),
724  & step(jclose), semi, rx, gi, ebt
725  86 FORMAT (' NEW CHAIN T NMJ NP STEPI STEPJ A RIJ GI EBT ',
726  & f9.3,i6,i4,1p,6e10.2)
727  CALL flush(6)
728 * Initiate chain regularization directly (B-B or B-S: see IMPACT).
729  jcomp = jclose
730  jcmax = 0
731  kspair = ipair
732  iphase = 8
733  ebch0 = eb
734 * Distinguish between case of single and binary intruder.
735  IF (jclose.LE.n) THEN
736  ks2 = 0
737  ELSE
738  ks2 = jclose - n
739 * Reduce c.m. body location and pair index if IPAIR terminated first.
740  IF (ks2.GT.ipair) jclose = jclose - 1
741  IF (ks2.GT.ipair) ks2 = ks2 - 1
742  END IF
743 * Initialize new ARchain.
744  CALL delay(1,ks2)
745  go to 100
746  END IF
747  END IF
748 *
749 * Check optional search criterion for multiple encounter or merger.
750  88 IF (kz(15).GT.0.AND.step(i).LT.dtmin) THEN
751  CALL impact(i)
752  END IF
753  go to 100
754 *
755 * Terminate regularization of current pair (IPAIR set in KSPAIR).
756  90 kspair = ipair
757 * Set indicator for calling KSTERM in MAIN (permits phase overlay).
758  iphase = 2
759 * Check case of hierarchical binary.
760  IF (name(i).LT.0) iphase = 7
761 *
762  100 RETURN
763 *
764  END