Nbody6
 All Files Functions Variables
merge2.f
Go to the documentation of this file.
1  SUBROUTINE merge2
2 *
3 *
4 * Merging 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 *
13 *
14 * COMMON variables for merge procedure
15 * ************************************
16 *
17 * ------------------------------------------------------------------
18 * CM Component masses of merged binary (2nd binary in CM(3&4)).
19 * HM Binding energy per unit mass of merged binary.
20 * KSTARM Roche stage indicator (KSTAR for 1st c.m.; ghost is OK).
21 * NAMEG Name of the associated ghost component.
22 * NAMEM Name of new c.m. (< 0) for identification of merger index.
23 * NMERG Total number of mergers.
24 * NMERGE Current number of merged binaries (maximum is MMAX).
25 * UM Regularized coordinates of merged binary.
26 * UMDOT Regularized velocity of merged binary.
27 * VREL Relative velocity of merged binary components.
28 * XREL Relative coordinates of merged binary components.
29 * ------------------------------------------------------------------
30 *
31 *
32 * Ensure that the hierarchy is retained as primary KS pair.
33  IF (name(n+kspair).GT.0) THEN
34  k = kspair
35  kspair = jcomp - n
36  jcomp = n + k
37  END IF
38 *
39 * Select deepest level as primary in case of two hierarchies.
40  IF (name(n+kspair).LT.0.AND.name(jcomp).LT.0) THEN
41 * Note possibility of a quartet (or even quintet) and triple.
42  IF (name(jcomp).LT.name(n+kspair)) THEN
43  k = kspair
44  kspair = jcomp - n
45  jcomp = n + k
46  END IF
47  END IF
48 *
49 * Set pair index & c.m. of inner binary and save outer component.
50  ipair = kspair
51  i = n + ipair
52  jcomp1 = jcomp
53  icomp1 = i
54  i1 = 2*ipair - 1
55  i2 = i1 + 1
56  semi = -0.5*body(i)/h(ipair)
57  ecc2 = (1.0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(body(i)*semi)
58  ecc = sqrt(ecc2)
59 *
60 * Produce diagnostics for [[B,S],B] quintuplet or sextuplet system.
61  IF (name(jcomp).GT.nzero) THEN
62  ri2 = 0.0
63  DO 2 k = 1,3
64  ri2 = ri2 + (x(k,i) - rdens(k))**2
65  2 CONTINUE
66  ri = sqrt(ri2)
67  eb = 0.5*body(i1)*body(i2)/semi
68  eb = eb*float(n-npairs)/zkin
69  eb = min(eb,999.9d0)
70  zm = (body(i) + body(jcomp))*smu
71  CALL findj(i,jg,im)
72 * Check possible extension of next look-up time (skip GR case).
73  IF (kz(19).GE.3.AND.kz(27).LT.3) THEN
74  jlist(1) = i1
75  jlist(2) = i2
76  jlist(3) = 2*(jcomp - n) - 1
77  jlist(4) = jlist(3) + 1
78  jlist(5) = jg
79  CALL newtev(5,ix)
80  END IF
81  which1 = ' QUINTUPLET'
82  IF (jg.GT.n) which1 = ' SEXTUPLET '
83  WRITE (6,5) which1, time+toff, zm, name(i1), name(i2),
84  & name(jcomp), ri, ecc, eb, semi,pcrit,gamma(ipair)
85  5 FORMAT (/,' NEW',a11,' T MT NM1 NM2 NM3 RI E0 EB0 A0 PC G0',
86  & f10.2,f6.2,3i6,2f6.2,f6.1,1p,3e10.2)
87  END IF
88 *
89 * Check for [[B,S],S] quartet or [[B,B],S] quintuplet.
90  IF (name(jcomp).GT.0.AND.name(jcomp).LE.nzero) THEN
91  CALL findj(i,jg,im)
92  IF (kz(19).GE.3.AND.kz(27).LT.3) THEN
93  jlist(1) = i1
94  jlist(2) = i2
95  jlist(3) = jcomp
96  jlist(4) = jg
97  CALL newtev(4,ix)
98  END IF
99  zm = (body(i) + body(jcomp))*smu
100  IF (jg.LE.n.AND.jcomp.LE.n) THEN
101  WRITE (6,7) time+toff, zm, name(i1), name(i2),
102  & name(jcomp), ecc, semi, pcrit, gamma(ipair)
103  7 FORMAT (/,' NEW QUARTET T MT NM1 NMG NM3 E0 A0 PC G0 ',
104  & f9.2,f6.2,3i6,f6.2,1p,3e10.2)
105  ELSE
106  WRITE (6,8) time+toff, zm, name(i1), name(i2),
107  & name(jcomp), ecc, semi, pcrit, gamma(ipair)
108  8 FORMAT (/,' NEW QUINTUP2 T MT NM1 NMG NM3 E0 A0 PC G0',
109  & f10.2,f6.2,3i6,f6.2,1p,3e10.2)
110  END IF
111  END IF
112 *
113 * Include diagnostics for double triple as [[B,S],[B,S]].
114  IF (name(jcomp).LT.0) THEN
115  jpair = jcomp - n
116  j1 = 2*jpair - 1
117  CALL findj(i1,ji,im)
118  CALL findj(j1,jj,jm)
119  IF (kz(19).GE.3.AND.kz(27).LT.3) THEN
120  jlist(1) = i1
121  jlist(2) = i2
122  jlist(3) = j1
123  jlist(4) = j1 + 1
124  jlist(5) = ji
125  jlist(6) = jj
126  CALL newtev(6,ix)
127  END IF
128  ai = -0.5*(cm(1,im) + cm(2,im))/hm(im)
129  aj = -0.5*(cm(1,jm) + cm(2,jm))/hm(jm)
130  zm = (body(i) + body(jcomp))*smu
131  gx = max(gamma(ipair),gamma(jpair))
132  WRITE (6,10) time+toff, zm, name(i1), name(i2), name(j1),
133  & name(2*jpair), ecc, ai, aj, r(ipair), r(jpair),
134  & pcrit, gx
135  10 FORMAT (/,' NEW HITRIP T MT NM E0 AI AJ RI RJ PC GX ',
136  & f9.2,f6.2,4i6,f6.2,1p,6e10.2)
137  END IF
138 *
139 * Ensure correct coordinates & velocities.
140  CALL resolv(ipair,3)
141 *
142 * Check optional diagnostics for hierarchy.
143  IF ((kz(18).EQ.1.OR.kz(18).EQ.3).AND.kstar(i).LE.20) THEN
144  CALL hiarch(ipair)
145  END IF
146 *
147 * Increase merger counter and set current merger index.
148  nmerg = nmerg + 1
149  nmerge = nmerge + 1
150  imerge = nmerge
151 *
152 * Save component masses and evaluate reduced mass of inner binary.
153  cm(1,imerge) = body(2*ipair-1)
154  cm(2,imerge) = body(2*ipair)
155  cm(3,imerge) = 0.0d0
156  zmu = body(2*ipair-1)*body(2*ipair)/body(n+ipair)
157 *
158 * Set current energy of any second perturbed binary (from XVPRED).
159  IF (jcomp.GT.n) THEN
160  CALL xvpred(jcomp,-1)
161  jpair = jcomp - n
162  IF (list(1,2*jpair-1).GT.0) h(jpair) = ht
163 * Ensure that outer binary components are resolved.
164  CALL ksres(jpair,j1,j2,0.0d0)
165 * Enforce non-zero membership for potential correction in NBPOT.
166  list(1,2*jpair-1) = 1
167  ELSE
168  CALL xvpred(jcomp,-1)
169  END IF
170 *
171 * Save the neighbours for correction procedure.
172  nnb = list(1,i)
173  DO 20 l = 1,nnb
174  j = list(l+1,i)
175  jpert(l) = j
176  20 CONTINUE
177 *
178 * Retain basic KS variables for explicit restart at merge termination.
179  hm(imerge) = h(ipair)
180  DO 25 k = 1,4
181  um(k,imerge) = u(k,ipair)
182  umdot(k,imerge) = udot(k,ipair)
183  25 CONTINUE
184 *
185 * Save stellar type.
186  kstarm(imerge) = kstar(i)
187 *
188 * Predict outer component to highest order if not on the block.
189  IF (time - t0(jcomp1).GT.0.0d0) THEN
190  CALL xvpred(jcomp1,-1)
191  END IF
192 *
193 * Set temporary KS components for hierarchical binary.
194  icomp = 2*ipair - 1
195  jcomp = icomp + 1
196 *
197 * Obtain potential energy with respect to inner components.
198  jlist(1) = icomp
199  jlist(2) = jcomp
200  CALL nbpot(2,nnb,pot1)
201 *
202 * Save relative configuration.
203  DO 30 k = 1,3
204  xrel(k,imerge) = x(k,icomp) - x(k,jcomp)
205  vrel(k,imerge) = x0dot(k,icomp) - x0dot(k,jcomp)
206 * Initialize primary velocity of JCOMP1 (needed in KSIN2).
207  x0dot(k,jcomp1) = xdot(k,jcomp1)
208  30 CONTINUE
209 *
210 * Include interaction of inner c.m. & neighbours before making new KS.
211  icomp = icomp1
212  jlist(1) = icomp
213  CALL nbpot(1,nnb,pot2)
214 *
215 * Copy mass of intruder to second KS component.
216  body(jcomp) = body(jcomp1)
217 *
218 * Replace name of #JCOMP with c.m. name (temporary for diagnostics).
219  name2 = name(jcomp)
220  name(jcomp) = name(jcomp1)
221 *
222 * Form new c.m. and initialize KS variables (JPERT safe first call!).
223  jcomp = jcomp1
224  CALL ksin2(1)
225 *
226 * See whether modifications due to second binary are needed.
227  pot3 = 0.0d0
228  pot4 = 0.0d0
229  IF (jcomp1.LE.n) go to 50
230 *
231 * Initialize unperturbed ghost binary of outer component.
232  t0(2*jpair-1) = 1.0e+06
233  list(1,2*jpair-1) = 0
234 *
235 * Apply tidal correction for outer binary perturbers.
236  jlist(1) = 2*jpair - 1
237  jlist(2) = 2*jpair
238  CALL nbpot(2,nnb,pot3)
239  jlist(1) = jcomp1
240  CALL nbpot(1,nnb,pot4)
241 *
242 * Update the merger energy to maintain conservation.
243  eb1 = body(2*jpair-1)*body(2*jpair)*h(jpair)/body(jcomp1)
244  emerge = emerge + eb1
245 *
246 * Save component masses and initialize ghost components.
247  cm(3,imerge) = body(2*jpair-1)
248  cm(4,imerge) = body(2*jpair)
249  body(2*jpair-1) = 0.0d0
250  body(2*jpair) = 0.0d0
251  list(1,jcomp) = 0
252 *
253 * Remove ghost from all neighbour lists.
254  50 jlist(1) = jcomp1
255  icm = n + kspair
256  CALL nbrem(icm,1,nnb)
257 *
258 * Specify JCOMP1 as ghost of zero mass.
259  body(jcomp1) = 0.0d0
260 *
261 * Initialize integration variables to prevent spurious predictions.
262  DO 60 k = 1,3
263  x0dot(k,jcomp1) = 0.0d0
264  xdot(k,jcomp1) = 0.0d0
265  f(k,jcomp1) = 0.0d0
266  fdot(k,jcomp1) = 0.0d0
267  d2(k,jcomp1) = 0.0d0
268  d3(k,jcomp1) = 0.0d0
269  d2r(k,jcomp1) = 0.0d0
270  d3r(k,jcomp1) = 0.0d0
271  60 CONTINUE
272 *
273 * Set large value of T0 which avoids integration of ghost particle.
274  t0(jcomp1) = 1.0e+06
275 * Set large X0 & X to avoid perturber selection (no escape removal).
276  x0(1,jcomp1) = 1.0e+06
277  x(1,jcomp1) = 1.0e+06
278 *
279 * Initialize c.m. & KS polynomials.
280  icomp = 2*ipair - 1
281  jcomp = icomp + 1
282  CALL ksin2(2)
283 *
284 * Define large negative c.m. name for identification & termination.
285  name(icm) = name(icm) - 3*nzero
286 *
287 * Set c.m. & ghost names for merger identification (include escape).
288  namem(imerge) = name(icm)
289  nameg(imerge) = name(jcomp1)
290  name(jcomp) = name2
291 *
292 * Copy stability limit for termination test A(1 - E) < R0 in KSINT.
293  r0(ipair) = pcrit
294 *
295 * Update merger energy to maintain conservation.
296  dphi = (pot2 - pot1) + (pot4 - pot3)
297  emerge = emerge + zmu*hm(imerge) + dphi
298 *
299 * Set IPHASE < 0 for new time-step list in routine INTGRT.
300  iphase = -1
301 *
302  RETURN
303 *
304  END