Nbody6
 All Files Functions Variables
nbint.f
Go to the documentation of this file.
1  SUBROUTINE nbint(I,IKS,IR,XI,XIDOT)
2 *
3 *
4 * Irregular integration.
5 * ----------------------
6 *
7  include 'common6.h'
8  common/chainc/ xc(3,ncmax),uc(3,ncmax),bodyc(ncmax),ich,
9  & listc(lmax)
10  REAL*8 xi(3),xidot(3),firr(3),freg(3),fd(3),fdum(3),dv(3)
11  SAVE tcall,istart
12  DATA istart /0/
13 *
14 *
15 * Initialize binary search time to current TTOT at start/restart.
16  IF (istart.EQ.0) THEN
17  istart = 1
18  tcall = ttot
19  END IF
20 *
21 * Check regularization criterion for single particles.
22  IF (step(i).LT.dtmin.AND.i.LE.n) THEN
23 * See whether dominant body can be regularized.
24  IF (iks.EQ.0) THEN
25  CALL search(i,iks)
26  END IF
27  END IF
28 *
29 * Include close encounter search for low-eccentric massive binaries.
30  IF (iks.EQ.0.AND.step(i).LT.8.0*dtmin.AND.tcall.LT.ttot) THEN
31 * Consider massive single bodies in absence of subsystems.
32  IF (i.LE.n.AND.body(i).GT.2.0*bodym.AND.nsub.EQ.0) THEN
33 *
34 * Obtain two-body elements and relative perturbation.
35  tcall = ttot + 0.01
36  jmin = 0
37  CALL orbit(i,jmin,semi,ecc,gi)
38 *
39  eb = -0.5*body(i)*body(jmin)/semi
40  IF (eb.LT.ebh.AND.gi.LT.0.25.AND.jmin.GE.ifirst) THEN
41  apo = semi*(1.0 + ecc)
42 * Check eccentricity (cf. max perturbation) and neighbour radius.
43  IF (ecc.LT.0.5.AND.apo.LT.0.02*rs(i)) THEN
44 * WRITE (6,3) NAME(I), NAME(JMIN), ECC, SEMI, EB
45 * 3 FORMAT (' KS TRY: NAM E A EB ',
46 * & 2I6,F7.3,1P,2E10.2)
47  iks = iks + 1
48  icomp = i
49  jcomp = jmin
50  tcall = ttot - 1.0d-04
51  END IF
52  END IF
53  END IF
54  END IF
55 *
56 * Obtain total force & first derivative.
57  DO 5 k = 1,3
58  xi(k) = x(k,i)
59  xidot(k) = xdot(k,i)
60  firr(k) = 0.0d0
61  fd(k) = 0.0d0
62  5 CONTINUE
63 *
64 * Assume small mass at centre for special case of no neighbours.
65  nnb0 = list(1,i)
66  IF (nnb0.EQ.0) THEN
67  ri2 = xi(1)**2 + xi(2)**2 + xi(3)**2
68  fij = 0.01*bodym/(ri2*sqrt(ri2))
69  rdot = 3.0*(xi(1)*xidot(1) + xi(2)*xidot(2) +
70  & xi(3)*xidot(3))/ri2
71  DO 10 k = 1,3
72  firr(k) = -fij*xi(k)
73  fd(k) = -(xidot(k) - rdot*xi(k))*fij
74  10 CONTINUE
75  IF (i.GT.n) ipair = i - n
76  go to 70
77  END IF
78 *
79 * Choose force loop for single particle or regularized c.m. body.
80  IF (i.LE.n) go to 20
81 *
82 * Set KS pair index.
83  ipair = i - n
84  i2 = 2*ipair
85  i1 = i2 - 1
86 *
87 * Adopt c.m. approximation for small total perturbation.
88  IF (list(1,i1).GT.0) THEN
89 * Obtain irregular force on perturbed c.m. body (including any chain).
90  CALL cmfirr(i,ipair,xi,xidot,firr,fd)
91  go to 70
92  END IF
93 *
94 * Copy c.m. coordinates & velocities for rare unperturbed intruder.
95  DO 15 k = 1,3
96  x(k,i1) = xi(k)
97  x(k,i2) = xi(k)
98  xdot(k,i1) = xidot(k)
99  xdot(k,i2) = xidot(k)
100  15 CONTINUE
101 *
102 * Set neighbour number & list index of the last single particle.
103  20 nnb1 = nnb0 + 1
104  nnb2 = nnb1
105  25 IF (list(nnb2,i).LE.n) go to 30
106  nnb2 = nnb2 - 1
107  IF (nnb2.GT.1) go to 25
108 * Include special case of only c.m. neighbours.
109  go to 40
110 *
111 * Sum over single particles (unperturbed case included).
112  30 DO 35 l = 2,nnb2
113  k = list(l,i)
114  a1 = x(1,k) - xi(1)
115  a2 = x(2,k) - xi(2)
116  a3 = x(3,k) - xi(3)
117  dv(1) = xdot(1,k) - xidot(1)
118  dv(2) = xdot(2,k) - xidot(2)
119  dv(3) = xdot(3,k) - xidot(3)
120  rij2 = a1*a1 + a2*a2 + a3*a3
121 *
122  dr2i = 1.0/rij2
123  dr3i = body(k)*dr2i*sqrt(dr2i)
124  drdv = 3.0*(a1*dv(1) + a2*dv(2) + a3*dv(3))*dr2i
125 *
126  firr(1) = firr(1) + a1*dr3i
127  firr(2) = firr(2) + a2*dr3i
128  firr(3) = firr(3) + a3*dr3i
129  fd(1) = fd(1) + (dv(1) - a1*drdv)*dr3i
130  fd(2) = fd(2) + (dv(2) - a2*drdv)*dr3i
131  fd(3) = fd(3) + (dv(3) - a3*drdv)*dr3i
132  35 CONTINUE
133 *
134 * See whether any c.m. neighbours should be included.
135  IF (nnb2.EQ.nnb1) go to 60
136 *
137  40 nnb3 = nnb2 + 1
138 * Set index for distinguishing c.m. or resolved components.
139  kdum = 0
140 *
141 * Sum over regularized c.m. neighbours.
142  DO 50 l = nnb3,nnb1
143  k = list(l,i)
144  a1 = x(1,k) - xi(1)
145  a2 = x(2,k) - xi(2)
146  a3 = x(3,k) - xi(3)
147  dv(1) = xdot(1,k) - xidot(1)
148  dv(2) = xdot(2,k) - xidot(2)
149  dv(3) = xdot(3,k) - xidot(3)
150  rij2 = a1*a1 + a2*a2 + a3*a3
151 *
152 * See whether c.m. approximation applies (ignore unperturbed case).
153  j = k - n
154  kdum = 2*j - 1
155  IF (rij2.GT.cmsep2*r(j)**2.OR.list(1,kdum).EQ.0) go to 48
156 *
157 * KDUM = 2*J - 1
158  k = kdum
159 * Sum over individual components of pair #J.
160  45 a1 = x(1,k) - xi(1)
161  a2 = x(2,k) - xi(2)
162  a3 = x(3,k) - xi(3)
163  dv(1) = xdot(1,k) - xidot(1)
164  dv(2) = xdot(2,k) - xidot(2)
165  dv(3) = xdot(3,k) - xidot(3)
166  rij2 = a1*a1 + a2*a2 + a3*a3
167 *
168 * Adopt c.m. approximation outside the effective perturber sphere.
169  48 dr2i = 1.0/rij2
170  dr3i = body(k)*dr2i*sqrt(dr2i)
171  drdv = 3.0*(a1*dv(1) + a2*dv(2) + a3*dv(3))*dr2i
172 *
173  firr(1) = firr(1) + a1*dr3i
174  firr(2) = firr(2) + a2*dr3i
175  firr(3) = firr(3) + a3*dr3i
176  fd(1) = fd(1) + (dv(1) - a1*drdv)*dr3i
177  fd(2) = fd(2) + (dv(2) - a2*drdv)*dr3i
178  fd(3) = fd(3) + (dv(3) - a3*drdv)*dr3i
179  IF (k.EQ.kdum) THEN
180  k = k + 1
181  go to 45
182  END IF
183  50 CONTINUE
184 *
185 * Include differential force treatment for regularized subsystem.
186  60 IF (nch.GT.0) THEN
187 * Distinguish between chain c.m. and any other particle.
188  IF (name(i).EQ.0) THEN
189  CALL chfirr(i,0,xi,xidot,firr,fd)
190  ELSE
191 * See if chain perturber list contains body #I.
192  np1 = listc(1) + 1
193  DO 65 l = 2,np1
194  j = listc(l)
195  IF (j.GT.i) go to 70
196  IF (j.EQ.i) THEN
197  CALL fchain(i,0,xi,xidot,firr,fd)
198  go to 70
199  END IF
200  65 CONTINUE
201  END IF
202  END IF
203 *
204 * Check option for external tidal field using predicted FREG.
205  70 dt = time - t0(i)
206  IF (kz(14).GT.0) THEN
207  DO 75 k = 1,3
208  freg(k) = fr(k,i) + frdot(k,i)*dt
209  75 CONTINUE
210  CALL xtrnlf(xi,xidot,firr,freg,fd,fdum,0)
211  END IF
212 *
213 * Include the corrector and set new T0, F, FDOT, D1, D2 & D3.
214  dtsq = dt**2
215  dt6 = 6.0d0/(dt*dtsq)
216  dt2 = 2.0d0/dtsq
217  dtsq12 = one12*dtsq
218  dt13 = one3*dt
219  t0(i) = time
220 *
221  DO 80 k = 1,3
222  df = fi(k,i) - firr(k)
223  fid = fidot(k,i)
224  sum = fid + fd(k)
225  at3 = 2.0d0*df + dt*sum
226  bt2 = -3.0d0*df - dt*(sum + fid)
227 *
228  x0(k,i) = xi(k) + (0.6d0*at3 + bt2)*dtsq12
229  x0dot(k,i) = xidot(k) + (0.75d0*at3 + bt2)*dt13
230 *
231 * X0(K,I) = X(K,I)
232 * X0DOT(K,I) = XDOT(K,I)
233 *
234  fi(k,i) = firr(k)
235  fidot(k,i) = fd(k)
236 * Use total force for irregular step (cf. Makino & Aarseth PASJ, 1992).
237  fdum(k) = firr(k) + fr(k,i)
238 *
239  d0(k,i) = firr(k)
240  d1(k,i) = fd(k)
241  d2(k,i) = (3.0d0*at3 + bt2)*dt2
242  d3(k,i) = at3*dt6
243 * NOTE: These are real derivatives!
244  80 CONTINUE
245 *
246 * Specify new time-step by standard criterion (STEPI version not good).
247  ttmp = tstep(fdum,fd,d2(1,i),d3(1,i),etai)
248  dt0 = ttmp
249 *
250 * IF (I.GT.N) THEN
251 * Check for hierarchical configuration but exclude small perturbations.
252 * IF (H(IPAIR).LT.-ECLOSE.AND.KZ(36).GT.0) THEN
253 * IF (GAMMA(IPAIR).GT.1.0E-04) THEN
254 * CALL KEPLER(I,TTMP)
255 * DT0 = TTMP
256 * END IF
257 * END IF
258 * END IF
259 *
260 * Include convergence test for large step (cf. Makino, Ap.J. 369, 200).
261  IF (ttmp.GT.stepj.AND.n.GT.1000) THEN
262  dv2 = 0.0
263  f2 = 0.0
264  DO 85 k = 1,3
265  dv2 = dv2 + (xdot(k,i) - x0dot(k,i))**2
266  f2 = f2 + firr(k)**2
267  85 CONTINUE
268 * Employ Jun's criterion to avoid over-shooting (cf. Book, 2.16).
269  dtj = step(i)*(1.0d-06*step(i)**2*f2/dv2)**0.1
270  ttmp = min(ttmp,dtj)
271  END IF
272 *
273 * Select discrete value (increased by 2, decreased by 2 or unchanged).
274  IF (ttmp.GT.2.0*step(i)) THEN
275  IF (dmod(time,2.0*step(i)).EQ.0.0d0) THEN
276  ttmp = min(2.0*step(i),smax)
277  ELSE
278  ttmp = step(i)
279  END IF
280  ELSE IF (ttmp.LT.step(i)) THEN
281  ttmp = 0.5*step(i)
282  IF (ttmp.GT.dt0) THEN
283  ttmp = 0.5*ttmp
284  END IF
285  ELSE
286  ttmp = step(i)
287  END IF
288 *
289 * Set new block step and update next time.
290  step(i) = ttmp
291  tnew(i) = step(i) + t0(i)
292 *
293 * See whether any KS candidates are in the same block as body #I.
294  IF (iks.GT.0.AND.i.EQ.icomp) THEN
295 * Accept same time, otherwise reduce STEP(ICOMP) and/or delay.
296  IF (t0(jcomp).EQ.t0(icomp)) THEN
297  icomp = min(icomp,jcomp)
298  jcomp = max(i,jcomp)
299  ELSE IF (t0(jcomp) + step(jcomp).LT.t0(icomp)) THEN
300  step(icomp) = 0.5d0*step(icomp)
301  tnew(icomp) = step(icomp) + t0(icomp)
302  iks = 0
303  ELSE
304  iks = 0
305  END IF
306  END IF
307 *
308 * See whether total force & derivative needs updating.
309  IF (ir.EQ.0) THEN
310 * Extrapolate regular force & first derivatives to obtain F & FDOT.
311  dtr = time - t0r(i)
312  DO 90 k = 1,3
313  f(k,i) = 0.5d0*(frdot(k,i)*dtr + fr(k,i) + firr(k))
314  fdot(k,i) = one6*(frdot(k,i) + fd(k))
315  90 CONTINUE
316  END IF
317 *
318 * Increase step counter and count perturbed c.m. steps.
319  nstepi = nstepi + 1
320  IF (i.GT.n) THEN
321  IF (list(1,2*ipair-1).GT.0) nstepb = nstepb + 1
322  END IF
323 *
324  RETURN
325 *
326  END