Nbody6
 All Files Functions Variables
regint.f
Go to the documentation of this file.
1  SUBROUTINE regint(I,XI,XIDOT)
2 *
3 *
4 * Regular 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),dv(3),fd(3),fdr(3)
11  REAL*8 frx(3),fdx(3)
12 *
13 *
14 * Set neighbour number, time-step & choice of central distance.
15  nnb0 = list(1,i)
16  dtr = time - t0r(i)
17  irskip = 0
18  IF (kz(39).EQ.0) THEN
19  ri2 = (xi(1) - rdens(1))**2 + (xi(2) - rdens(2))**2 +
20  & (xi(3) - rdens(3))**2
21  rh2 = rscale**2
22  ELSE
23  ri2 = xi(1)**2 + xi(2)**2 + xi(3)**2
24  rh2 = 9.0*rscale**2
25  END IF
26 *
27 * Obtain irregular & regular force and determine current neighbours.
28  rs2 = rs(i)**2
29 * Take volume between inner and outer radius equal to basic sphere.
30  rcrit2 = 1.59*rs2
31 * Set radial velocity factor for the outer shell.
32  vrfac = -0.1*rs2/dtr
33 * Start count at 2 and subtract 1 at the end to avoid ILIST(NNB+1).
34  nnb = 1
35 *
36 * Initialize scalars for forces & derivatives.
37  DO 5 k = 1,3
38  firr(k) = 0.0d0
39  freg(k) = 0.0d0
40  fd(k) = 0.0
41  fdr(k) = 0.0
42  5 CONTINUE
43 *
44 * Choose appropriate force loop for single particle or c.m. body.
45  IF (i.GT.n) THEN
46 * Treat unperturbed KS in the single particle approximation.
47  i1 = 2*(i - n) - 1
48  IF (list(1,i1).GT.0) THEN
49 * Obtain total force on c.m. particle.
50  CALL cmfreg(i,rs2,rcrit2,vrfac,nnb,xi,xidot,firr,freg,fd,fdr)
51  go to 20
52  END IF
53  END IF
54 *
55 * Perform fast force loop over single particles.
56  DO 10 j = ifirst,n
57 * Note that skip here may be better for optimizing compiler.
58  IF (j.EQ.i) go to 10
59  a1 = x(1,j) - xi(1)
60  a2 = x(2,j) - xi(2)
61  a3 = x(3,j) - xi(3)
62 * Predicted coordinates avoids spurious force differences.
63  dv(1) = xdot(1,j) - xidot(1)
64  dv(2) = xdot(2,j) - xidot(2)
65  dv(3) = xdot(3,j) - xidot(3)
66 *
67  rij2 = a1*a1 + a2*a2 + a3*a3
68  drdv = a1*dv(1) + a2*dv(2) + a3*dv(3)
69 *
70 * First see whether the distance exceeds the outer shell radius.
71  IF (rij2.GT.rcrit2) go to 8
72 *
73 * Retain particles with small steps (large derivative corrections).
74  IF (rij2.GT.rs2.AND.step(j).GT.10.0*smin) THEN
75 * Accept member if maximum penetration factor exceeds 8 per cent.
76  IF (drdv.GT.vrfac) go to 8
77  END IF
78 *
79 * Increase neighbour counter and obtain current irregular force.
80  nnb = nnb + 1
81  ilist(nnb) = j
82  dr2i = 1.0/rij2
83  drdv = 3.0*drdv*dr2i
84  dr3i = body(j)*dr2i*sqrt(dr2i)
85  firr(1) = firr(1) + a1*dr3i
86  firr(2) = firr(2) + a2*dr3i
87  firr(3) = firr(3) + a3*dr3i
88  fd(1) = fd(1) + (dv(1) - a1*drdv)*dr3i
89  fd(2) = fd(2) + (dv(2) - a2*drdv)*dr3i
90  fd(3) = fd(3) + (dv(3) - a3*drdv)*dr3i
91  go to 10
92 *
93 * Obtain the regular force.
94  8 dr2i = 1.0/rij2
95  drdv = 3.0*drdv*dr2i
96  dr3i = body(j)*dr2i*sqrt(dr2i)
97  freg(1) = freg(1) + a1*dr3i
98  freg(2) = freg(2) + a2*dr3i
99  freg(3) = freg(3) + a3*dr3i
100  fdr(1) = fdr(1) + (dv(1) - a1*drdv)*dr3i
101  fdr(2) = fdr(2) + (dv(2) - a2*drdv)*dr3i
102  fdr(3) = fdr(3) + (dv(3) - a3*drdv)*dr3i
103  10 CONTINUE
104 *
105 * Add any contributions from regularized c.m. particles.
106  IF (npairs.GT.0) THEN
107  CALL cmfreg(i,rs2,rcrit2,vrfac,nnb,xi,xidot,firr,freg,fd,fdr)
108  END IF
109 *
110 * Include treatment for regularized subsystem.
111  IF (nch.GT.0) THEN
112 * Distinguish between chain c.m. and any other particle.
113  IF (name(i).EQ.0) THEN
114  CALL chfirr(i,1,xi,xidot,firr,fd)
115  ELSE
116 * Search the chain perturber list for #I.
117  np1 = listc(1) + 1
118  DO 15 l = 2,np1
119  j = listc(l)
120  IF (j.GT.i) go to 20
121  IF (j.EQ.i) THEN
122  CALL fchain(i,1,xi,xidot,firr,fd)
123  go to 20
124  END IF
125  15 CONTINUE
126  END IF
127  END IF
128 *
129 * Check optional interstellar clouds.
130  20 IF (kz(13).GT.0) THEN
131  CALL fcloud(i,freg,fdr,1)
132  END IF
133 *
134 * See whether an external force should be added.
135  IF (kz(14).GT.0) THEN
136 * Save current values for deriving work done by tides (#14 = 3).
137  IF (kz(14).EQ.3) THEN
138  DO 22 k = 1,3
139  frx(k) = freg(k)
140  fdx(k) = fdr(k)
141  22 CONTINUE
142  END IF
143 *
144 * Obtain the tidal perturbation (force and first derivative).
145  CALL xtrnlf(xi,xidot,firr,freg,fd,fdr,1)
146 *
147 * Form rate of tidal energy change during last regular step.
148  IF (kz(14).EQ.3) THEN
149  wdot = 0.0
150  w2dot = 0.0
151 * W3DOT = 0.0
152 * Integrate Taylor series for V*P using final values.
153  DO 24 k = 1,3
154  px = freg(k) - frx(k)
155  dpx = fdr(k) - fdx(k)
156  wdot = wdot + xidot(k)*px
157  w2dot = w2dot + (freg(k) + firr(k))*px + xidot(k)*dpx
158 * W3DOT = W3DOT + 2.0*(FREG(K) + FIRR(K))*DPX +
159 * & (FDR(K) + FD(K))*PX
160 * Second-order term derived by Douglas Heggie (Aug/03).
161  24 CONTINUE
162 * Accumulate tidal energy change for general galactic potential.
163 * ETIDE = ETIDE - BODY(I)*((ONE6*W3DOT*DTR - 0.5*W2DOT)*DTR
164 * & + WDOT)*DTR
165 * Note Taylor series at end of interval with negative argument.
166  etide = etide + body(i)*(0.5*w2dot*dtr - wdot)*dtr
167  END IF
168  END IF
169 *
170  nnb = nnb - 1
171 * Check for zero neighbour number.
172  IF (nnb.EQ.0) THEN
173 * Assume small mass at centre to avoid zero irregular force.
174  r2 = xi(1)**2 + xi(2)**2 + xi(3)**2
175  fij = 0.01*bodym/(r2*sqrt(r2))
176  rdot = 3.0*(xi(1)*xidot(1) + xi(2)*xidot(2) +
177  & xi(3)*xidot(3))/r2
178  DO 25 k = 1,3
179  firr(k) = firr(k) - fij*xi(k)
180  fd(k) = fd(k) - (xidot(k) - rdot*xi(k))*fij
181  25 CONTINUE
182 * Modify neighbour sphere gradually but allow encounter detection.
183  IF (rdot.GT.0.0) THEN
184  rs(i) = max(0.75*rs(i),0.1*rscale)
185  ELSE
186  rs(i) = 1.1*rs(i)
187  END IF
188  nbvoid = nbvoid + 1
189  irskip = 1
190 * Skip another full N loop.
191  go to 50
192  END IF
193 *
194 * Restrict neighbour number < NNBMAX to permit one normal addition.
195  IF (nnb.LT.nnbmax) go to 40
196 *
197 * Reduce search radius by cube root of 90 % volume factor.
198  30 nnb2 = 0.9*nnbmax
199  a1 = float(nnb2)/float(nnb)
200  IF (rs(i).GT.5.0*rscale) THEN
201  a1 = min(5.0*a1,0.9d0)
202  irskip = 1
203  END IF
204  rs2 = rs2*a1**0.66667
205  rcrit2 = 1.59*rs2
206  rs(i) = sqrt(rs2)
207  nnb1 = 0
208 *
209  DO 35 l = 1,nnb
210  j = ilist(l+1)
211  IF (l + nnb2.GT.nnb + nnb1) go to 32
212 * Sum of neighbours (NNB1) & those left (NNB+1-L) set to NNB2.
213  a1 = x(1,j) - xi(1)
214  a2 = x(2,j) - xi(2)
215  a3 = x(3,j) - xi(3)
216  dv(1) = xdot(1,j) - xidot(1)
217  dv(2) = xdot(2,j) - xidot(2)
218  dv(3) = xdot(3,j) - xidot(3)
219 *
220  rij2 = a1*a1 + a2*a2 + a3*a3
221  dr2i = 1.0/rij2
222  dr3i = body(j)*dr2i*sqrt(dr2i)
223  drdv = a1*dv(1) + a2*dv(2) + a3*dv(3)
224  IF (rij2.GT.rcrit2) go to 34
225 *
226  IF (rij2.GT.rs2.AND.step(j).GT.smin) THEN
227  IF (drdv.GT.vrfac.AND.j.LE.n) go to 34
228 * Retain any c.m. because of complications in force correction.
229  END IF
230 *
231  32 nnb1 = nnb1 + 1
232  jlist(nnb1+1) = j
233  go to 35
234 *
235 * Subtract neighbour force included above and add to regular force.
236  34 drdv = 3.0*drdv*dr2i
237  firr(1) = firr(1) - a1*dr3i
238  firr(2) = firr(2) - a2*dr3i
239  firr(3) = firr(3) - a3*dr3i
240  fd(1) = fd(1) - (dv(1) - a1*drdv)*dr3i
241  fd(2) = fd(2) - (dv(2) - a2*drdv)*dr3i
242  fd(3) = fd(3) - (dv(3) - a3*drdv)*dr3i
243  freg(1) = freg(1) + a1*dr3i
244  freg(2) = freg(2) + a2*dr3i
245  freg(3) = freg(3) + a3*dr3i
246  fdr(1) = fdr(1) + (dv(1) - a1*drdv)*dr3i
247  fdr(2) = fdr(2) + (dv(2) - a2*drdv)*dr3i
248  fdr(3) = fdr(3) + (dv(3) - a3*drdv)*dr3i
249  35 CONTINUE
250 *
251  DO 38 l = 2,nnb1+1
252  ilist(l) = jlist(l)
253  38 CONTINUE
254  nnb = nnb1
255  nbfull = nbfull + 1
256 * See whether to reduce NNB further.
257  IF (nnb.GE.nnbmax) go to 30
258 *
259 * Stabilize NNB between ZNBMIN & ZNBMAX by square root of contrast.
260  40 a3 = alpha*sqrt(float(nnb)*rs(i))/rs2
261  nbp = a3
262  a3 = min(a3,znbmax)
263 * Reduce predicted membership slowly outside half-mass radius.
264  IF (ri2.GT.rh2) THEN
265  a3 = a3*rscale/sqrt(ri2)
266  ELSE
267  a3 = max(a3,znbmin)
268  END IF
269  a4 = a3/float(nnb)
270 * Include inertial factor to prevent resonance oscillations in RS.
271  IF ((a3 - float(nnb0))*(a3 - float(nnb)).LT.0.0) a4 = sqrt(a4)
272 *
273 * Modify volume ratio by radial velocity factor outside the core.
274  IF (ri2.GT.rc2.AND.kz(39).EQ.0.AND.ri2.LT.9.0*rh2) THEN
275  ridot = (xi(1) - rdens(1))*xidot(1) +
276  & (xi(2) - rdens(2))*xidot(2) +
277  & (xi(3) - rdens(3))*xidot(3)
278  a4 = a4*(1.0 + ridot*dtr/ri2)
279  END IF
280 *
281 * See whether neighbour radius of c.m. body should be increased.
282  IF (i.GT.n.AND.ri2.LT.rh2) THEN
283 * Set perturber range (soft binaries & H > 0 have short duration).
284  a2 = 100.0*body(i)/abs(h(i-n))
285 * Stabilize NNB on ZNBMAX if too few perturbers.
286  IF (a2.GT.rs(i)) THEN
287  a3 = max(1.0 - float(nnb)/znbmax,0.0d0)
288 * Modify volume ratio by approximate square root factor.
289  a4 = 1.0 + 0.5*a3
290  END IF
291  END IF
292 *
293 * Limit change of RS for small steps (RSFAC = MAX(25/TCR,1.5*VC/RSMIN).
294  a5 = min(rsfac*dtr,0.25d0)
295 * Restrict volume ratio by inertial factor of 25 per cent either way.
296  a4 = a4 - 1.0
297  IF (a4.GT.a5) THEN
298  a4 = a5
299  ELSE
300  a4 = max(a4,-a5)
301  END IF
302 *
303 * Modify neighbour sphere radius by volume factor.
304  IF (irskip.EQ.0) THEN
305  a3 = one3*a4
306  a1 = 1.0 + a3 - a3*a3
307 * Second-order cube root expansion (maximum volume error < 0.3 %).
308  IF (rs(i).GT.5.0*rscale) a1 = sqrt(a1)
309 * Prevent reduction of small NNB if predicted value exceeds ZNBMIN.
310  IF (nnb.LT.znbmin.AND.nbp.GT.znbmin) THEN
311  IF (a1.LT.1.0) a1 = 1.05
312  END IF
313 * Skip modification for small changes (avoids oscillations in RS).
314  IF (abs(a1 - 1.0d0).GT.0.003) THEN
315  rs(i) = a1*rs(i)
316  END IF
317  END IF
318 *
319 * Calculate the radial velocity with respect to at most 3 neighbours.
320  IF (nnb.LE.3.AND.ri2.GT.rh2) THEN
321  a1 = 2.0*rs(i)
322 *
323  DO 45 l = 1,nnb
324  j = ilist(l+1)
325  rij = sqrt((xi(1) - x(1,j))**2 + (xi(2) - x(2,j))**2 +
326  & (xi(3) - x(3,j))**2)
327  rsdot = ((xi(1) - x(1,j))*(xidot(1) - xdot(1,j)) +
328  & (xi(2) - x(2,j))*(xidot(2) - xdot(2,j)) +
329  & (xi(3) - x(3,j))*(xidot(3) - xdot(3,j)))/rij
330 * Find smallest neighbour distance assuming constant regular step.
331  a1 = min(a1,rij + rsdot*dtr)
332  45 CONTINUE
333 *
334 * Increase neighbour sphere if all members are leaving inner region.
335  rs(i) = max(a1,1.1*rs(i))
336 * Impose minimum size to include wide binaries (cf. NNB = 0 condition).
337  rs(i) = max(0.1*rscale,rs(i))
338  END IF
339 *
340 * Check minimum neighbour sphere since last output (skip NNB = 0).
341  IF (list(1,i).GT.0) rsmin = min(rsmin,rs(i))
342 *
343 * Check optional procedures for adding neighbours.
344  IF ((kz(37).EQ.1.AND.listv(1).GT.0).OR.kz(37).GT.1) THEN
345  CALL checkl(i,nnb,xi,xidot,rs2,firr,freg,fd,fdr)
346  END IF
347 *
348 * Find loss or gain of neighbours at the same time.
349  50 nbloss = 0
350  nbgain = 0
351 *
352 * Accumulate tidal energy change for general galactic potential.
353  IF (kz(14).EQ.3) THEN
354 * Note: Taylor series at end of interval with negative argument.
355 * ETIDE = ETIDE - BODY(I)*((ONE6*W3DOT*DTR - 0.5*W2DOT)*DTR +
356 * & WDOT)*DTR
357  etide = etide + body(i)*(0.5*w2dot*dtr - wdot)*dtr
358 * Note: integral of Taylor series for V*P using final values.
359  END IF
360 *
361 * Check case of zero old membership (NBGAIN = NNB specifies net gain).
362  IF (nnb0.EQ.0) THEN
363  nbgain = nnb
364 * Copy new members for fpcorr.f.
365  DO 52 l = 1,nnb
366  jlist(l) = ilist(l+1)
367  52 CONTINUE
368  go to 70
369  END IF
370 *
371  jmin = 0
372  l = 2
373  lg = 2
374 * Set termination value in ILIST(NNB+2) and save last list member.
375  ilist(nnb+2) = ntot + 1
376  ilist(1) = list(nnb0+1,i)
377 *
378 * Compare old and new list members in locations L & LG.
379  56 IF (list(l,i).EQ.ilist(lg)) go to 58
380 *
381 * Now check whether inequality means gain or loss.
382  IF (list(l,i).GE.ilist(lg)) THEN
383  nbgain = nbgain + 1
384  jlist(nnb0+nbgain) = ilist(lg)
385 * Number of neighbour losses can at most be NNB0.
386  l = l - 1
387 * The same location will be searched again after increasing L below.
388  ELSE
389  nbloss = nbloss + 1
390  j = list(l,i)
391  jlist(nbloss) = j
392 * Check SMIN step indicator (rare case permits fast skip below).
393  IF (step(j).LT.0.1*dtmin) jmin = j
394  lg = lg - 1
395  END IF
396 *
397 * See whether the last location has been checked.
398  58 IF (l.LE.nnb0) THEN
399  l = l + 1
400  lg = lg + 1
401 * Last value of second search index is NNB + 2 which holds NTOT + 1.
402  go to 56
403  ELSE IF (lg.LE.nnb) THEN
404  lg = lg + 1
405  list(l,i) = ntot + 1
406 * Last location of list holds termination value (saved in ILIST(1)).
407  go to 56
408  END IF
409 *
410 * See whether any old neighbour with small step should be retained.
411  IF (jmin.EQ.0) go to 70
412 *
413  k = 1
414  60 IF (nnb.GT.nnbmax.OR.i.GT.n) go to 70
415  j = jlist(k)
416 * Skip single regularized component (replaced by c.m.) and also c.m.
417  IF (step(j).GT.0.1*dtmin.OR.j.LT.ifirst.OR.j.GT.n) go to 68
418 * Retain old neighbour inside 1.4*RS to avoid large correction terms.
419  rij2 = (xi(1) - x(1,j))**2 + (xi(2) - x(2,j))**2 +
420  & (xi(3) - x(3,j))**2
421  IF (rij2.GT.2.0*rs2) go to 68
422 *
423  l = nnb + 1
424  62 IF (ilist(l).LT.j) go to 64
425  ilist(l+1) = ilist(l)
426  l = l - 1
427  IF (l.GT.1) go to 62
428 *
429 * Save index of body #J and update NNB & NBLOSS.
430  64 ilist(l+1) = j
431  nnb = nnb + 1
432  nbloss = nbloss - 1
433 * Restore last old neighbour in case NBLOSS = 0 at end of search.
434  list(nnb0+1,i) = ilist(1)
435  nbsmin = nbsmin + 1
436 *
437 * Perform correction to irregular and regular force components.
438  a1 = x(1,j) - xi(1)
439  a2 = x(2,j) - xi(2)
440  a3 = x(3,j) - xi(3)
441  dv(1) = xdot(1,j) - xidot(1)
442  dv(2) = xdot(2,j) - xidot(2)
443  dv(3) = xdot(3,j) - xidot(3)
444 *
445  rij2 = a1*a1 + a2*a2 + a3*a3
446  dr2i = 1.0/rij2
447  dr3i = body(j)*dr2i*sqrt(dr2i)
448  drdv = 3.0*(a1*dv(1) + a2*dv(2) + a3*dv(3))*dr2i
449 *
450  firr(1) = firr(1) + a1*dr3i
451  firr(2) = firr(2) + a2*dr3i
452  firr(3) = firr(3) + a3*dr3i
453  fd(1) = fd(1) + (dv(1) - a1*drdv)*dr3i
454  fd(2) = fd(2) + (dv(2) - a2*drdv)*dr3i
455  fd(3) = fd(3) + (dv(3) - a3*drdv)*dr3i
456  freg(1) = freg(1) - a1*dr3i
457  freg(2) = freg(2) - a2*dr3i
458  freg(3) = freg(3) - a3*dr3i
459  fdr(1) = fdr(1) - (dv(1) - a1*drdv)*dr3i
460  fdr(2) = fdr(2) - (dv(2) - a2*drdv)*dr3i
461  fdr(3) = fdr(3) - (dv(3) - a3*drdv)*dr3i
462 *
463 * Remove body #J from JLIST unless it is the last or only member.
464  IF (k.GT.nbloss) go to 70
465  DO 66 l = k,nbloss
466  jlist(l) = jlist(l+1)
467  66 CONTINUE
468 * Index of last body to be moved up is L = NBLOSS + 1.
469  k = k - 1
470 * Check the same location again since a new body has appeared.
471  68 k = k + 1
472 * Last member to be considered is in location K = NBLOSS.
473  IF (k.LE.nbloss) go to 60
474 *
475 * Form time-step factors and update regular force time.
476  70 dtsq = dtr**2
477  dt6 = 6.0d0/(dtr*dtsq)
478  dt2 = 2.0d0/dtsq
479  dtsq12 = one12*dtsq
480  dtr13 = one3*dtr
481  t0r(i) = time
482 *
483 * Suppress the corrector for large time-step ratios (experimental).
484  IF (step(i).LT.5.0d0*dtmin.AND.dtr.GT.50.0*step(i)) THEN
485  dtr13 = 0.0
486  dtsq12 = 0.0
487  END IF
488 *
489 * Include the corrector and save F, FI, FR & polynomial derivatives.
490  DO 75 k = 1,3
491 * Subtract change of neighbour force to form actual first derivative.
492  dfr = fr(k,i) - (firr(k) - fi(k,i)) - freg(k)
493  fdr0 = fdr(k) - (fidot(k,i) - fd(k))
494 *
495  frd = frdot(k,i)
496  sum = frd + fdr0
497  at3 = 2.0d0*dfr + dtr*sum
498  bt2 = -3.0d0*dfr - dtr*(sum + frd)
499 *
500  x0(k,i) = x0(k,i) + (0.6d0*at3 + bt2)*dtsq12
501  x0dot(k,i) = x0dot(k,i) + (0.75d0*at3 + bt2)*dtr13
502 *
503 * X0(K,I) = X(K,I)
504 * X0DOT(K,I) = XDOT(K,I)
505 *
506  fi(k,i) = firr(k)
507  fr(k,i) = freg(k)
508  f(k,i) = 0.5d0*(freg(k) + firr(k))
509  fidot(k,i) = fd(k)
510  frdot(k,i) = fdr(k)
511  fdot(k,i) = one6*(fdr(k) + fd(k))
512 *
513  d0(k,i) = firr(k)
514  d0r(k,i) = freg(k)
515 * D0R(K,I) = FREG(K) - (FI(K,I) - FIRR(K))
516 * D1R(K,I) = FDR0
517 * Use actual first derivatives (2nd derivs only consistent in FPCORR).
518  d1(k,i) = fd(k)
519  d1r(k,i) = fdr(k)
520 * Set second & third derivatives based on old neighbours (cf. FPCORR).
521  d2r(k,i) = (3.0d0*at3 + bt2)*dt2
522  d3r(k,i) = at3*dt6
523  75 CONTINUE
524 *
525 * Correct force polynomials due to neighbour changes (KZ(38) or I > N).
526  IF (kz(38).GT.0) THEN
527  CALL fpcorr(i,nbloss,nbgain,xi,xidot)
528  ELSE IF (i.GT.n) THEN
529  IF (gamma(i-n).GT.gmax) THEN
530  CALL fpcorr(i,nbloss,nbgain,xi,xidot)
531  END IF
532  END IF
533 *
534 * Copy new neighbour list if different from old list.
535  IF (nbloss + nbgain.GT.0) THEN
536  list(1,i) = nnb
537  DO 80 l = 2,nnb+1
538  list(l,i) = ilist(l)
539  80 CONTINUE
540  nbflux = nbflux + nbloss + nbgain
541  END IF
542 *
543 * Check for boundary reflection (RI2 < 0 denotes new polynomials set).
544 * IF (KZ(29).GT.0) THEN
545 * RI2 = X(1,I)**2 + X(2,I)**2 + X(3,I)**2
546 * IF (RI2.GT.RSPH2) THEN
547 * CALL REFLCT(I,RI2)
548 * IF (RI2.LT.0.0) GO TO 120
549 * END IF
550 * END IF
551 *
552 * Obtain new regular integration step using composite expression.
553 * STEPR = (ETAR*(F*F2DOT + FDOT**2)/(FDOT*F3DOT + F2DOT**2))**0.5.
554  fr2 = freg(1)**2 + freg(2)**2 + freg(3)**2
555  w1 = fdr(1)**2 + fdr(2)**2 + fdr(3)**2
556  w2 = d2r(1,i)**2 + d2r(2,i)**2 + d2r(3,i)**2
557  w3 = d3r(1,i)**2 + d3r(2,i)**2 + d3r(3,i)**2
558 *
559 * Form new step by relative criterion.
560  w0 = (sqrt(fr2*w2) + w1)/(sqrt(w1*w3) + w2)
561  w0 = etar*w0
562  ttmp = sqrt(w0)
563  dt0 = ttmp
564 *
565 * Determine new regular step.
566 * TTMP = TSTEP(FREG,FDR,D2R(1,I),D3R(1,I),ETAR)
567 *
568 * Adopt FAC*MIN(FREG,FIRR) (or tidal force) for convergence test.
569  fac = etar
570  IF (kz(14).EQ.0) THEN
571  fi2 = firr(1)**2 + firr(2)**2 + firr(3)**2
572  df2 = fac**2*min(fr2,fi2)
573 * Ignore irregular force criterion if no neighbours.
574  IF (nnb.EQ.0) df2 = fac**2*fr2
575  ELSE IF (kz(14).EQ.1) THEN
576  w0 = (tidal(1)*xi(1))**2
577  df2 = fac**2*max(fr2,w0)
578  ELSE
579  df2 = fac**2*fr2
580  END IF
581 *
582 * Obtain regular force change using twice the predicted step.
583  dtc = 2.0*ttmp
584  s2 = 0.5*dtc
585  s3 = one3*dtc
586  w0 = 0.0
587  DO 105 k = 1,3
588  w2 = ((d3r(k,i)*s3 + d2r(k,i))*s2 + fdr(k))*dtc
589  w0 = w0 + w2**2
590  105 CONTINUE
591 *
592 * See whether regular step can be increased by factor 2.
593  IF (w0.LT.df2) THEN
594  ttmp = dtc
595  END IF
596 *
597 * Impose a smooth step reduction inside compact core.
598  IF (nc.LT.50.AND.ri2.LT.rc2) THEN
599  ttmp = ttmp*min(1.0d0,0.5d0*(1.0d0 + ri2*rc2in))
600  END IF
601 *
602 * Select discrete value (increased by 2, decreased by 2 or unchanged).
603  IF (ttmp.GT.2.0*stepr(i)) THEN
604  IF (dmod(time,2.0*stepr(i)).EQ.0.0d0) THEN
605  ttmp = min(2.0*stepr(i),smax)
606  ELSE
607  ttmp = stepr(i)
608  END IF
609  ELSE IF (ttmp.LT.stepr(i)) THEN
610  ttmp = 0.5*stepr(i)
611 * Allow a second reduction to prevent spurious contributions.
612  IF (ttmp.GT.dt0) THEN
613  ttmp = 0.5*ttmp
614  END IF
615  ELSE
616  ttmp = stepr(i)
617  END IF
618 *
619 * Set new regular step and reduce irregular step if STEPR < STEP.
620  stepr(i) = ttmp
621  110 IF (ttmp.LT.step(i)) THEN
622  step(i) = 0.5d0*step(i)
623  tnew(i) = t0(i) + step(i)
624  niconv = niconv + 1
625  go to 110
626  END IF
627 *
628 * Reduce irregular step on switching from zero neighbour number.
629  IF (nnb0.EQ.0.AND.nnb.GT.0) THEN
630  step(i) = 0.25*step(i)
631  tnew(i) = t0(i) + step(i)
632  END IF
633 *
634  nstepr = nstepr + 1
635 *
636  RETURN
637 *
638  END