Nbody6
 All Files Functions Variables
slow.f
Go to the documentation of this file.
1  SUBROUTINE slow(Y)
2 *
3 *
4 * Slow-down treatment of chain binary.
5 * ------------------------------------
6 *
7  include 'commonc.h'
8  include 'common2.h'
9  LOGICAL kslow,kcoll,kcase
10  REAL*8 ksch,ksnew,vi(nmx3),vc(nmx3),rc1(3),rc2(3)
11  common/slow1/ tk2(0:nmx),ejump,ksch(nmx),kslow,kcoll
12  common/ccoll2/ qk(nmx4),pk(nmx4),rij(nmx,nmx),SIZE(nmx),vstar1,
13  & ecoll1,rcoll,qperi,istar(nmx),icoll,isync,ndiss1
14  common/ksave/ k1,k2
15  common/slow3/ gcrit,kz26
16  common/chreg/ timec,tmax,rmaxc,cm(10),namec(6),nstep1,kz27,kz30
17  SAVE
18 *
19 *
20 * Check for switching off slow-down at start of iteration.
21  if (gcrit.eq.0.0d0) then
22  ksnew = 1.0
23  go to 90
24  end if
25 *
26 * Set logical variable to avoid multiple copies of QK & PK.
27  kcase = .false.
28 *
29 * Perform perturbation check if slow-down not active.
30  IF (.NOT.kslow) THEN
31 * Determine chain index and largest inverse distance.
32  rm = 0.0
33  DO 1 i = 1,n-1
34  IF (rinv(i).GT.rm) THEN
35  rm = rinv(i)
36  i1 = i
37  END IF
38  1 CONTINUE
39 *
40 * Check carefully two possible binaries (eccentricity effect).
41  IF (n.GT.3) THEN
42  kcase = .true.
43 * Save QK & PK and copy current configuration for EREL & TRANSK.
44  DO 5 i = 1,n-1
45  ks = 4*(i - 1)
46  DO 4 j = 1,4
47  qk(ks+j) = q(ks+j)
48  pk(ks+j) = p(ks+j)
49  4 CONTINUE
50  5 CONTINUE
51 * Evaluate first semi-major axis from non-singular variables.
52  k1 = iname(i1)
53  k2 = iname(i1+1)
54  CALL erel(i1,eb,semi)
55 * Determine index of second smallest separation.
56  ri2 = 0.0
57  DO 10 i = 1,n-1
58  IF (i.NE.i1.AND.rinv(i).GT.ri2) THEN
59  i2 = i
60  ri2 = rinv(i)
61  END IF
62  10 CONTINUE
63 * Obtain second semi-major axis (errors only affect decision-making).
64  k1 = iname(i2)
65  k2 = iname(i2+1)
66  CALL erel(i2,eb,semi2)
67 * Switch to #I2 if second binary is smaller or first pair is not bound.
68  IF (semi2.GT.0.0) THEN
69  IF (semi2.LT.semi.OR.semi.LT.0.0) THEN
70  i1 = i2
71  semi = semi2
72  rm = rinv(i2)
73  END IF
74 * Exit in case of two hyperbolic two-body motions.
75  ELSE IF (semi.LT.0.0) THEN
76  go to 100
77  END IF
78  END IF
79 *
80 * Sum the perturbing forces m/r^3 next to #i1 (two terms if i1 = 2).
81  sum = 0.0
82  do 15 i = 1,n-1
83  if (iabs(i-i1).eq.1) then
84  j = i
85  if (i.gt.i1) j = i + 1
86  sum = sum + mc(j)*rinv(i)**3
87  end if
88  15 continue
89 *
90 * Include one more contribution for two consecutive perturbers.
91  if (i1.eq.1) then
92  lj = 3*i1
93  do k = 1,3
94  rc2(k) = xc(k+lj+3) + xc(k+lj)
95  end do
96 * Add perturbation from second subsequent member (i = i1 + 2).
97  rj = sqrt(rc2(1)**2 + rc2(2)**2 + rc2(3)**2)
98  j = i1 + 3
99  sum = sum + mc(j)/rj**3
100  else if (i1.ge.3) then
101  lj = 3*(i1 - 2)
102  do k = 1,3
103  rc2(k) = xc(k+lj-3) + xc(k+lj)
104  end do
105 * Add the previous perturbation (neglected in do 15 loop).
106  rj = sqrt(rc2(1)**2 + rc2(2)**2 + rc2(3)**2)
107  j = i1 - 2
108  sum = sum + mc(j)/rj**3
109  end if
110 *
111 * Skip further search if relative perturbation exceeds limit.
112  mb = mc(i1) + mc(i1+1)
113  gamma = 2.0*sum/(mb*rm**3)
114  IF (gamma.GT.gcrit) THEN
115  go to 100
116  END IF
117  END IF
118 *
119 * Determine chain index for slow-down binary (otherwise from above).
120  IF (kslow.AND.n.GT.3) THEN
121  DO 20 i = 1,n-1
122  IF (ksch(i).GT.1.0d0) i1 = i
123  20 CONTINUE
124 * See whether a closer particle pair is present (factor of 2).
125  i2 = 0
126  r1 = 1.0/rinv(i1)
127  DO 60 i = 1,n-1
128  ri2 = 1.0/rinv(i)
129  IF (i.NE.i1.AND.r1.LT.0.5*ri2) THEN
130  i2 = i
131  IF (r1.LT.0.1*semi.AND.ri2.LT.0.1*semi) THEN
132  go to 100
133  END IF
134 * Compare closest separation with current slow-down binary.
135  IF (ri2.LT.0.5*semi.AND.r1.GT.0.5*semi) THEN
136 * Evaluate semi-major axis directly (cf. small RI2 in EREL).
137  l = 3*(n-2)
138  DO 25 k = 1,3
139  vi(k) = -wc(k)/mc(1)
140  vi(l+k+3) = wc(l+k)/mc(n)
141  25 CONTINUE
142  DO 35 ii = 2,n-1
143  l = 3*(ii-1)
144  DO 30 k = 1,3
145  vi(l+k) = (wc(l+k-3) - wc(l+k))/mc(ii)
146  30 CONTINUE
147  35 CONTINUE
148  DO 40 j = 1,3*(n-1)
149  vc(j) = vi(j+3) - vi(j)
150  40 CONTINUE
151  l = 3*(i1-1)
152  r2 = xc(l+1)**2 + xc(l+2)**2 + xc(l+3)**2
153  w2 = vc(l+1)**2 + vc(l+2)**2 + vc(l+3)**2
154  semi = 2.0d0/r1 - w2/(mc(i1) + mc(i1+1))
155  semi = 1.0d0/semi
156  go to 80
157  ELSE IF (ri2.LT.min(semi,r1).AND.r1.GT.0.5*semi) THEN
158 * Determine second binary by regular expression (R1 not too small).
159  DO 50 ii = 1,n-1
160  ks = 4*(ii - 1)
161  DO 45 j = 1,4
162  qk(ks+j) = q(ks+j)
163  pk(ks+j) = p(ks+j)
164  45 CONTINUE
165  50 CONTINUE
166  k1 = iname(i2)
167  k2 = iname(i2+1)
168  CALL erel(i2,eb2,semi2)
169  IF (semi2.GT.0.0.AND.semi2.LT.semi) THEN
170  ksnew = 1.0d0
171  go to 90
172  ELSE
173 * Continue with the current binary (ie. small change in perturbation).
174  go to 80
175  END IF
176  END IF
177  END IF
178  60 CONTINUE
179  END IF
180 *
181 * Obtain regular semi-major axis for missing cases (including N = 3).
182  IF (.NOT.kcase) THEN
183 * Save QK & PK and copy current configuration for EREL & TRANSK.
184  DO 70 i = 1,n-1
185  ks = 4*(i - 1)
186  DO 65 j = 1,4
187  qk(ks+j) = q(ks+j)
188  pk(ks+j) = p(ks+j)
189  65 CONTINUE
190  70 CONTINUE
191 *
192 * Evaluate semi-major axis from non-singular variables.
193  k1 = iname(i1)
194  k2 = iname(i1+1)
195  CALL erel(i1,eb,semi)
196 * Exit if no current binary (set KSLOW = .false. just in case).
197  IF (semi.LE.0.0d0) THEN
198  kslow = .false.
199  tk2(0) = 0.0
200  go to 100
201  END IF
202  END IF
203 *
204 * Check for switching to smaller binary (exchange leads to escape).
205  IF (kslow.AND.n.GT.3) THEN
206  IF (i2.GT.0) THEN
207 * Evaluate second semi-major axis (K1 & K2 can be over-written).
208  k1 = iname(i2)
209  k2 = iname(i2+1)
210  CALL erel(i2,eb2,semi2)
211  IF (semi2.GT.0.0.AND.semi2.LT.semi) THEN
212 * Switch off the present pair to prepare re-activation.
213  ksnew = 1.0
214  go to 90
215  END IF
216  END IF
217  END IF
218 *
219 * Sum the perturbations (on either side and non-dominant terms).
220  80 IF (kslow) THEN
221  sum = 0.0
222  do i = 1,n-1
223 *
224 * Include full perturbation (non-symmetrical i1 for n = 4 or n > 4).
225  if ((n.eq.4.and.i1.ne.2).or.n.gt.4) then
226  do k = 1,3
227  rc1(k) = 0.0
228  end do
229 * Save vector sum on either side of #i1 excluding closest neighbour.
230  if (i.lt.i1-1) then
231  do j = i+1,i1-1
232  lj = 3*(j-1)
233  do k = 1,3
234  rc1(k) = rc1(k) + xc(k+lj)
235  end do
236  end do
237 * Check alternative case of subsequent distant members (i > i1 + 1).
238  else if (i.gt.i1+1) then
239  do j = i1+1,i-1
240  lj = 3*(j-1)
241  do k = 1,3
242  rc1(k) = rc1(k) + xc(k+lj)
243  end do
244  end do
245  end if
246  end if
247 *
248 * Set chain offset and mass reference index.
249  l = 3*(i-1)
250  j = i
251  if (i.gt.i1) j = i + 1
252 * Use actual separation if perturber is next to dominant binary.
253  if (iabs(i-i1).eq.1) then
254  sum = sum + mc(j)*rinv(i)**3
255  else if (i.ne.i1) then
256 * Include chain vector to yield full distance to binary.
257  do k = 1,3
258  rc2(k) = rc1(k) + xc(k+l)
259  end do
260 * Add contribution from more distant member.
261  rj = sqrt(rc2(1)**2 + rc2(2)**2 + rc2(3)**2)
262  sum = sum + mc(j)/rj**3
263  end if
264  end do
265  END IF
266 *
267 * Form relative perturbation at maximum apocentre.
268  rap = 2.0*semi
269  pert = 2.0*sum*rap**3/(mc(i1) + mc(i1+1))
270 *
271 * Specify the slow-down factor.
272  IF (pert.LT.gcrit) THEN
273  ksnew = sqrt(gcrit/pert)
274  ELSE
275  ksnew = 1.0
276  END IF
277 *
278 * --------------------------------------------------------------
279 * Implement discrete scheme (suppressed).
280 * i = i1
281 * rat = ksch(i)
282 * rat = ksnew/rat
283 * Check for significant changes (by 2) in slow-down factor.
284 * if (rat.gt.0.5.and.rat.le.2.0) then
285 * ksnew = Ksch(i)
286 * else if (rat.gt.2.0) then
287 * ksnew = 2*Ksch(i)
288 * Allow an extra factor of 32 at initialization.
289 * if (.not.KSLOW) then
290 * if (rat.ge.64.0) ksnew = 32*ksnew
291 * end if
292 * else if (rat.le.0.5) then
293 * ksnew = Ksch(i)/2
294 * if (rat.le.0.125) ksnew = Ksch(i)/4
295 * ksnew = max(ksnew,1.0D0)
296 * end if
297 * --------------------------------------------------------------
298 *
299 * Check slow-down switch and include any change in binding energy.
300  90 if (ksnew.ne.ksch(i1)) then
301  if (ksnew.eq.1.0d0) then
302  kslow = .false.
303  else
304  kslow = .true.
305  end if
306 * Add change in binary energy and save new slow-down index.
307  eb = -0.5d0*mc(i1)*mc(i1+1)/semi
308  deb = eb*(1.0/ksnew - 1.0/ksch(i1))
309  ksch(i1) = ksnew
310  else
311  deb = 0.0
312  end if
313 *
314 * Accumulate total energy difference due to slow-down.
315  ejump = ejump + deb
316 * Ensure zero EJUMP without slow-down to avoid small residual.
317  IF (.NOT.kslow) ejump = 0.0
318 *
319 * Form modified mass factors.
320  do i = 1,n
321  tk1(i) = -1.0/mc(i)
322  end do
323  do i = 1,n-1
324  if (ksch(i).ne.1.0d0) then
325  tk1(i) = tk1(i)/ksch(i)
326  tk1(i+1) = tk1(i+1)/ksch(i)
327  end if
328  end do
329  DO i = 1,n-1
330  tkk(i) = 0.5d0*(-tk1(i) - tk1(i+1))
331  mkk(i) = mc(i)*mc(i+1)/ksch(i)
332  END DO
333  do i = 1,n-1
334  m12 = mc(i) + mc(i+1)
335  dt12 = 0.5d0*(1.0d0 - 1.0d0/ksch(i))/m12
336  if (i.gt.1) tkk(i-1) = tkk(i-1) + dt12
337  if (i.lt.n-1) tkk(i+1) = tkk(i+1) + dt12
338  if (i.gt.1.and.i.lt.n-1) tk2(i) = -2.0d0*dt12
339  end do
340  tk2(0) = 0.0
341  tk2(n) = 0.0
342 *
343  100 RETURN
344 *
345  END