Nbody6
 All Files Functions Variables
start3.f
Go to the documentation of this file.
1  SUBROUTINE start3(ISUB)
2 *
3 *
4 * Initialization & restart of triple system.
5 * ------------------------------------------
6 *
7  include 'common6.h'
8  REAL*8 m
9  common/azreg/ time3,tmax,q(8),p(8),r1,r2,r3,energy,m(3),x3(3,3),
10  & xdot3(3,3),cm(10),c11,c12,c19,c20,c24,c25,
11  & nstep3,name3(3),kz15,kz27
12  common/close/ rij(4,4),rcoll,qperi,SIZE(4),ecoll3,ip(4)
13  common/azcoll/ rk(3),qk(8),pk(8),icall,icoll,ndiss3
14  common/clump/ bodys(ncmax,5),t0s(5),ts(5),steps(5),rmaxs(5),
15  & names(ncmax,5),isys(5)
16 *
17 *
18 * Decide between new run, termination or collision (= 0, > 0, < 0).
19  IF (isub.NE.0) THEN
20  iterm = isub
21  isub = iabs(isub)
22  go to 100
23  END IF
24 *
25 * Define sequence with dominant binary component as reference body.
26  name3(1) = jclose
27 * Perturber index from routine IMPACT (ICOMP & JCOMP set in KSTERM).
28  IF (body(icomp).LT.body(jcomp)) THEN
29  k = icomp
30  icomp = jcomp
31  jcomp = k
32  END IF
33  name3(2) = jcomp
34  name3(3) = icomp
35 *
36  DO 2 k = 1,7
37  cm(k) = 0.0d0
38  2 CONTINUE
39 *
40 * Transform to the local c.m. reference frame.
41  DO 4 l = 1,3
42  j = name3(l)
43  m(l) = body(j)
44  cm(7) = cm(7) + m(l)
45  DO 3 k = 1,3
46  x3(k,l) = x(k,j)
47  xdot3(k,l) = xdot(k,j)
48  cm(k) = cm(k) + m(l)*x3(k,l)
49  cm(k+3) = cm(k+3) + m(l)*xdot3(k,l)
50  3 CONTINUE
51  4 CONTINUE
52 *
53 * Set c.m. coordinates & velocities for triple system.
54  DO 5 k = 1,6
55  cm(k) = cm(k)/cm(7)
56  5 CONTINUE
57 *
58 * Specify initial conditions for three-body regularization.
59  DO 8 l = 1,3
60  DO 7 k = 1,3
61  x3(k,l) = x3(k,l) - cm(k)
62  xdot3(k,l) = xdot3(k,l) - cm(k+3)
63  7 CONTINUE
64  8 CONTINUE
65 *
66 * Calculate internal energy and include in total subsystem energy.
67  CALL trans3(0)
68  esub = esub + 0.5d0*energy
69 *
70 * Save global indices & particle attributes of subsystem.
71  DO 10 l = 1,3
72  j = name3(l)
73  jlist(l) = j
74  name3(l) = name(j)
75  SIZE(l) = radius(j)
76  ip(l) = kstar(j)
77  10 CONTINUE
78 *
79 * Make perturber list for potential energy correction.
80  icm = jlist(1)
81  go to 200
82 *
83 * Obtain potential energy of resolved subsystem & perturbers.
84  20 CALL nbpot(3,np,pot1)
85 *
86 * Define subsystem indicator (ISYS = 1, 2, 3 for triple, quad, chain).
87  isys(nsub+1) = 1
88 *
89 * Form ghosts and initialize c.m. motion in ICOMP (= JLIST(1) = I).
90  CALL subsys(3,cm)
91 *
92 * Include interaction of subsystem c.m. & perturbers for net effect.
93  CALL nbpot(1,np,pot2)
94 *
95 * Form square of c.m. velocity correction due to differential force.
96  vi2 = x0dot(1,icomp)**2 + x0dot(2,icomp)**2 + x0dot(3,icomp)**2
97  corr = 1.0 + 2.0*(pot2 - pot1)/(cm(7)*vi2)
98  IF (corr.LE.0.0d0) corr = 0.0
99 *
100 * Modify c.m. velocity by net tidal energy correction.
101  DO 30 k = 1,3
102  x0dot(k,icomp) = sqrt(corr)*x0dot(k,icomp)
103  30 CONTINUE
104 *
105 * Remove ghosts from perturber neighbour lists.
106  CALL nbrem(icm,3,np)
107 *
108 * Set maximum integration interval equal to c.m. step.
109  tmax = step(icomp)
110 *
111 * Copy total energy and output & capture option for routine TRIPLE.
112  cm(8) = be(3)
113  kz15 = kz(15)
114  kz27 = kz(27)
115 *
116 * Assign new subsystem index and begin triple regularization.
117  isub = nsub
118  ntrip = ntrip + 1
119  go to 180
120 *
121 * Prepare KS regularization and direct integration of third body.
122  100 imin = 1
123  IF (r2.LT.r1) imin = 2
124 *
125 * Specify global names of the KS candidates and least dominant body.
126  nam1 = name3(3)
127  nam2 = name3(imin)
128  nam3 = name3(3-imin)
129 *
130 * Identify current global indices by searching all single particles.
131  i = 0
132  DO 102 j = ifirst,n
133  IF (name(j).EQ.nam1) icomp = j
134  IF (name(j).EQ.nam2) jcomp = j
135  IF (name(j).EQ.nam3) i = j
136  102 CONTINUE
137 *
138 * Identify global index for c.m. body.
139  icm = i
140  IF (body(icomp).GT.0.0d0) icm = icomp
141  IF (body(jcomp).GT.0.0d0) icm = jcomp
142 *
143 * Quantize the elapsed interval since last step.
144  time2 = t0s(isub) + time3 - tprev
145  dt8 = (tblock - tprev)/8.0d0
146 *
147 * Adopt the nearest truncated step (at most 8 subdivisions).
148  dt2 = time2
149  IF (time2.GT.0.0d0) THEN
150  CALL stepk(dt2,dtn2)
151  dtn = nint(dtn2/dt8)*dt8
152  ELSE
153 * Choose negative step if pericentre time < TPREV (cf. iteration).
154  dt2 = -dt2
155  CALL stepk(dt2,dtn2)
156  dtn = -nint(dtn2/dt8)*dt8
157  END IF
158 *
159 * Update time for new polynomial initializations (also for CMBODY).
160  time = tprev + dtn
161  time = min(tblock,time)
162 *
163 * Predict current coordinates & velocities before termination.
164  CALL xvpred(icm,0)
165 *
166 * Copy c.m. coordinates & velocities.
167  DO 105 k = 1,3
168  cm(k) = x(k,icm)
169  cm(k+3) = xdot(k,icm)
170  105 CONTINUE
171 *
172 * Re-determine the perturber list.
173  go to 200
174 *
175 * Obtain potential energy for the c.m. subsystem & JPERT(NP).
176  110 jlist(1) = icm
177  CALL nbpot(1,np,pot1)
178 *
179 * Save global indices of three-body system.
180  jlist(1) = i
181  jlist(2) = icomp
182  jlist(3) = jcomp
183 *
184 * Set configuration pointers for escaper & KS candidates.
185  jlist(4) = 3 - imin
186  jlist(5) = 3
187  jlist(6) = imin
188 *
189 * Place new coordinates in the original locations.
190  DO 120 l = 1,3
191  j = jlist(l)
192 * Compare global name & subsystem name to restore the mass.
193  DO 112 k = 1,3
194  IF (name(j).EQ.names(k,isub)) THEN
195  body(j) = bodys(k,isub)
196  t0(j) = time
197  END IF
198  112 CONTINUE
199  ll = jlist(l+3)
200  DO 115 k = 1,3
201  x(k,j) = x3(k,ll) + cm(k)
202  115 CONTINUE
203  120 CONTINUE
204 *
205 * Obtain potential energy of subsystem & perturbers at the end.
206  CALL nbpot(3,np,pot2)
207 *
208 * Form square of c.m. velocity correction due to differential force.
209  vi2 = cm(4)**2 + cm(5)**2 + cm(6)**2
210  corr = 1.0 + 2.0*(pot2 - pot1)/(cm(7)*vi2)
211  IF (corr.LE.0.0d0) corr = 0.0
212 *
213 * Modify c.m. velocity by net tidal energy correction.
214  DO 122 k = 1,3
215  cm(k+3) = sqrt(corr)*cm(k+3)
216  122 CONTINUE
217 *
218 * Transform to global velocities using corrected c.m. velocity.
219  DO 130 l = 1,3
220  j = jlist(l)
221  ll = jlist(l+3)
222  DO 125 k = 1,3
223  xdot(k,j) = xdot3(k,ll) + 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 - 0.5d0*energy - ecoll3
252 *
253 * Ensure ICOMP < JCOMP for KS treatment.
254  k = icomp
255  icomp = min(icomp,jcomp)
256  jcomp = max(k,jcomp)
257 *
258 * Select consistent neighbours for single body and new c.m.
259  rs0 = rs(icm)
260  CALL nblist(i,rs0)
261  CALL nblist(icomp,rs0)
262 *
263 * Add any other ghosts to perturber list for replacement of #ICM.
264  DO 160 jsub = 1,nsub
265  DO 155 l = 1,4
266  IF (names(l,jsub).EQ.0) go to 155
267  DO 154 j = 1,n
268  IF (name(j).EQ.names(l,jsub).AND.body(j).LE.0.0d0) THEN
269  np = np + 1
270  jpert(np) = j
271  go to 155
272  END IF
273  154 CONTINUE
274  155 CONTINUE
275  160 CONTINUE
276 *
277 * Replace ICM in relevant neighbour lists by all subsystem members.
278  CALL nbrest(icm,3,np)
279 *
280 * Check for stellar collision (only needs coordinates & velocities).
281  IF (iterm.LT.0) THEN
282  jlist(1) = icomp
283  jlist(2) = jcomp
284  jlist(3) = i
285 * Initialize fictitious fourth body for general treatment.
286  jlist(4) = 0
287  go to 170
288  END IF
289 *
290 * Exclude the dominant interaction for c.m. approximation (large FDOT).
291  IF (r3**2.GT.cmsep2*min(r1,r2)**2) THEN
292  nnb = 1
293  ELSE
294  nnb = 3
295  END IF
296 *
297 * Set dominant F & FDOT on ICOMP & JCOMP for body #I in FPOLY2.
298  CALL fclose(icomp,nnb)
299  CALL fclose(jcomp,nnb)
300 *
301 * Initialize force polynomials & time-steps for body #I.
302  CALL fpoly1(i,i,0)
303  CALL fpoly2(i,i,0)
304 *
305 * Perform KS regularization of dominant components.
306  CALL ksreg
307 *
308 * Check minimum two-body distance.
309  dmin3 = min(dmin3,rcoll)
310 *
311 * Update net binary energy change.
312  sbcoll = sbcoll + cm(9)
313 *
314 * Update number of DIFSY calls, tidal dissipations & collision energy.
315  170 nstept = nstept + nstep3
316  ndiss = ndiss + ndiss3
317  ecoll = ecoll + ecoll3
318  e(10) = e(10) + ecoll3
319 *
320 * Check for subsystem at last COMMON dump (no restart with NSUB > 0).
321  IF (nsub.EQ.0.AND.kz(2).GE.1) THEN
322  IF (time - tdump.LT.time3) THEN
323  tdump = time
324  CALL mydump(1,2)
325  END IF
326  END IF
327 *
328 * Set phase indicator = -1 to ensure new time-step list in INTGRT.
329  180 iphase = -1
330 *
331  RETURN
332 *
333 * Form the current perturber list.
334  200 rp2 = rs(icm)**2
335  np = 0
336  nnb2 = list(1,icm) + 1
337 *
338 * Loop over all single particles & c.m. but skip subsystem members.
339  DO 210 l = 2,nnb2
340  j = list(l,icm)
341  rij2 = (x(1,j) - cm(1))**2 + (x(2,j) - cm(2))**2 +
342  & (x(3,j) - cm(3))**2
343  IF (rij2.LT.rp2) THEN
344  IF (j.EQ.i) go to 210
345  IF (j.EQ.icomp.OR.j.EQ.jcomp) go to 210
346  np = np + 1
347  jpert(np) = j
348  END IF
349  210 CONTINUE
350 *
351  IF (isub.EQ.0) go to 20
352  go to 110
353 *
354  END