Nbody6
 All Files Functions Variables
reset2.f
Go to the documentation of this file.
1  SUBROUTINE reset2
2 *
3 *
4 * Termination of double hierarchy.
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),iflag(mmax)
11  CHARACTER*11 which1
12  REAL*8 xx(3,3),vv(3,3)
13  SAVE istab
14  DATA istab /0/
15 *
16 *
17 * Set index of disrupted pair and save output parameters.
18  ipair = kspair
19  i = n + ipair
20  e1 = body(2*ipair-1)*body(2*ipair)*h(ipair)/body(i)
21  g1 = gamma(ipair)
22  r1 = r(ipair)
23  semi1 = -0.5*body(i)/h(ipair)
24  ecc2 = (1.0 - r1/semi1)**2 + tdot2(ipair)**2/(body(i)*semi1)
25  ecc = sqrt(ecc2)
26 *
27 * Locate current position in the merger table.
28  imerge = 0
29  DO 1 k = 1,nmerge
30  IF (namem(k).EQ.name(i)) imerge = k
31  1 CONTINUE
32 *
33 * Produce info on quintuplet or sextuplet (NAMEG has outer c.m. name).
34  IF (nameg(imerge).GT.nzero) THEN
35  ri2 = 0.0
36  DO 2 k = 1,3
37  ri2 = ri2 + (x(k,i) - rdens(k))**2
38  2 CONTINUE
39  ri = sqrt(ri2)
40  pmin = semi1*(1.0 - ecc)
41  eb = 0.5*body(2*ipair-1)*body(2*ipair)/semi1
42  eb = eb*float(n-npairs)/zkin
43  zm = body(i)*smu
44  which1 = ' QUINTUPLET'
45 * Determine previous merger index with NAME greater by 2*NZERO.
46  jm = 1
47  DO 3 k = 1,nmerge
48  IF (namem(k).EQ.name(i) + 2*nzero) jm = k
49  3 CONTINUE
50 * Use previous merger index to identify ghost from earlier level.
51  DO 4 j = ifirst,ntot
52  IF (name(j).EQ.nameg(jm)) jg = j
53  4 CONTINUE
54  IF (jg.GT.n) which1 = ' SEXTUPLET '
55  WRITE (6,5) which1, time+toff, zm, name(2*ipair-1),
56  & name(2*ipair), nameg(imerge), ri, ecc, eb, semi1,
57  & pmin, g1
58  5 FORMAT (' END',a11,' T MT NM1 NM2 NM3 RI E1 EB1 A1 PM G1 ',
59  & f9.2,f6.2,3i6,2f6.2,f6.1,1p,3e10.2)
60  END IF
61 *
62 * Check presence of [[B,S],S] quartet or [[B,B],S] quintuplet.
63  IF (nameg(imerge).GT.0.AND.nameg(imerge).LE.nzero) THEN
64  i1 = 2*ipair - 1
65  i2 = 2*ipair
66  CALL findj(i1,jg,im)
67  zm = body(i)*smu
68  IF (name(i2).LE.nzero) THEN
69  WRITE (6,9) time+toff, zm, name(i1), name(i2), name(jg),
70  & ecc, semi1, r1, g1
71  9 FORMAT (/,' END QUARTET T MT NM1 NMG NM3 E1 A1 R1 G1 ',
72  & f9.2,f6.2,3i6,f6.2,1p,3e10.2)
73  ELSE
74  WRITE (6,11) time+toff, zm, name(i1), name(i2), name(jg),
75  & ecc, semi1, r1, g1
76  11 FORMAT (/,' END QUINTUP2 T MT NM1 NMG NM3 E1 A1 R1 G1',
77  & f10.2,f6.2,3i6,f6.2,1p,3e10.2)
78  END IF
79  END IF
80 *
81 * Include diagnostics for double triple ([[B,S],[B,S]]).
82  IF (nameg(imerge).LT.0) THEN
83  i1 = 2*ipair - 1
84  j1 = 2*ipair
85 * Obtain current merger index (note NAMEM(JM) = NAMEG(JG)).
86  CALL findj(j1,jj,jg)
87  im = 1
88  jm = 1
89 * Determine original merger indices of the two inner binaries.
90  DO 16 k = 1,nmerge
91  IF (namem(k).EQ.name(i) + 2*nzero) im = k
92  IF (namem(k).EQ.nameg(jg)) jm = k
93  16 CONTINUE
94  ai = -0.5*(cm(1,im) + cm(2,im))/hm(im)
95  aj = -0.5*(cm(1,jm) + cm(2,jm))/hm(jm)
96  pmin = semi1*(1.0 - ecc)
97  zm = body(i)*smu
98  WRITE (6,18) time+toff, zm, name(i1), nameg(im), -namem(jm),
99  & nameg(jm), ecc, ai, aj, r(ipair), semi1, pmin,
100  & g1
101  18 FORMAT (/,' END HITRIP T MT NM E1 AI AJ R1 A1 PM G1 ',
102  & f9.2,f6.2,4i6,f6.2,1p,6e10.2)
103  END IF
104 *
105 * Check optional diagnostics for hierarchy.
106  IF ((kz(18).EQ.1.OR.kz(18).EQ.3).AND.kstar(imerge).LE.20) THEN
107  CALL hiarch(ipair)
108  END IF
109 *
110 * Save neighbours for correction procedure.
111  nnb = list(1,i) + 1
112  DO 6 l = 2,nnb
113  j = list(l,i)
114  jpert(l) = j
115  6 CONTINUE
116 *
117 * Ensure that c.m. coordinates are known to highest order.
118  CALL xvpred(i,0)
119 *
120 * Predict neighbour coordinates & velocities (XDOT used by FPOLY1).
121  DO 7 l = 2,nnb
122  j = jpert(l)
123  CALL xvpred(j,0)
124  7 CONTINUE
125 *
126 * Obtain current coordinates & velocities and specify KS components.
127  CALL resolv(ipair,2)
128  icomp = 2*ipair - 1
129  jcomp = icomp + 1
130 *
131 * Initialize mass, coordinates and velocities for new c.m. body.
132  body(i) = body(icomp)
133  DO 8 k = 1,3
134  x(k,i) = x(k,icomp)
135  xdot(k,i) = xdot(k,icomp)
136  x0dot(k,i) = xdot(k,icomp)
137  8 CONTINUE
138 *
139 * Add outer component to perturber list.
140  jpert(1) = jcomp
141 *
142 * Sum first part of potential energy correction due to tidal effect.
143  jlist(1) = icomp
144  CALL nbpot(1,nnb,pot1)
145 *
146 * Find correct location of ghost particle using identification name.
147  icm = i
148  jcomp1 = jcomp
149  DO 10 i = 1,ntot
150  IF (body(i).EQ.0.0d0.AND.name(i).EQ.nameg(imerge)) jcomp1 = i
151  10 CONTINUE
152 *
153 * Regularize two-body configuration if JCOMP1 cannot be identified.
154  IF (jcomp.EQ.jcomp1) THEN
155  WRITE (6,12) imerge, nameg(imerge), jcomp
156  12 FORMAT (/,5x,'WARNING! RESET2 JCOMP NOT IDENTIFIED ',
157  & ' IM =',i3,' NAMEG =',i6,' JCOMP =',i6)
158  go to 100
159  END IF
160 *
161 * Initialize basic variables for ghost and new c.m. (JCOMP -> JCOMP1).
162  j1 = jcomp1
163  j = jcomp
164  13 t0(j1) = time
165  body(j1) = body(j)
166  DO 14 k = 1,3
167  x(k,j1) = x(k,j)
168  x0(k,j1) = x(k,j)
169  xdot(k,j1) = xdot(k,j)
170  x0dot(k,j1) = xdot(k,j)
171  14 CONTINUE
172  IF (j.EQ.jcomp) THEN
173  j1 = icm
174  j = icomp
175  go to 13
176  END IF
177 *
178 * Restore masses, coordinates & velocities of inner binary.
179  body(icomp) = cm(1,imerge)
180  body(jcomp) = cm(2,imerge)
181  zm = -body(icomp)/(body(icomp) + body(jcomp))
182 *
183 * Begin with second component since ICOMP holds new c.m. variables.
184  i = jcomp
185  DO 20 kcomp = 1,2
186  DO 15 k = 1,3
187  x(k,i) = x(k,icomp) + zm*xrel(k,imerge)
188  x0dot(k,i) = x0dot(k,icomp) + zm*vrel(k,imerge)
189  xdot(k,i) = x0dot(k,i)
190 * Note that XDOT is only needed for improved polynomials of JCOMP.
191  15 CONTINUE
192  i = icomp
193  zm = body(jcomp)/(body(icomp) + body(jcomp))
194  20 CONTINUE
195 *
196 * Copy KS variables for inner binary (small TDOT2 near apo/peri).
197  i1 = 2*ipair - 1
198  t0(i1) = time
199  list(1,i1) = 1
200  h(ipair) = hm(imerge)
201  r(ipair) = 0.0d0
202  tdot2(ipair) = 0.0d0
203  DO 30 k = 1,4
204  u(k,ipair) = um(k,imerge)
205  u0(k,ipair) = u(k,ipair)
206  udot(k,ipair) = umdot(k,imerge)
207  r(ipair) = r(ipair) + u(k,ipair)**2
208  30 CONTINUE
209 *
210 * Add #JCOMP1 to all neighbour lists containing ICM.
211  jlist(1) = jcomp1
212  CALL nbrest(icm,1,nnb)
213 *
214 * Add #ICM to neighbour list of #JCOMP1.
215  jlist(1) = jcomp1
216  jpert(1) = icm
217  CALL nbrest(icm,1,1)
218 *
219 * Form new neighbour list for #ICM (NB! NBREST(JCOMP1,1,1) skips).
220  rs0 = rs(icm)
221  CALL nblist(icm,rs0)
222 *
223 * Get new list for #JCOMP1 (NNB = 1 from MERGE2/NBREM for ghost c.m.).
224  rs0 = rs(jcomp1)
225  CALL nblist(jcomp1,rs0)
226 *
227 * Initialize force polynomial for outer component using resolved c.m.
228  CALL fpoly1(jcomp1,jcomp1,0)
229  CALL fpoly2(jcomp1,jcomp1,0)
230 *
231 * Initialize c.m. polynomials and activate inner binary.
232  CALL ksin2(3)
233 *
234 * Restore original name of inner hierarchy (c.m. NAME set in MERGE2).
235  name(icm) = name(icm) + 2*nzero
236 *
237 * Locate position of inner binary in merger table (IM = 1 for safety).
238  im = 1
239  DO 35 k = 1,nmerge
240  IF (namem(k).EQ.name(icm)) im = k
241  35 CONTINUE
242 *
243 * Determine inclination between inner relative motion and outer orbit.
244  rv = 0.0
245  DO 40 k = 1,3
246  xx(k,1) = xrel(k,im)
247  xx(k,2) = 0.0
248  xx(k,3) = x(k,jcomp)
249  vv(k,1) = vrel(k,im)
250  vv(k,2) = 0.0
251  vv(k,3) = xdot(k,jcomp)
252  rv = rv + xrel(k,im)*vrel(k,im)
253  40 CONTINUE
254  CALL inclin(xx,vv,x(1,icm),xdot(1,icm),angle)
255 *
256 * Evaluate stability parameter from current elements (minimum ECC).
257  semi0 = -0.5*body(icomp)/hm(im)
258  semi = -0.5*body(icm)/h(ipair)
259  ecc2 = (1.0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(body(icm)*semi)
260  ecc1 = sqrt(ecc2)
261 * Form inner eccentricity from two-body elements.
262  rin = sqrt(xrel(1,im)**2 + xrel(2,im)**2 + xrel(3,im)**2)
263  ecc2 = (1.0 - rin/semi0)**2 + rv**2/(body(icomp)*semi0)
264  ecc = sqrt(ecc2)
265 *
266 * Evaluate the general stability function.
267  IF (ecc1.LT.1.0) THEN
268  nst = nstab(semi0,semi,ecc,ecc1,angle,cm(1,imerge),
269  & cm(2,imerge),body(jcomp))
270  IF (nst.EQ.0) THEN
271  pcrit = 0.98*semi*(1.0 - ecc1)
272  pcr = stability(cm(1,imerge),cm(2,imerge),body(jcomp),
273  & ecc,ecc1,angle)*semi0
274  istab = istab + 1
275  IF (istab.LT.50) THEN
276  WRITE (6,41) ecc, ecc1, semi0, semi, pcrit, pcr
277  41 FORMAT (' STABLE2 E E1 A A1 PCR PC99 ',
278  & 2f7.3,1p,4e10.2)
279  END IF
280  ELSE
281  pcrit = 1.02*semi*(1.0 - ecc1)
282  END IF
283  ELSE
284  pcrit = stability(cm(1,imerge),cm(2,imerge),body(jcomp),
285  & ecc,ecc1,angle)*semi0
286  END IF
287 *
288 * Set critical pericentre distance for stability check.
289  r0(ipair) = pcrit
290 *
291 * Form new force polynomials for perturbers inside 10*R1 (skip J > N).
292  nnb1 = list(1,icm) + 1
293  rs2 = (10.0*r1)**2
294  DO 46 l = 2,nnb1
295  j = list(l,icm)
296  rij2 = 0.0
297  DO 42 k = 1,3
298  rij2 = rij2 + (x(k,icm) - x(k,j))**2
299  42 CONTINUE
300  IF (rij2.LT.rs2.AND.j.NE.jcomp1.AND.j.LE.n) THEN
301  t0(j) = time
302  DO 45 k = 1,3
303  x0dot(k,j) = xdot(k,j)
304  45 CONTINUE
305  CALL fpoly1(j,j,0)
306  CALL fpoly2(j,j,0)
307  END IF
308  46 CONTINUE
309 *
310 * Rename perturber list for routine NBPOT.
311  jpert(1) = jcomp1
312 *
313 * Restore Roche stage indicator (any ghost c.m. is OK).
314  kstar(icm) = kstarm(imerge)
315 *
316 * See whether the outer component is a single or composite particle.
317  pot3 = 0.0d0
318  pot4 = 0.0d0
319  IF (jcomp1.LE.n) go to 50
320 *
321 * Restore component masses for outer binary.
322  jpair = jcomp1 - n
323  body(2*jpair-1) = cm(3,imerge)
324  body(2*jpair) = cm(4,imerge)
325 *
326 * Update look-up times & radii and check possible Roche condition.
327  IF (kz(19).GE.3) THEN
328  IF (kstar(jcomp1).GT.0.AND.kstar(jcomp1).LE.10) THEN
329  CALL trflow(jpair,dtr)
330  tev(jcomp1) = time + dtr
331  tmdot = min(tev(jcomp1),tmdot)
332  tmdot = min(tev(2*jpair),tmdot)
333  END IF
334  END IF
335 *
336 * Obtain coordinates & velocities of unperturbed binary components.
337  CALL resolv(jpair,1)
338 *
339 * Select new perturbers and initialize polynomials for KS motion.
340  CALL kslist(jpair)
341  CALL kspoly(jpair,1)
342 *
343 * Apply tidal correction for outer binary perturbers.
344  jlist(1) = 2*jpair - 1
345  jlist(2) = 2*jpair
346  CALL nbpot(2,nnb,pot3)
347  jlist(1) = jcomp1
348  CALL nbpot(1,nnb,pot4)
349 *
350 * Update the merger energy.
351  eb2 = body(2*jpair-1)*body(2*jpair)*h(jpair)/body(jcomp1)
352  emerge = emerge - eb2
353 *
354  e2 = e1/eb2
355  eb2 = eb2/be(3)
356  dp = pot4 - pot3
357  IF (kz(15).GT.1) THEN
358  WRITE (6,48) jpair, h(jpair), body(2*jpair-1),
359  & body(2*jpair), e2, eb2, r(jpair), gamma(jpair),
360  & dp
361  48 FORMAT (' END OUTER MERGER',i5,' H =',f7.1,' M =',2f7.4,
362  & ' E1 =',f6.3,' EB2 =',f6.3,' RB2 =',1pe8.1,
363  & ' G2 =',e8.1,' DP =',e8.1)
364  END IF
365 *
366 * Include interaction of body #ICOMP & JCOMP with perturbers.
367  50 jlist(1) = icomp
368  jlist(2) = jcomp
369  CALL nbpot(2,nnb,pot2)
370 *
371 * Form square of c.m. velocity correction due to tidal effects.
372 * VI2 = X0DOT(1,ICM)**2 + X0DOT(2,ICM)**2 + X0DOT(3,ICM)**2
373  dphi = (pot2 - pot1) + (pot4 - pot3)
374 * CORR = 1.0 + 2.0*DPHI/(BODY(ICM)*VI2)
375 * IF (CORR.LE.0.0D0) CORR = 0.0
376 *
377 * Adjust c.m. velocity by net tidal energy correction.
378 * DO 60 K = 1,3
379 * X0DOT(K,ICM) = SQRT(CORR)*X0DOT(K,ICM)
380 * 60 CONTINUE
381 *
382 * Modify the merger energy to maintain conservation.
383  eb = body(2*ipair-1)*body(2*ipair)*h(ipair)/body(icm)
384  emerge = emerge - eb + dphi
385 *
386  e1 = e1/eb
387  eb = eb/be(3)
388  IF (kz(15).GT.1) THEN
389  WRITE (6,65) imerge, time+toff, body(2*ipair-1),
390  & body(2*ipair), r1, semi1, eb, e1,
391  & gamma(ipair), g1, nnb-1
392  65 FORMAT (' END MERGE2',i3,' T =',f8.2,' M =',2f7.4,
393  & ' R1 =',1pe8.1,' A1 =',e8.1,' EB =',0pf6.3,
394  & ' E1 =',f6.3,' GB =',1pe8.1,' G =',0pf6.3,
395  & ' NB =',i3)
396  END IF
397 *
398 * Check Roche look-up time.
399  IF (kstar(icm).GT.0.AND.kstar(icm).LE.10) THEN
400  k = icm - n
401  CALL trflow(k,dtr)
402  tev(icm) = min(tev(icm),time + dtr)
403  tmdot = min(tev(icm),tmdot)
404  END IF
405 *
406 * Reduce merger counter and update tables (unless last or only pair).
407  70 nmerge = nmerge - 1
408  DO 80 l = imerge,nmerge
409  l1 = l + 1
410  hm(l) = hm(l1)
411  nameg(l) = nameg(l1)
412  namem(l) = namem(l1)
413  kstarm(l) = kstarm(l1)
414  DO 74 k = 1,3
415  xrel(k,l) = xrel(k,l1)
416  vrel(k,l) = vrel(k,l1)
417  74 CONTINUE
418  DO 75 k = 1,4
419  cm(k,l) = cm(k,l1)
420  um(k,l) = um(k,l1)
421  umdot(k,l) = umdot(k,l1)
422  75 CONTINUE
423  80 CONTINUE
424 *
425 * Examine merger list for possible escapers (retain up to 3 levels).
426  DO 90 l = 1,nmerge
427  DO 85 j = 1,npairs
428  IF (namem(l).EQ.name(n+j).OR.
429  & namem(l).EQ.name(n+j) + 2*nzero.OR.
430  & namem(l).EQ.name(n+j) + 4*nzero) go to 90
431  85 CONTINUE
432 * Remove tables for any merger not identified.
433  imerge = l
434  go to 70
435  90 CONTINUE
436 *
437 * Set IPHASE < 0 for new time-step list in INTGRT.
438  iphase = -1
439 *
440  100 RETURN
441 *
442  END