Nbody6
 All Files Functions Variables
chterm.f
Go to the documentation of this file.
1  SUBROUTINE chterm(ISUB)
2 *
3 *
4 * Termination of chain system.
5 * ----------------------------
6 *
7  include 'common6.h'
8  parameter(nmx=10,nmx2=2*nmx,nmx3=3*nmx,nmx4=4*nmx,
9  & nmx8=8*nmx,nmxm=nmx*(nmx-1)/2)
10  REAL*8 m,mass,mc,mij,mkk,r2(nmx,nmx),ang(3)
11  INTEGER ij(nmx)
12  common/chain1/ xch(nmx3),vch(nmx3),m(nmx),
13  & zz(nmx3),wc(nmx3),mc(nmx),
14  & xi(nmx3),pi(nmx3),mass,rinv(nmxm),rsum,mkk(nmx),
15  & mij(nmx,nmx),tkk(nmx),tk1(nmx),iname(nmx),nn
16  common/chainc/ xc(3,ncmax),uc(3,ncmax),bodyc(ncmax),ich,
17  & listc(lmax)
18  common/chreg/ timec,tmax,rmaxc,cm(10),namec(6),nstep1,kz27,kz30
19  common/clump/ bodys(ncmax,5),t0s(5),ts(5),steps(5),rmaxs(5),
20  & names(ncmax,5),isys(5)
21  common/ccoll2/ qk(nmx4),pk(nmx4),rik(nmx,nmx),SIZE(nmx),vstar1,
22  & ecoll1,rcoll,qperi,istar(nmx),icoll,isync,ndiss1
23  common/incond/ x4(3,nmx),xdot4(3,nmx)
24  common/echain/ ech
25  common/ksave/ k1,k2
26 *
27 *
28 * Decide between standard termination or collision (ISUB > 0 or < 0).
29  IF (isub.NE.0) THEN
30  iterm = isub
31  isub = iabs(isub)
32  END IF
33 *
34 * Prepare KS regularization(s) and direct integration of other bodies.
35  CALL r2sort(ij,r2)
36 * Save dominant pair (note change from MIN(R2) to M/R2; SJA 08/2009).
37  i1 = ij(1)
38  i2 = ij(2)
39  i3 = ij(3)
40  i5 = 0
41  i6 = 0
42  IF (nch.EQ.2) i3 = i2
43  IF (nch.LE.3) THEN
44  i4 = i3
45  i5 = i3
46  r2(i2,i4) = r2(i2,i3)
47  r2(i3,i4) = 0.0
48  ELSE IF (nch.EQ.4) THEN
49  i4 = ij(4)
50  ELSE IF (nch.GT.4) THEN
51 * Determine indices for second closest pair.
52  rx1 = 1.0
53  rx0 = r2(i1,i2)
54  DO 2 j1 = 1,nch
55 * Avoid choosing close pair I1-I2.
56  IF (j1.EQ.i1.OR.j1.EQ.i2) go to 2
57  DO 1 j2 = j1+1,nch
58  IF (j2.EQ.i1.OR.j2.EQ.i2) go to 1
59  IF (r2(j1,j2).LT.rx1.AND.r2(j1,j2).GT.rx0) THEN
60  rx1 = r2(j1,j2)
61  i3 = j1
62  i4 = j2
63  END IF
64  1 CONTINUE
65  2 CONTINUE
66 * Identify remaining single particle(s) by exclusion.
67  DO 3 i = 1,nch
68  IF (i.EQ.i1.OR.i.EQ.i2.OR.i.EQ.i3.OR.i.EQ.i4) go to 3
69  IF (i5.EQ.0) THEN
70  i5 = i
71  ELSE
72  i6 = i
73  END IF
74  3 CONTINUE
75  END IF
76  IF (i5.EQ.0) i5 = i4
77  IF (i6.EQ.0) i6 = i4
78 *
79  IF (kz(30).GT.1) THEN
80  WRITE (6,4) sqrt(r2(i1,i2)), sqrt(r2(i1,i3)),sqrt(r2(i2,i3)),
81  & sqrt(r2(i2,i4)), sqrt(r2(i3,i4))
82  4 FORMAT (' CHTERM: RIJ (1-2 1-3 2-3 2-4 3-4) ',1p,5e9.1)
83  END IF
84 *
85  jlist(6) = namec(i1)
86  jlist(7) = namec(i2)
87  jlist(8) = namec(i3)
88  jlist(9) = namec(i4)
89  jlist(10) = namec(i5)
90  IF (nch.EQ.6) jlist(11) = namec(i6)
91 *
92 * Specify chain phase indicator and restore original name to c.m. body.
93  iphase = 8
94  name(ich) = name0
95 *
96 * Identify current global indices by searching all particles.
97  icm = 0
98  DO 10 j = ifirst,ntot
99  DO 5 l = 1,nch
100  IF (name(j).EQ.jlist(l+5)) THEN
101  jlist(l) = j
102  IF (body(j).GT.0.0d0) icm = j
103  END IF
104  5 CONTINUE
105  10 CONTINUE
106 *
107 * Modify identification list for special cases NCH = 2 & NCH = 3.
108  IF (nch.EQ.2) THEN
109  i3 = i1
110  i4 = i2
111  jlist(3) = jlist(1)
112  jlist(4) = jlist(2)
113  ELSE IF (nch.EQ.3) THEN
114  jlist(4) = jlist(3)
115  END IF
116 *
117 * Ensure ICOMP < JCOMP for KS regularization.
118  icomp = min(jlist(1),jlist(2))
119  jcomp = max(jlist(1),jlist(2))
120 *
121 * Copy final coordinates & velocities to standard variables.
122  lk = 0
123  rx2 = 0.0
124  DO 20 l = 1,nch
125  ri2 = 0.0
126  vi2 = 0.0
127  DO 15 k = 1,3
128  lk = lk + 1
129  x4(k,l) = xch(lk)
130  xdot4(k,l) = vch(lk)
131  ri2 = ri2 + x4(k,l)**2
132  vi2 = vi2 + xdot4(k,l)**2
133  15 CONTINUE
134 * Save velocity & NAME of distant escaper for c.m. diagnostics.
135  IF (ri2.GT.rx2) THEN
136  rx2 = ri2
137  vx2 = vi2
138  namex = namec(l)
139  END IF
140  20 CONTINUE
141 *
142 * Quantize the elapsed interval since last step (note: none at TBLOCK).
143  time2 = t0s(isub) + timec - tprev
144  dt8 = (tblock - tprev)/8.0d0
145 *
146 * Adopt the nearest truncated step (at most 8 subdivisions).
147  dt2 = time2
148  IF (time2.GT.0.0d0) THEN
149  CALL stepk(dt2,dtn2)
150  dtn = nint(dtn2/dt8)*dt8
151  ELSE
152 * Choose negative step if pericentre time < TPREV (cf. iteration).
153  dt2 = -dt2
154  CALL stepk(dt2,dtn2)
155  dtn = -nint(dtn2/dt8)*dt8
156  END IF
157 *
158 * Update time for new polynomial initializations (also for CMBODY).
159  time = tprev + dtn
160 * Adopt block time in rare case of DT2 reduced to SMAX.
161  IF (dtn.LE.smax) time = tblock
162 *
163 * Set TIMEC for collision to preserve TIME2 after T0S = TIME in SUBSYS.
164  IF (iterm.LT.0.0) timec = t0s(isub) + timec - time
165 *
166 * Restrict TIME to current block-step and save for use by KSPERI.
167  time = min(tblock,time)
168  time0 = time
169 *
170 * Predict current X & XDOT for c.m. and neighbours to order F3DOT.
171  CALL xvpred(icm,-1)
172  nnb1 = list(1,ich) + 1
173  DO 25 l = 2,nnb1
174 * Note possibility T0(J) = TIME in routine REDUCE would skip on -2.
175  j = list(l,ich)
176  CALL xvpred(j,-1)
177  25 CONTINUE
178 *
179 * Copy c.m. coordinates & velocities.
180  DO 30 k = 1,3
181  cm(k) = x(k,icm)
182  cm(k+3) = xdot(k,icm)
183  30 CONTINUE
184 *
185 * Set configuration pointers for KS candidates & distant bodies.
186  jlist(7) = i1
187  jlist(8) = i2
188  jlist(9) = i3
189  jlist(10) = i4
190  jlist(11) = i5
191  jlist(12) = i6
192 *
193 * Place new coordinates in the original locations.
194  DO 40 l = 1,nch
195  j = jlist(l)
196 * Compare global name & subsystem name to restore the mass & T0.
197  DO 32 k = 1,nch
198  IF (name(j).EQ.namec(k)) THEN
199  body(j) = bodyc(k)
200  t0(j) = time
201  END IF
202  32 CONTINUE
203 * Transform to global coordinates & velocities using c.m. values.
204  ll = jlist(l+6)
205  vj2 = 0.0
206  DO 35 k = 1,3
207  x(k,j) = x4(k,ll) + cm(k)
208  xdot(k,j) = xdot4(k,ll) + cm(k+3)
209  x0(k,j) = x(k,j)
210  x0dot(k,j) = xdot(k,j)
211  vj2 = vj2 + xdot(k,j)**2
212  35 CONTINUE
213  40 CONTINUE
214 *
215 * Print diagnostics for high-velocity chain escaper after COLL or COAL.
216  IF (ech.GT.0.5*bodym*16.0*eclose) THEN
217  WRITE (6,41) namex, nch, ech, sqrt(vx2)*vstar
218  41 FORMAT (' CHAIN DISRUPT NMX NCH ECH VX ',i7,i4,f8.4,f6.1)
219  END IF
220 *
221 * Set one ghost index for skipping neighbour list search.
222  jg = icomp
223  IF (jg.EQ.icm) jg = jcomp
224 *
225 * Search all neighbour lists for splitting chain c.m. into components.
226  nnb1 = 0
227  DO 50 j = ifirst,ntot
228  nnb2 = list(1,j) + 1
229  DO 45 l = 2,nnb2
230  IF (list(l,j).GT.icm) go to 50
231 * Form list of all ICM without one ghost and predict X & XDOT to FDOT.
232  IF (list(l,j).EQ.icm) THEN
233  DO 42 k = 2,nnb2
234  IF (list(k,j).EQ.jg) go to 50
235  42 CONTINUE
236 * Skip subsystem members.
237  DO 44 k = 1,nch
238  IF (j.EQ.jlist(k)) go to 50
239  44 CONTINUE
240 * Copy all relevant indices to JPERT (note limit of 5*LMAX).
241  IF (nnb1.LT.5*lmax) THEN
242  nnb1 = nnb1 + 1
243  jpert(nnb1) = j
244  CALL xvpred(j,0)
245  END IF
246  END IF
247  45 CONTINUE
248  50 CONTINUE
249 *
250 * Check for stellar collision (only needs coordinates & velocities).
251  IF (iterm.LT.0) THEN
252  jlist(1) = icomp
253  jlist(2) = jcomp
254 * See whether re-labelling is required (indices I1 - I4 still local).
255  IF (r2(i1,i4).LT.r2(i1,i3).OR.r2(i3,i4).LT.r2(i1,i3)) THEN
256  IF (r2(i1,i4).LT.r2(i3,i4)) THEN
257 * Switch body #I3 & I4 to give new dominant pair I1 & I3.
258  i = jlist(4)
259  jlist(4) = jlist(3)
260  jlist(3) = i
261  ELSE
262 * Set JLIST(7) < 0 to denote that body #I3 & I4 will be new KS pair.
263  jlist(7) = -1
264  END IF
265  END IF
266 * Define collision distance and indicator for routine CMBODY.
267  dminc = rcoll
268  iphase = 9
269  go to 100
270  END IF
271 *
272 * Update subsystem COMMON variables unless last or only case.
273  IF (isub.LT.nsub) THEN
274  DO 60 l = isub,nsub
275  DO 55 k = 1,6
276  bodys(k,l) = bodys(k,l+1)
277  names(k,l) = names(k,l+1)
278  55 CONTINUE
279  t0s(l) = t0s(l+1)
280  ts(l) = ts(l+1)
281  steps(l) = steps(l+1)
282  rmaxs(l) = rmaxs(l+1)
283  isys(l) = isys(l+1)
284  60 CONTINUE
285  END IF
286 *
287 * Assign new neighbours for dominant KS and any other members.
288  rs0 = rs(icm)
289  CALL nblist(icomp,rs0)
290  DO 65 l = 3,nch
291  j = jlist(l)
292  CALL nblist(j,rs0)
293  65 CONTINUE
294 *
295 * Replace ICM in neighbour lists by all subsystem members.
296  CALL nbrest(icm,nch,nnb1)
297 *
298 * Exclude the dominant interaction for c.m. approximation (large FDOT).
299  IF (min(r2(i1,i3),r2(i2,i4)).GT.cmsep2*r2(i1,i2)) THEN
300  jlist(1) = jlist(3)
301  jlist(2) = jlist(4)
302  nnb = 2
303  IF (nch.EQ.3) jlist(2) = max(ifirst + 4,jlist(1) + 1)
304  ELSE
305  nnb = nch
306  END IF
307 *
308  IF (nch.EQ.2) THEN
309  jlist(1) = max(ifirst + 4,jcomp + 1)
310  nnb = 1
311  END IF
312 *
313 * Set dominant F & FDOT on body #ICOMP & JCOMP for #I3 & I4 in FPOLY2.
314  IF (nch.GE.3) THEN
315  CALL fclose(icomp,nnb)
316  CALL fclose(jcomp,nnb)
317  END IF
318 *
319 * See whether a second binary is present (NCH >= 4 & RB < RMIN).
320  ks2 = 0
321  IF (nch.GE.4.AND.r2(i3,i4).LT.rmin2) THEN
322 * Accept second KS pair if well separated (> 2) from smallest binary.
323  IF (r2(i3,i4).LT.0.25*min(r2(i1,i3),r2(i2,i4))) THEN
324  IF (max(jlist(3),jlist(4)).LE.n) THEN
325  ks2 = 1
326  END IF
327  END IF
328  END IF
329 *
330 * Specify global indices of least dominant bodies (I4 = I3 if NCH = 3).
331  i3 = jlist(3)
332  i4 = jlist(4)
333  IF (nch.GE.5) i5 = jlist(5)
334  IF (nch.EQ.6) i6 = jlist(6)
335 *
336 * Copy name of two first perturbers in case of exchange in KSREG.
337  IF (listc(1).GT.0) THEN
338  namc2 = name(listc(2))
339  IF (listc(1).GT.1) THEN
340  namc3 = name(listc(3))
341  END IF
342  END IF
343 *
344 * Save global indices of #I3 & I4 for initialization of second KS pair.
345  IF (ks2.GT.0) THEN
346  name3 = name(i3)
347  name4 = name(i4)
348  jlist(1) = icomp
349  jlist(2) = jcomp
350  CALL fclose(i3,4)
351  CALL fclose(i4,4)
352  END IF
353 *
354 * Initialize force polynomials & time-steps for body #I3 & #I4.
355  IF (nch.GT.2.AND.ks2.EQ.0) THEN
356  CALL fpoly1(i3,i3,0)
357  IF (nch.GE.4) THEN
358  CALL fpoly1(i4,i4,0)
359  CALL fpoly2(i4,i4,0)
360  IF (nch.GE.5.AND.i5.GT.0) THEN
361  CALL fpoly1(i5,i5,0)
362  CALL fpoly2(i5,i5,0)
363  IF (nch.EQ.6) THEN
364  CALL fpoly1(i6,i6,0)
365  CALL fpoly2(i6,i6,0)
366  END IF
367  END IF
368  END IF
369  CALL fpoly2(i3,i3,0)
370  step3 = step(i3)
371  step4 = step(i4)
372 *
373 * Check re-initialization of dormant super-hard binary (second pair).
374  IF (max(i3,i4).GT.n) THEN
375  i = max(i3,i4)
376  CALL renew(i)
377  END IF
378 * Include initialization of #I5/I6 before KS calls change address.
379  ELSE IF (nch.GE.5.AND.i5.GT.0) THEN
380  CALL fpoly1(i5,i5,0)
381  CALL fpoly2(i5,i5,0)
382  IF (nch.EQ.6) THEN
383  CALL fpoly1(i6,i6,0)
384  CALL fpoly2(i6,i6,0)
385  END IF
386  END IF
387 *
388 * Perform KS regularization of dominant components (ICOMP < JCOMP).
389  IF (jcomp.LE.n) THEN
390  CALL ksreg
391 * Produce diagnostics on ejection velocity of escaper > 2*VSTAR.
392  IF (sqrt(vx2)*vstar.GT.2.0*vstar) THEN
393  j1 = 2*npairs - 1
394  j2 = j1 + 1
395  a1 = -0.5*body(ntot)/h(npairs)
396  pb = days*a1*sqrt(abs(a1)/body(ntot))
397  vi2 = xdot(1,ntot)**2 + xdot(2,ntot)**2 + xdot(3,ntot)**2
398  vcm = sqrt(vi2)*vstar
399  WRITE (6,68) name(j1), name(j2), body(j1)*smu,
400  & body(j2)*smu, vcm, a1, pb
401  68 FORMAT (' CHAIN BINARY NM M1 M2 VCM A PB ',
402  & 2i7,2f6.1,f7.1,1p,2e10.2)
403  END IF
404 * Restore TIME in case modified by routine KSPERI.
405  time = time0
406 * Include optional rectification by standard KS procedure.
407  IF (kz(26).GT.2) THEN
408  rx = 0.0
409 * Determine chain index of closest pair.
410  DO 70 k = 1,nch-1
411  IF (rinv(k).GT.rx) THEN
412  rx = rinv(k)
413  im = k
414  END IF
415  70 CONTINUE
416 * Set dominant pair mass references as IM & IM+1.
417  k1 = iname(im)
418  k2 = iname(im+1)
419 * Obtain regular values of binding energy and semi-major axis.
420  CALL erel(im,eb,semi)
421  hi = -0.5*body(ntot)/semi
422 * Evaluate total energy explicitly.
423  CALL const(xch,vch,m,nch,ener2,ang,gam)
424  err = (ech - ener2)/ech
425  eb1 = eb/ech
426  gi = gamma(npairs)
427 * Rectify KS elements of dominant and weakly perturbed binary.
428  IF (eb1.GT.0.9.AND.gi.LT.gmax) THEN
429  IF (abs(err).GT.1.0d-08) THEN
430  WRITE (6,72) nch, nstep1, eb1, semi, gi, err
431  72 FORMAT (' CHAIN RECTIFY NCH # EB/E A G ERR ',
432  & i4,i6,f6.2,1p,3e10.2)
433  END IF
434 * Copy new value of H and employ standard rectification.
435  h(npairs) = hi
436  CALL ksrect(npairs)
437  END IF
438  END IF
439  ELSE
440 * Initialize components separately and re-activate dormant binary.
441  CALL nblist(jcomp,rs0)
442  CALL fpoly1(jcomp,jcomp,0)
443  CALL renew(jcomp)
444  CALL fpoly1(icomp,icomp,0)
445  CALL fpoly2(icomp,icomp,0)
446  CALL fpoly2(jcomp,jcomp,0)
447  END IF
448 *
449 * Check for high-velocity ejection with outwards hyperbolic motion.
450  rij2 = 0.0
451  vij2 = 0.0
452  rdot = 0.0
453  IF (nch.EQ.3) THEN
454  DO 74 k = 1,3
455  rij2 = rij2 + (x(k,i3) - x(k,ntot))**2
456  vij2 = vij2 + (xdot(k,i3) - xdot(k,ntot))**2
457  rdot = rdot + (x(k,i3) - x(k,ntot))*
458  & (xdot(k,i3) - xdot(k,ntot))
459  74 CONTINUE
460  hi = 0.5*vij2 - (body(i3) + body(ntot))/sqrt(rij2)
461  IF (hi.GT.0.0.AND.rdot.GT.0.0) THEN
462  CALL hivel(i3)
463  END IF
464  END IF
465 *
466  IF (kz(30).GT.1) THEN
467  WRITE (6,75) time+toff, i3, i4, ntot, step3, step4,step(ntot)
468  75 FORMAT (' CHTERM: T I3 I4 NT DT ',f10.4,3i6,1p,3e9.1)
469  END IF
470 *
471 * Initialize second KS pair if separation is small.
472  IF (ks2.GT.0) THEN
473  i3 = 0
474 * Re-determine the global indices after updating in KSREG.
475  DO 80 i = ifirst,n
476  IF (name(i).EQ.name3) i3 = i
477  IF (name(i).EQ.name4) THEN
478  i4 = i
479  IF (i3.GT.0) go to 85
480  END IF
481  80 CONTINUE
482 *
483 * Define components and perform new regularization.
484  85 icomp = min(i3,i4)
485  jcomp = max(i3,i4)
486  CALL ksreg
487  time = time0
488  IF (sqrt(vx2)*vstar.GT.2.0*vstar) THEN
489  j1 = 2*npairs - 1
490  j2 = j1 + 1
491  a1 = -0.5*body(ntot)/h(npairs)
492  pb = days*a1*sqrt(abs(a1)/body(ntot))
493  vi2 = xdot(1,ntot)**2 + xdot(2,ntot)**2 + xdot(3,ntot)**2
494  vcm = sqrt(vi2)*vstar
495  WRITE (6,88) name(j1), name(j2), body(j1)*smu,
496  & body(j2)*smu, vcm, a1, pb
497  88 FORMAT (' CHAIN BINARY2 NM M1 M2 VCM A PB ',
498  & 2i7,2f6.1,f7.1,1p,2e10.2)
499  END IF
500  IF (kz(30).GT.1) THEN
501  WRITE (6,90) i3, i4, step(ntot), stepr(ntot),
502  & r(npairs), h(npairs), gamma(npairs)
503  90 FORMAT (' CHTERM: SECOND BINARY I3 I4 DT DTR R H G ',
504  & 2i5,1p,5e9.1)
505  END IF
506  END IF
507 *
508 * Replace any old perturbers that have been exchanged in KSREG.
509  IF (listc(1).GT.0.AND.listc(2).LT.ifirst) THEN
510  DO 92 j = ifirst,n
511  IF (name(j).EQ.namc2) THEN
512  listc(2) = j
513  ELSE IF (listc(1).GT.1.AND.name(j).EQ.namc3) THEN
514  listc(3) = j
515  END IF
516  92 CONTINUE
517  END IF
518 *
519 * Initialize the perturbers to improve derivatives (same locations).
520  nbc1 = listc(1) + 1
521  DO 95 l = 2,nbc1
522  j = listc(l)
523  IF (j.GT.n) go to 95
524  DO 94 k = 1,3
525  x0dot(k,j) = xdot(k,j)
526  94 CONTINUE
527  CALL fpoly1(j,j,0)
528  CALL fpoly2(j,j,0)
529  95 CONTINUE
530 *
531 * Check minimum two-body distance.
532  dminc = min(dminc,rcoll)
533 *
534 * Update net binary energy change.
535  chcoll = chcoll + cm(9)
536 *
537 * Update number of DIFSY calls, tidal dissipations & collision energy.
538  nstepc = nstepc + nstep1
539  ndiss = ndiss + ndiss1
540  nsync = nsync + max(isync,0)
541  ecoll = ecoll + ecoll1
542  egrav = egrav + ecoll1
543  e(10) = e(10) + ecoll1
544 *
545 * Reduce subsystem counter and initialize membership & internal energy.
546  nsub = nsub - 1
547  nch = 0
548  ech = 0.0
549 *
550 * Check for subsystem at last COMMON dump (no restart with NSUB > 0).
551  IF (nsub.EQ.0.AND.kz(2).GT.1) THEN
552  IF (time - tdump.LT.timec) THEN
553  tdump = time
554  CALL mydump(1,2)
555  END IF
556  END IF
557 *
558 * Set new value of TPREV (note: ensure TIME > TPREV in INTGRT).
559  100 IF (iterm.GE.0) tprev = time - 16.0*dt8
560 *
561  RETURN
562 *
563  END