Nbody6
 All Files Functions Variables
start4.f
Go to the documentation of this file.
1  SUBROUTINE start4(ISUB)
2 *
3 *
4 * Initialization & restart of four-body system.
5 * ---------------------------------------------
6 *
7  include 'common6.h'
8  REAL*8 m,mij
9  common/creg/ m(4),x4(12),xdot4(12),p(12),q(12),time4,energy,
10  & epsr2,xr(9),w(9),rr(6),ta(6),mij(6),cm(10),rmax4,
11  & tmax,ds,tstep,eps,nstep4,name4(4),kz15,kz27,nreg,nfn
12  common/config/ r2(4,4),i1,i2,i3,i4
13  common/close/ rij(4,4),rcoll,qperi,SIZE(4),ecoll3,ip(4)
14  common/ccoll/ qk(12),pk(12),icall,icoll,ndiss4
15  common/clump/ bodys(ncmax,5),t0s(5),ts(5),steps(5),rmaxs(5),
16  & names(ncmax,5),isys(5)
17 *
18 *
19 * Decide between new run, termination or collision (= 0, > 0, < 0).
20  IF (isub.NE.0) THEN
21  iterm = isub
22  isub = iabs(isub)
23  go to 100
24  END IF
25 *
26  DO 2 k = 1,7
27  cm(k) = 0.0d0
28  2 CONTINUE
29 *
30 * Transform to the local c.m. reference frame.
31  lk = 0
32  DO 4 l = 1,4
33  j = 2*npairs + l
34 * Copy four-body system from the first single particle locations.
35  m(l) = body(j)
36  cm(7) = cm(7) + m(l)
37  DO 3 k = 1,3
38  lk = lk + 1
39  x4(lk) = x(k,j)
40  xdot4(lk) = xdot(k,j)
41  cm(k) = cm(k) + m(l)*x4(lk)
42  cm(k+3) = cm(k+3) + m(l)*xdot4(lk)
43  3 CONTINUE
44  4 CONTINUE
45 *
46 * Set c.m. coordinates & velocities of four-body system.
47  DO 5 k = 1,6
48  cm(k) = cm(k)/cm(7)
49  5 CONTINUE
50 *
51 * Specify initial conditions for chain regularization.
52  lk = 0
53  DO 8 l = 1,4
54  DO 7 k = 1,3
55  lk = lk + 1
56  x4(lk) = x4(lk) - cm(k)
57  xdot4(lk) = xdot4(lk) - cm(k+3)
58  7 CONTINUE
59  8 CONTINUE
60 *
61 * Calculate internal energy and include in total subsystem energy.
62  CALL newsys(x4,xdot4,m,4,energy,gam)
63  esub = esub + energy
64 *
65 * Save global indices & particle attributes of subsystem.
66  DO 10 l = 1,4
67  j = 2*npairs + l
68  jlist(l) = j
69  name4(l) = name(j)
70  SIZE(l) = radius(j)
71  ip(l) = kstar(j)
72  10 CONTINUE
73 *
74 * Make perturber list for potential energy correction.
75  icm = jlist(1)
76  go to 200
77 *
78 * Obtain potential enegy of resolved subsystem & perturbers.
79  20 CALL nbpot(4,np,pot1)
80 *
81 * Define subsystem indicator (ISYS = 1, 2, 3 for triple, quad, chain).
82  isys(nsub+1) = 2
83 *
84 * Form ghosts and initialize c.m. motion in ICOMP (= JLIST(1)).
85  CALL subsys(4,cm)
86 *
87 * Include interaction of subsystem c.m. & perturbers for net effect.
88  CALL nbpot(1,np,pot2)
89 *
90 * Form square of c.m. velocity correction due to differential force.
91  vi2 = x0dot(1,icomp)**2 + x0dot(2,icomp)**2 + x0dot(3,icomp)**2
92  corr = 1.0 + 2.0*(pot2 - pot1)/(cm(7)*vi2)
93  IF (corr.LE.0.0d0) corr = 0.0
94 *
95 * Modify c.m. velocity by net tidal energy correction.
96  DO 30 k = 1,3
97  x0dot(k,icomp) = sqrt(corr)*x0dot(k,icomp)
98  30 CONTINUE
99 *
100 * Remove ghosts from perturber neighbour lists.
101  CALL nbrem(icm,4,np)
102 *
103 * Set maximum integration interval equal to c.m. step.
104  tmax = step(icomp)
105 *
106 * Copy total energy and output & capture options for routine QUAD.
107  cm(8) = be(3)
108  kz15 = kz(15)
109  kz27 = kz(27)
110 *
111 * Assign new subsystem index and begin four-body regularization.
112  isub = nsub
113  nquad = nquad + 1
114  go to 180
115 *
116 * Prepare KS regularization and direct integration of two bodies.
117  100 jlist(5) = name4(i1)
118  jlist(6) = name4(i2)
119  jlist(7) = name4(i3)
120  jlist(8) = name4(i4)
121 *
122 * Indentify current global index by searching all single particles.
123  DO 102 j = ifirst,n
124  DO 101 l = 1,4
125  IF (name(j).EQ.jlist(l+4)) THEN
126  jlist(l) = j
127  END IF
128  101 CONTINUE
129  102 CONTINUE
130 *
131 * Ensure ICOMP < JCOMP for KS regularization.
132  icomp = min(jlist(1),jlist(2))
133  jcomp = max(jlist(1),jlist(2))
134 *
135 * Identify global index of c.m. body.
136  icm = 0
137  DO 104 l = 1,4
138  j = jlist(l)
139  IF (body(j).GT.0.0d0) icm = j
140  104 CONTINUE
141 *
142 * Quantize the elapsed interval since last step.
143  time2 = t0s(isub) + time4 - 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  time = min(tblock,time)
161 *
162 * Predict current coordinates & velocities before termination.
163  CALL xvpred(icm,0)
164 *
165 * Coopy c.m. coordinates & velocities.
166  DO 105 k = 1,3
167  cm(k) = x(k,icm)
168  cm(k+3) = xdot(k,icm)
169  105 CONTINUE
170 *
171 * Re-determine the perturber list.
172  go to 200
173 *
174 * Obtain potential energy of the c.m. subsystem & JPERT(NP).
175  110 i = jlist(1)
176  jlist(1) = icm
177  CALL nbpot(1,np,pot1)
178 *
179 * Set configuration pointers for KS candidates & distant bodies.
180  jlist(5) = i1
181  jlist(6) = i2
182  jlist(7) = i3
183  jlist(8) = i4
184 *
185 * Place new coordinates in the original locations.
186  jlist(1) = i
187  lk = 0
188  DO 120 l = 1,4
189  j = jlist(l)
190 * Compare global name & subsystem name to restore the mass.
191  DO 112 k = 1,4
192  IF (name(j).EQ.names(k,isub)) THEN
193  body(j) = bodys(k,isub)
194  END IF
195  112 CONTINUE
196  ll = jlist(l+4)
197  DO 115 k = 1,3
198  lk = lk + 1
199  x(k,j) = x4(lk) + cm(k)
200  115 CONTINUE
201  120 CONTINUE
202 *
203 * Obtain potential energy of subsystem & perturbers at the end.
204  CALL nbpot(4,np,pot2)
205 *
206 * Form square of c.m. velocity correction due to differential force.
207  vi2 = cm(4)**2 + cm(5)**2 + cm(6)**2
208  corr = 1.0 + 2.0*(pot2 - pot1)/(cm(7)*vi2)
209  IF (corr.LE.0.0d0) corr = 0.0
210 *
211 * Modify c.m. velocity by net tidal energy correction.
212  DO 122 k = 1,3
213  cm(k+3) = sqrt(corr)*cm(k+3)
214  122 CONTINUE
215 *
216 * Transform to global velocities using corrected c.m. velocity.
217  lk = 0
218  DO 130 l = 1,4
219  j = jlist(l)
220  ll = jlist(l+4)
221  DO 125 k = 1,3
222  lk = lk + 1
223  xdot(k,j) = xdot4(lk) + cm(k+3)
224  x0dot(k,j) = xdot(k,j)
225  125 CONTINUE
226  130 CONTINUE
227 *
228 * Predict coordinates & velocities of perturbers to order FDOT.
229  DO 140 l = 1,np
230  j = jpert(l)
231  CALL xvpred(j,0)
232  140 CONTINUE
233 *
234 * Update subsystem COMMON variables unless last or only case.
235  IF (isub.LT.nsub) THEN
236  DO 150 l = isub,nsub
237  DO 145 k = 1,6
238  bodys(k,l) = bodys(k,l+1)
239  names(k,l) = names(k,l+1)
240  145 CONTINUE
241  t0s(l) = t0s(l+1)
242  ts(l) = ts(l+1)
243  steps(l) = steps(l+1)
244  rmaxs(l) = rmaxs(l+1)
245  isys(l) = isys(l+1)
246  150 CONTINUE
247  END IF
248 *
249 * Reduce subsystem counter and subtract internal binding energy.
250  nsub = nsub - 1
251  esub = esub - energy - ecoll3
252 *
253 * Select neighbours for single bodies and new c.m. (JLIST is *2).
254  rs0 = rs(icm)
255  j3 = jlist(3)
256  j4 = jlist(4)
257  CALL nblist(j3,rs0)
258  CALL nblist(j4,rs0)
259  CALL nblist(icomp,rs0)
260 *
261 * Add any other ghosts to perturber list for replacement of #ICM.
262  DO 160 jsub = 1,nsub
263  DO 155 l = 1,4
264  IF (names(l,jsub).EQ.0) go to 155
265  DO 154 j = 1,n
266  IF (name(j).EQ.names(l,jsub).AND.body(j).LE.0.0d0) THEN
267  np = np + 1
268  jpert(np) = j
269  go to 155
270  END IF
271  154 CONTINUE
272  155 CONTINUE
273  160 CONTINUE
274 *
275 * Replace ICM in perturber neighbour lists by all subsystem members.
276  CALL nbrest(icm,4,np)
277 *
278 * Check for stellar collision (only needs coordinates & velocities).
279  IF (iterm.LT.0) THEN
280  jlist(1) = icomp
281  jlist(2) = jcomp
282 *
283 * See whether relabelling is required (indices I1 - I4 still local).
284  IF (r2(i1,i4).LT.r2(i1,i3).OR.r2(i3,i4).LT.r2(i1,i3)) THEN
285  IF (r2(i1,i4).LT.r2(i3,i4)) THEN
286 * Switch body #I3 & I4 to give new dominant pair I1 & I3.
287  i = jlist(4)
288  jlist(4) = jlist(3)
289  jlist(3) = i
290  ELSE
291 * Set JLIST(5) < 0 to denote that body #I3 & I4 will be new KS pair.
292  jlist(5) = -1
293  END IF
294  END IF
295  go to 170
296  END IF
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(i3)
301  jlist(2) = jlist(i4)
302  nnb = 2
303  ELSE
304  nnb = 4
305  END IF
306 *
307 * Set dominant F & FDOT on body #ICOMP & JCOMP for #I3 & I4 in FPOLY2.
308  CALL fclose(icomp,nnb)
309  CALL fclose(jcomp,nnb)
310 *
311 * Specify global indices of least dominant bodies.
312  i3 = jlist(3)
313  i4 = jlist(4)
314 *
315 * Initialize force polynomials & time-steps for body #I3 & #I4.
316  CALL fpoly1(i3,i3,0)
317  CALL fpoly1(i4,i4,0)
318 *
319  CALL fpoly2(i3,i3,0)
320  CALL fpoly2(i4,i4,0)
321 *
322 * Perform KS regularization of dominant components (ICOMP < JCOMP).
323  CALL ksreg
324 *
325 * Check minimum two-body distance.
326  dmin4 = min(dmin4,rcoll)
327 *
328 * Update net binary energy change.
329  bbcoll = bbcoll + cm(9)
330 *
331 * Update number of DIFSY calls, tidal dissipations & collision energy.
332  170 nstepq = nstepq + nstep4
333  ndiss = ndiss + ndiss4
334  ecoll = ecoll + ecoll3
335  e(10) = e(10) + ecoll3
336 *
337 * Check for subsystem at last COMMON dump (no restart with NSUB > 0).
338  IF (nsub.EQ.0.AND.kz(2).GE.1) THEN
339  IF (time - tdump.LT.time4) THEN
340  tdump = time
341  CALL mydump(1,2)
342  END IF
343  END IF
344 *
345 * Set phase indicator = -1 to ensure new time-step list in INTGRT.
346  180 iphase = -1
347 *
348  RETURN
349 *
350 * Form the current perturber list.
351  200 rp2 = rs(icm)**2
352  np = 0
353  nnb2 = list(1,icm) + 1
354 *
355 * Loop over all single particles & c.m. but skip subsystem members.
356  DO 210 l = 2,nnb2
357  j = list(l,icm)
358  rij2 = (x(1,j) - cm(1))**2 + (x(2,j) - cm(2))**2 +
359  & (x(3,j) - cm(3))**2
360  IF (rij2.LT.rp2) THEN
361  DO 205 k = 1,4
362  IF (j.EQ.jlist(k)) go to 210
363  205 CONTINUE
364  np = np + 1
365  jpert(np) = j
366  END IF
367  210 CONTINUE
368 *
369  IF (isub.EQ.0) go to 20
370  go to 110
371 *
372  END