Nbody6
 All Files Functions Variables
chain.f
Go to the documentation of this file.
1  SUBROUTINE chain(ISUB)
2 *
3 *
4 * Perturbed chain regularization.
5 * -------------------------------
6 *
7 * Method of Mikkola & Aarseth, Celestial Mechanics 47, 375.
8 * .........................................................
9 *
10  include 'commonc.h'
11  include 'common2.h'
12  REAL*8 g0(3),y(nmx8),r2(nmx,nmx),ksch
13  INTEGER ij(nmx),iold(nmx)
14  LOGICAL check,kslow,kcoll,stopb,icase
15  common/slow1/ tk2(0:nmx),ejump,ksch(nmx),kslow,kcoll
16  common/slow2/ stepl,stopb
17  common/chreg/ timec,tmax,rmaxc,cm(10),namec(6),nstep1,kz27,kz30
18  common/clump/ bodys(nmx,5),t0s(5),ts(5),steps(5),rmaxs(5),
19  & names(nmx,5),isys(5)
20  common/intfac/ lx,le,lp,lv,lt,j10,nhalf2
21  common/echain/ ech
22  common/cpert/ rgrav,gpert,ipert,npert
23  common/calls/ tpr,tkpr,step,ider,icall,nfn,nreg,iter,imcirc
24  common/ccoll2/ qk(nmx4),pk(nmx4),rij(nmx,nmx),SIZE(nmx),vstar1,
25  & ecoll1,rcoll,qperi,istar(nmx),icoll,isync,ndiss1
26  common/bssave/ ep(4),dsc,facm,tfac,itfac,jc
27  common/ebsave/ ebs
28  common/ksave/ k1,k2
29  common/slow3/ gcrit,kz26
30  common/swcall/ ncall
31  EXTERNAL chmod
32  SAVE
33 *
34 *
35 * Main variables for chain regularization
36 * ***************************************
37 *
38 * -------------------------------------------------------------------
39 * CHTIME Local physical time (from CHAIN code).
40 * CM(1-7) Coordinates, velocities & total mass of system.
41 * CM(8) Total energy of N-body system (copied in routine CHINIT).
42 * CM(9) Binary energy change (copied to CHCOLL in CHTERM).
43 * ECH Total energy of perturbed system (N-body interface).
44 * ENERGY Total energy.
45 * EPS Tolerance for DIFSY1 (1.0E-10 is recommended).
46 * ICALL Pericentre indicator (activated if R(IMIN) < EPSR2**0.5).
47 * ICOLL Collision indicator (activated in routine DERQP).
48 * IDER Indicator for setting time derivativative on first call).
49 * IPERT Perturbation indicator (=0: once per call; =1: every call).
50 * I1-I4 Indices of final configuration (I1 & I2 is closest binary).
51 * KSCH Slow-down factor (= 1 when inactive).
52 * KSLOW Slow-down indicator (= .true. when active).
53 * KZ27 Tidal dissipation option.
54 * KZ30 Diagnostic option (copied in CHINIT; full output if > 1).
55 * M Particle mass (CM(7) = sum M(I), I = 1,N).
56 * NAMEC Particle identity (initialized to global name).
57 * NFN Number of function calls.
58 * NPERT Number of perturbers (for diagnostic output).
59 * NREG Number of regularization switches.
60 * NSTEP1 Number of DIFSY1 calls.
61 * P Regularized momenta.
62 * Q Regularized coordinates.
63 * RINV Inverse chain distances.
64 * RCOLL Minimum two-body separation.
65 * RGRAV Gravitational radius ((sum M(I)*M(J))/ABS(ENERGY)).
66 * RMAXS Maximum size of unperturbed configuration.
67 * RMAXC Maximum of unperturbed & initial size (+ 20 per cent).
68 * RSUM Sum of all chain distances.
69 * STEPS Current integration interval (set by routine INTGRT).
70 * TCR Local crossing time ((sum M(I))**2.5/ABS(2*ENERGY)**1.5).
71 * TIMEC Local physical time in scaled units (= CHTIME).
72 * TMAX Maximum integration time (based on c.m. step).
73 * TS Global physical time at each return to N-body integration.
74 * X Particle coordinates (X(1), X(2), X(3) is X, Y, Z).
75 * V Velocity components.
76 * -------------------------------------------------------------------
77 *
78 *
79 * Save termination indicator and check for restart.
80  iterm = isub
81  IF (isub.GT.0) THEN
82 * Choose small step for termination (routine PERMIT & INTGRT).
83  IF (steps(isub).EQ.0.0d0) THEN
84  step = 1.0d-06*step
85  go to 20
86  END IF
87 * Update maximum prediction interval at start of every call.
88  CALL tchain(isub,tsmin)
89  steps(isub) = tsmin
90 * Synchronize next time interval with subsystem step.
91  tmax = timec + steps(isub)
92  go to 20
93  END IF
94 *
95 * Copy initial conditions from N-body COMMON and prepare chain.
96  timec = 0.0d0
97  CALL chinit(isub)
98 *
99 * Initialize diagnostic & COMMON variables.
100  rcoll = 100.0
101  qperi = 100.0
102  zmi1 = 100.0
103  zmi2 = 100.0
104  zmi3 = 100.0
105  ecoll1 = 0.0d0
106  dsc = 1.0
107  DO 2 j = 1,10
108  DO 1 k = 1,10
109  rij(j,k) = 100.0
110  1 CONTINUE
111  2 CONTINUE
112  icall = 0
113  icoll = 0
114  isync = 0
115  ncirc = 0
116  nstepx = 0
117  ider = 0
118  ipert = 1
119  iter = 0
120  imcirc = 0
121  ifail = 0
122  next = 0
123  ndiss1 = 0
124  jc = 0
125  nstep1 = 0
126  nreg = 0
127  nfn = 0
128  ncall = 0
129  kcase = 0
130  j10 = 10
131  iprev = 0
132  icase = .false.
133  names(nmx,5) = 0
134 *
135 * Specify the tolerance for DIFSY1.
136 * EPS = 1.0E-10
137  eps = 1.0e-12
138 *
139 * Initialize subsystem time and set dummy variable for DIFSY.
140  chtime = 0.0d0
141  stime = 0.0d0
142 *
143 * Initialize chain regularization (new case or modified chain).
144  10 neq = 8*n
145  IF (kcase.GT.0) THEN
146  tmax = timec + steps(isub)
147  chtime = timec
148  next = nstep1 + 1
149  zmi1 = 100.0
150  zmi2 = 100.0
151  zmi3 = 100.0
152  rcoll = 100.0
153  END IF
154  y(neq) = chtime
155 *
156 * Define indices for DIFSY1 calls to derivative routine.
157  nc = n - 1
158  lx = 4*nc + 1
159  le = 4*nc + 4
160  lp = 4*nc + 5
161  lv = 8*nc + 5
162  lt = 8*nc + 8
163  nhalf2 = (neq/2)*2
164 *
165 * Ensure routine XTPERT uses correct subsystem index (ISYS(5) is free).
166  isys(5) = isub
167 *
168 * Evaluate constants of the motion and save TPR.
169  CALL const(x,v,m,n,ener0,g0,alag)
170  tpr = 1.0/alag
171  tpr0 = tpr + 1.0
172 *
173 * Select chain indices and transform from chain to KS variables.
174  CALL SELECT
175  CALL transq
176 *
177 * Evaluate inverse distances (for routine INVERT).
178  DO 11 i = 1,n-1
179  ks1 = 4*(i - 1) + 1
180  rk = q(ks1)**2 + q(ks1+1)**2 + q(ks1+2)**2 + q(ks1+3)**2
181  rinv(i) = 1.0/rk
182  11 CONTINUE
183 *
184 * Define total energy (including any c.m. motion & external potential).
185  energy = ener0 - 0.5d0*mass*(cmv(1)**2 + cmv(2)**2 + cmv(3)**2)
186 *
187 * Copy whole input array (for DIFSY call).
188  CALL ycopy(y)
189 *
190 * Find sum of mass products and set current gravitational radius.
191  sum = 0.0d0
192  DO 14 l = 1,n-1
193  DO 12 k = l+1,n
194  sum = sum + mij(l,k)
195  12 CONTINUE
196  14 CONTINUE
197  rgrav = sum/abs(energy)
198 *
199 * Impose initial size as upper limit in case of weakly bound systems.
200  rgrav = min(rgrav,0.5*rsum)
201  rslow = 0.5*rgrav
202 * Check slow-down option for chain (KZ26 copied in CHINIT).
203  IF (kz26.LT.2) rslow = 0.0
204 *
205 * Set current crossing time and specify initial step using T' = 1/L.
206  tcr = mass**2.5/abs(2.0d0*energy)**1.5
207 *
208 * Determine the smallest two-body time-scale from parabolic orbit.
209  CALL r2sort(ij,r2)
210  i1 = ij(1)
211  i2 = ij(2)
212  rm = sqrt(r2(i1,i2))
213  vp2 = 2.0*(m(i1) + m(i2))/rm
214  tp = rm/sqrt(vp2)
215  tstar = min(tp,tcr)
216  stepin = eps**0.2*tstar*alag
217  small = 1.0d-04*stepin
218 * Save step in case of backwards last iteration or collision.
219  saveit = stepin
220 * RUNAV = 0.0001*RGRAV
221 *
222 * (Re-)Initialize slow-down variables and minimum stellar radius.
223  rstar = 0.0
224  do i = 1,n-1
225  ksch(i) = 1.0d0
226  rstar = max(SIZE(i),rstar)
227  end do
228  rstar = max(SIZE(n),rstar)
229 * Set pericentre search distance (twice 4*r_max; discrete intervals).
230  rstar = 8.0*rstar ! no longer used.
231  kslow = .false.
232  kcoll = .false.
233  ejump = 0.0d0
234  stopb = .false.
235 * Specify slow-down criterion with conservative value for N > 4.
236  gcrit = 1.0d-05
237  IF (n.EQ.4) gcrit = 5.0d-06
238  IF (n.GT.4) gcrit = 1.0d-06
239 *
240  IF (kz30.GT.1) THEN
241  WRITE (6,15) n, npert, energy, rsum, rgrav, tcr, rmaxs(isub),
242  & (namec(i),i=1,n)
243  15 FORMAT (' NEW CHAIN N NP E RSUM RGRAV TCR RMAXS NAM ',
244  & 2i4,f10.5,1p,4e9.1,0p,6i6)
245  END IF
246  IF (rsum.LT.1.0e-05) eps = 1.0d-12
247 *
248 * Evaluate binary energies (initially or after membership change).
249  CALL recoil(0)
250 *
251 * Assign integration step (initial case or modified membership).
252  IF (kcase.EQ.0) THEN
253  step = stepin
254  ELSE
255  ipert = 1
256  kcase = 0
257  step = 0.01*stepin
258  CALL tchain(isub,tsmin)
259 * Copy current value in case STEPS = 0 which leads to termination.
260  steps(isub) = tsmin
261 * Ensure termination if only two particles left or positive energy.
262  IF (n.EQ.2.OR.energy.GT.0.0) go to 70
263  END IF
264 *
265 * Perform regularized integration until termination or modification.
266  20 rmaxc = 3.0*rgrav
267 * Replace perturbed boundary radius with MIN(3*RGRAV,2*RSUM).
268  rmaxc = min(rmaxc,2.0*rsum)
269  21 step0 = step
270  time0 = chtime
271 *
272 * Modify DIFSY integration order according to relative perturbation.
273 * JDIF = 10 - (1.0D+04*GPERT)**0.4
274 * IF (JDIF.LT.J10) THEN
275 * J10 = MAX(4,JDIF)
276 * ELSE IF (JDIF.GT.J10) THEN
277 * J10 = MIN(10,J10+1)
278 * END IF
279 *
280 * Obtain maximum step from two-body iteration based on c-functions.
281  IF (icoll.EQ.0.AND.icall.EQ.0.AND.steps(isub).GT.0.0d0) THEN
282 * Set slightly large interval to avoid small steps and check STEPS.
283  dt = 1.01*tmax - chtime
284  IF (dt.LT.0.0) dt = abs(dt)
285  dt = min(1.01*steps(isub),dt)
286 * Perform iteration of DTAU based on dominant two-body motion.
287  CALL invert(dt,dtau)
288 * Increase step by 10% during expansion and check convergence value.
289  IF (tpr.GT.tpr0) dtau = 1.1*dtau
290 * Avoid choosing a negative step from previous backward integration.
291  IF (step.GT.0.0) THEN
292  step = min(dtau,step)
293  ELSE
294  step = min(dtau,saveit)
295  END IF
296  END IF
297 *
298 * Copy restart step after tidal dissipation or collision/failure.
299  IF (icase.AND.iprev.GT.0) THEN
300  step = saveit
301  step0 = step
302  iprev = 0
303  icase = .false.
304  END IF
305 *
306 * Check time-step rejuvenation during standard integration.
307  IF (mod(nstep1,1000).EQ.0.AND.icall.EQ.0.AND.icoll.EQ.0) THEN
308  CALL const(x,v,m,n,ener0,g0,alag)
309  small = 1.0d-04*eps**0.2*tstar*alag
310 * Increase by modest factor to avoid danger of zero step.
311  IF (step.LT.small) step = small
312 * WRITE (6,18) NSTEP1, SMALL, STEP, (1.0/RINV(K),K=1,N-1)
313 * 18 FORMAT (' STEP INCREASE # SM S R ',I8,1P,2E9.1,2X,5E10.2)
314  ELSE
315 * Suppress ICALL after 10 unsuccessful attempts.
316  IF (ifail.GT.10) icall = 0
317  END IF
318 *
319 * Advance the solution one step.
320  CALL difsy1(neq,eps,step,stime,y)
321 *
322 * Copy current physical time and save COMMON variables.
323  chtime = y(neq)
324  timec = chtime
325  CALL ysave(y)
326  CALL transx
327  nstep1 = nstep1 + 1
328  tpr0 = tpr
329 *
330 * Save new step during standard integration for subsequent restart.
331  22 IF (icall.EQ.0.AND.icoll.EQ.0) THEN
332  saveit = step
333  icase = .true.
334 * Reset collision indicator and activate IPREV on failed iteration.
335  ELSE IF (icoll.LT.0.AND..NOT.kcoll) THEN
336  icoll = 0
337  iprev = 1
338  icase = .true.
339  kcoll = .false.
340  next = nstep1 + 1
341  step = saveit
342 * Restore the original configuration and copy to input array.
343  IF (imcirc.GT.0) isync = 1
344  tpr = tkpr
345  DO 24 i = 1,n-1
346  ks = 4*(i - 1)
347  DO 23 j = 1,4
348  q(ks+j) = qk(ks+j)
349  p(ks+j) = pk(ks+j)
350  23 CONTINUE
351  24 CONTINUE
352  CALL ycopy(y)
353  go to 21
354  END IF
355 *
356 * Check slow-down procedure for small minimum distance.
357  IF (kz26.GE.2.AND.(.not.stopb.or.kslow)) then
358  rm = 0.0
359  DO 25 i = 1,n-1
360  IF (rm.LT.rinv(i)) THEN
361  rm = rinv(i)
362  im = i
363  END IF
364  25 CONTINUE
365  rm = 1.0/rm
366  IF (rm.LT.rslow.AND.icoll.EQ.0) THEN
367  CALL slow(y)
368  END IF
369  END IF
370 *
371 * Include possible emergency stop of slow-down (DERQP; inactive).
372  if (stopb) then
373  CALL ycopy(y)
374  stopb = .false.
375  it = 0
376  28 CALL difsy1(neq,eps,stepl,stime,y)
377  if (.not.stopb) then
378  CALL ysave(y)
379  else if (it.LT.5) then
380  it = it + 1
381  go to 28
382  end if
383  CALL slow(y)
384  go to 20
385  end if
386 *
387  IF (step.GT.10.0*step0.AND.icoll.EQ.0) THEN
388  step = 10.0*abs(step0)
389  ELSE IF (step.EQ.0.0d0) THEN
390  WRITE (6,*) ' Stepsize = 0!', char(7)
391  stop
392  END IF
393 *
394  IF (kz30.GT.2) THEN
395  WRITE (6,30) step, tmax-chtime, gpert, (1.0/rinv(k),k=1,n-1)
396  30 FORMAT (' CHAIN: STEP TM-CHT G R ',1p,8e9.1)
397  END IF
398 *
399 * Determine two-body distances for stability test and collision search.
400  IF (chtime.GT.time0.AND.jc.EQ.0) THEN
401  rm = 0.0
402  rc2 = 0.0
403  DO 35 k = 1,n-1
404 * Find minimum separation for stability test and save chain index.
405  IF (rinv(k).GT.rm) THEN
406  rm = rinv(k)
407  ix = k
408  END IF
409  rc2 = rc2 + rinv(k)**2
410  35 CONTINUE
411 * Update sum of 1/R^2 during forward integration (R^2 before 12/99).
412  rm = 1.0/rm
413  zmi = rc2
414  zmi1 = zmi2
415  zmi2 = zmi3
416  zmi3 = zmi
417 * RUNAV = 0.9*RUNAV + 0.1*RM
418 * Set search distance for closest separation.
419  i1 = iname(ix)
420  i2 = iname(ix+1)
421 * Replace sum of radii by maximum value (08/08).
422 * SX = SIZE(I1) + SIZE(I2)
423  sx = max(SIZE(i1),SIZE(i2))
424 * Turn off circular orbit indicators for overlapping stars.
425  IF (imcirc.GT.0) THEN
426 * Prevent repeated switching (old test RM < 2*MAX(SIZE); 07/08).
427  IF (rm.LT.sx) THEN
428  sx = 0.5*rm
429  ncirc = 0
430  isync = 0
431  imcirc = 0
432  ELSE
433  sx = 0.0
434  END IF
435  END IF
436 * See whether another star is a more likely collider (compare SIZE/R).
437  ry = rm
438  DO 36 k = 1,n-1
439  IF (k.NE.ix) THEN
440  j1 = iname(k)
441  j2 = iname(k+1)
442  sy = max(SIZE(j1),SIZE(j2))
443  IF (sy*rinv(k).GT.sx/ry) THEN
444  sx = sy
445  ry = 1.0/rinv(k)
446  END IF
447  END IF
448  36 CONTINUE
449 * Include factor of 3 for criterion QPMIN < 4*MAX(SIZE(K1),SIZE(K2)).
450  sx = 3.0*sx
451  gx = 0.0
452 * Compare sum of 1/R**3 neighbouring perturbations with dominant term.
453  DO 38 k = 1,n-1
454  IF (iabs(k-ix).EQ.1) THEN
455  gx = gx + rinv(k)**3
456  END IF
457  38 CONTINUE
458 * Adopt safety ratio 0.01 for initializing close interaction search.
459  pmax = 0.01*rinv(ix)**3
460  ELSE
461  gx = 0.0
462  pmax = 1.0
463  END IF
464 *
465 * Switch on search indicator during iteration or just after pericentre.
466  IF (icoll.LT.0) icall = 1
467  IF (ry.LT.sx.AND.nstep1.GT.next) THEN
468  IF (zmi3.LT.zmi2.AND.zmi2.GT.zmi1) THEN
469 * Skip pericentre determination on large perturbation (avoids looping).
470  IF (isync.EQ.0.AND.ncirc.LT.25.AND.gx.LT.pmax) THEN
471  icall = 1
472  IF (kz26.GE.2.AND.kslow) THEN
473  gsave = gcrit
474  gcrit = 0.0
475  ncirc = ncirc + 1
476  CALL slow(y)
477  gcrit = gsave
478  END IF
479  END IF
480  END IF
481  END IF
482 *
483 * Restore indicators for further attempts after 100 steps.
484  IF (isync.GT.0.AND.nstep1.GT.nstepx + 100) THEN
485  ncirc = 0
486  isync = 0
487  imcirc = 0
488  nstepx = nstep1
489  END IF
490 *
491 * Check for collision, capture or circularization during last step.
492  IF (icoll.LE.0) go to 50
493  iprev = icoll
494 *
495 * Restore the minimum configuration from DERQP.
496  tpr = tkpr
497  DO 45 i = 1,n-1
498  ks = 4*(i - 1)
499  DO 40 j = 1,4
500  q(ks+j) = qk(ks+j)
501  p(ks+j) = pk(ks+j)
502  40 CONTINUE
503  45 CONTINUE
504 *
505 * Delay next search by one step to avoid the same pericentre.
506  next = nstep1 + 1
507  icall = 0
508 *
509 * Transform to physical variables.
510  CALL transx
511 *
512 * Distinguish between tidal energy loss and collision (ICOLL holds IM).
513  im = icoll
514  icoll = 0
515 *
516 * Check for capture, circularization or collision.
517  IF (qperi.LT.4.0*max(SIZE(k1),SIZE(k2))) THEN
518 *
519 * Distinguish between tidal energy loss and collision.
520  j1 = k1
521  j2 = k2
522  IF (SIZE(k2).GT.SIZE(k1)) THEN
523  j1 = k2
524  j2 = k1
525  END IF
526 *
527 * Avoid tidal dissipation for circularized orbit.
528  icirc = 0
529  IF (kz27.GT.0.AND.ebs.LT.0.0) THEN
530  semi = -0.5*m(k1)*m(k2)/ebs
531  ecc = 1.0 - qperi/semi
532  IF (ecc.GT.0.0022) icirc = 1
533  END IF
534 *
535 * Adopt collision criterion of Kochanek (Ap.J. 385, 684, 1992)
536  IF (kz27.LE.2) THEN
537  fac = 0.5*(m(j1) + m(j2))/m(j1)
538  rcr = 1.7*fac**0.3333*SIZE(j1)
539  ELSE
540  clight = 3.0d+05/vstar1
541  rcr = 6.0*(m(j1) + m(j2))/clight**2
542  END IF
543 *
544  IF (qperi.LT.rcr) THEN
545 * Obtain global coordinates & velocities (ITERM < 0 denotes collision).
546  ifail = 0
547  iterm = -1
548  isub = -isub
549  CALL chterm(isub)
550 *
551 * Combine internal energy and external c.m. kinetic energy.
552  hc = energy + 0.5d0*mass*(cm(4)**2 + cm(5)**2 + cm(6)**2)
553 *
554 * Determine dominant two-body energy from non-singular expressions.
555  CALL erel(im,ebs,semi)
556 *
557 * Form composite body and terminate with new KS or continue chain.
558  CALL cmbody(hc,n)
559 *
560 * Include chain restart (denoted by N < 0) for CE with > 3 members.
561  IF (n.LT.0) THEN
562  n = -n
563  CALL chinit(isub)
564  go to 10
565  END IF
566 * Terminate for one single body after mass-less collision/coalescence.
567  IF (n.EQ.0) THEN
568  iterm = -1
569  go to 100
570  END IF
571 * Re-initialize with reduced membership NCH.
572  n = n - 1
573  CALL chinit(isub)
574 * Activate indicator for new chain (terminates at once with N = 2).
575  kcase = 1
576  go to 10
577  ELSE IF (icirc.GT.0.AND.kz27.GT.0.AND.isync.EQ.0) THEN
578 * Modify regularized variables due to tidal dissipation.
579  en0 = energy
580  CALL qpmod(im,iterm)
581 *
582 * Check reduction (ITERM = -2) or termination (ITERM = -1; N <= 4).
583  IF (iterm.EQ.-2) THEN
584  IF (n.EQ.2) go to 70
585  go to 10
586  END IF
587 * Re-determine the gravitational radius (cf. escape delay criterion).
588  rgrav = sum/abs(energy)
589  CALL ycopy(y)
590  IF (iterm.EQ.-1.AND.n.LE.4) THEN
591  chtime = y(neq)
592  CALL ysave(y)
593  CALL transx
594  timec = chtime
595  ech = energy
596  go to 70
597  END IF
598 * Suppress ICALL before integration after ten zero energy changes.
599  IF (abs(energy-en0).EQ.0.0d0) THEN
600  ifail = ifail + 1
601  END IF
602  go to 21
603  ELSE
604 * Allow for case of too long circularization time.
605  isync = 1
606  icoll = -1
607  kcoll = .false.
608  go to 22
609  END IF
610  ELSE
611 * Include case of large pericentre (first or second distance).
612  IF (imcirc.GT.0) THEN
613  isync = 1
614  icoll = -1
615  go to 22
616  END IF
617  CALL ycopy(y)
618  go to 21
619  END IF
620 *
621 * Check switching condition (Note: RINV now updated after switch).
622  50 isw = 0
623  CALL swcond(check)
624  IF (check) THEN
625  DO 52 i = 1,n
626  iold(i) = iname(i)
627  52 CONTINUE
628  CALL switch(y)
629  isw = 1
630 * See whether the chain name list is unchanged (even if CHECK =.TRUE.).
631  DO 54 i = 1,n
632  IF (iold(i).NE.iname(i)) isw = isw + 1
633  54 CONTINUE
634  IF (isw.GT.1) nreg = nreg + 1
635 * Update slow-down if relevant (avoids perturbation jump).
636  IF (kslow.AND.isw.GT.1) THEN
637  CALL slow(y)
638  END IF
639  END IF
640 *
641 * Check termination or strong perturbation (T > TMAX or GPERT > 0.01).
642  IF (chtime.GT.tmax.OR.gpert.GT.0.01) THEN
643  CALL ysave(y)
644  CALL transx
645  ech = energy
646  timec = chtime
647  IF (kz30.GT.2) THEN
648  WRITE (6,55) nstep1, t0s(isub)+timec, tmax-timec,
649  & (1.0/rinv(k),k=1,n-1)
650  55 FORMAT (' CHAIN: # T DTR R ',i5,f10.4,1p,6e9.1)
651  END IF
652 * Avoid checking after switch (just in case).
653  IF (isw.LE.1) THEN
654 * Update RGRAV in case of compact initial size.
655  rgrav = min(sum/abs(energy),0.5*rsum)
656  CALL chmod(isub,kcase)
657  IF (kcase.GT.0) THEN
658  CALL recoil(1)
659  go to 10
660  END IF
661  IF (kcase.LT.0) go to 60
662  END IF
663  IF (isw.EQ.0.AND.n.EQ.3) THEN
664  IF (rsum.GT.4.0*rm) THEN
665  CALL chstab(iterm)
666  IF (iterm.LT.0) go to 70
667  END IF
668  END IF
669  IF (chtime.LT.tmax.AND.steps(isub).GT.0.0d0) go to 21
670  go to 60
671  ELSE
672 * Exit chain integration if time limit exceeded or STEPS = 0.
673  IF (chtime.GT.tmax) go to 60
674 * Terminate on large step number and tiny STEP > 0.
675  IF (nstep1.GT.100000.AND.step.LT.small.AND.step.GT.0.0) THEN
676  WRITE (6,58) nstep1, n, step, (1.0/rinv(k),k=1,n-1)
677  58 FORMAT (' ENFORCED CHAIN # N STEP R',i8,i4,1p,9e9.1)
678  go to 70
679  END IF
680  IF (steps(isub).GT.0.0d0) go to 21
681  END IF
682 *
683 * Check hierarchical stability condition for triple or quad.
684  60 IF (n.EQ.3) THEN
685  IF (rsum.GT.4.0*rm) THEN
686  CALL chstab(iterm)
687  IF (iterm.LT.0) go to 70
688  END IF
689  ELSE IF (n.EQ.4) THEN
690  IF (rm.LT.0.1*rsum) THEN
691 * Find largest separation to distinguish triple or quad case.
692  rx = 1.0d+10
693  DO 65 k = 1,n-1
694  rx = min(rx,rinv(k))
695  65 CONTINUE
696  rx = 1.0/rx
697 * Check case of two binaries or degenerate triple and small binary.
698  IF (rx.GT.0.7*rsum) THEN
699  CALL cstab4(iterm)
700  IF (iterm.EQ.0.AND.ix.NE.2) THEN
701  CALL cstab2(iterm)
702  END IF
703 * Skip small middle distance (done by CSTAB3 called from CSTAB4).
704  ELSE IF (rm.LT.0.01*rsum.AND.ix.NE.2) THEN
705  CALL cstab2(iterm)
706  END IF
707  IF (iterm.LT.0) go to 70
708  END IF
709 * Reduce five/six-body system to triple if biggest binary < 0.04*RSUM.
710  ELSE IF (n.GE.5) THEN
711  CALL cstab5(iterm)
712  IF (iterm.LT.0) go to 70
713  END IF
714 *
715 * See whether temporary or actual termination (continue if N > 5).
716  IF (n.GT.5.OR.(kcase.EQ.0.AND.steps(isub).GT.0.0d0)) go to 100
717  IF (kcase.LT.0) go to 70
718  IF (timec.GT.tmax.AND.rsum.LT.rmaxc) THEN
719  IF (steps(isub).GT.0.0d0.OR.n.GT.4) go to 100
720  END IF
721 *
722 * Check for dominant binary.
723  70 CALL recoil(2)
724 *
725 * Set zero step to define termination (just in case).
726  steps(isub) = 0.0d0
727 *
728 * Copy Q & P for TRANSK via CHTERM & EREL (KZ26 > 2: rectification).
729  IF (kz26.GT.2) THEN
730  DO 75 i = 1,n-1
731  ks = 4*(i - 1)
732  DO 72 j = 1,4
733  qk(ks+j) = q(ks+j)
734  pk(ks+j) = p(ks+j)
735  72 CONTINUE
736  rik = q(ks+1)**2 + q(ks+2)**2 + q(ks+3)**2 + q(ks+4)**2
737  rinv(i) = 1.0/rik
738  75 CONTINUE
739  END IF
740 *
741 * Transform to global variables and begin new KS (I1 & I2), I3 & I4.
742  CALL chterm(isub)
743 *
744 * Activate termination index for routine INTGRT.
745  iterm = -1
746 *
747  IF (kz30.GT.1.AND.qperi.LT.1.0) THEN
748  WRITE (6,80) rij(1,2), rij(1,3), rij(2,3), rcoll, qperi
749  80 FORMAT (' END CHAIN RIJ RCOLL QPERI ',1p,5e10.1)
750  END IF
751 *
752 * Update current time unless termination and set subsystem index.
753  100 IF (iterm.GE.0) ts(isub) = t0s(isub) + timec
754  isub = iterm
755 *
756  RETURN
757 *
758  END