Nbody6
ntint.f
Go to the documentation of this file.
1  SUBROUTINE ntint(I)
2 *
3 *
4 * Single star integration.
5 * ------------------------
6 *
7  include 'common6.h'
8  REAL*8 xi(3),xidot(3),firr(3),fd(3)
9 *
10 *
11 * Predict position and velocity up to end of STEP.
12  iter = 0
13  1 j = i
14  s = step(j)
15  s1 = 1.5*s
16  s2 = 2.0*s
17  x(1,j) = ((fdot(1,j)*s + f(1,j))*s +x0dot(1,j))*s + x0(1,j)
18  x(2,j) = ((fdot(2,j)*s + f(2,j))*s +x0dot(2,j))*s + x0(2,j)
19  x(3,j) = ((fdot(3,j)*s + f(3,j))*s +x0dot(3,j))*s + x0(3,j)
20  xdot(1,j) = (fdot(1,j)*s1 + f(1,j))*s2 + x0dot(1,j)
21  xdot(2,j) = (fdot(2,j)*s1 + f(2,j))*s2 + x0dot(2,j)
22  xdot(3,j) = (fdot(3,j)*s1 + f(3,j))*s2 + x0dot(3,j)
23 *
24 * Copy X & XDOT to scalars and initialize FIRR & FD.
25  DO 5 k = 1,3
26  xi(k) = x(k,i)
27  xidot(k) = xdot(k,i)
28  firr(k) = 0.0d0
29  fd(k) = 0.0d0
30  5 CONTINUE
31 *
32 * Evaluate the galactic force and first derivative.
33  CALL xtrnlt(xi,xidot,firr,fd)
34 *
35 * Include the corrector and set new T0, F, FDOT, D1, D2 & D3.
36  dt = step(i)
37  dtsq = dt**2
38  dt6 = 6.0d0/(dt*dtsq)
39  dt2 = 2.0d0/dtsq
40  dtsq12 = one12*dtsq
41  dt13 = one3*dt
42  t0(i) = t0(i) + step(i)
43 *
44  DO 10 k = 1,3
45  df = fi(k,i) - firr(k)
46  fid = fidot(k,i)
47  sum = fid + fd(k)
48  at3 = 2.0d0*df + dt*sum
49  bt2 = -3.0d0*df - dt*(sum + fid)
50 *
51  x0(k,i) = xi(k) + (0.6d0*at3 + bt2)*dtsq12
52  x0dot(k,i) = xidot(k) + (0.75d0*at3 + bt2)*dt13
53 *
54 * Update the corrected values (OK for test particles).
55  x(k,i) = x0(k,i)
56  xdot(k,i) = x0dot(k,i)
57 *
58  fi(k,i) = firr(k)
59  fidot(k,i) = fd(k)
60  f(k,i) = 0.5d0*firr(k)
61  fdot(k,i) = one6*fd(k)
62 *
63 * Form derivatives even though not needed for commensurate times.
64  d1(k,i) = fd(k)
65  d2(k,i) = (3.0d0*at3 + bt2)*dt2
66  d3(k,i) = at3*dt6
67 * NOTE: These are real derivatives!
68  10 CONTINUE
69 *
70 * Specify new time-step by standard criterion.
71  ttmp = tstep(firr,fd,d2(1,i),d3(1,i),etai)
72  dt0 = ttmp
73  ttime = t0(i)
74 *
75 * Select discrete value (increased by 2, decreased by 2 or unchanged).
76  IF (ttmp.GT.2.0*step(i)) THEN
77  ttmp = min(2.0*step(i),1.0d0)
78  ELSE IF (ttmp.LT.step(i)) THEN
79  ttmp = 0.5*step(i)
80  ELSE
81  ttmp = step(i)
82  END IF
83 *
84 * Do not permit current TIME to be exceeded.
85  IF (t0(i) + ttmp.GT.time) THEN
86  ttmp = time - t0(i)
87  END IF
88 *
89 * Set new block step and update next time.
90  step(i) = ttmp
91  tnew(i) = step(i) + t0(i)
92 *
93 * Increase step counter.
94  nstail = nstail + 1
95 *
96 * See whether to continue until end of large block-step.
97  IF (tnew(i).LT.time) THEN
98  iter = iter + 1
99  IF (iter.LT.10) go to 1
100  WRITE (6,20) i, time, step(i), firr
101  20 FORMAT (' SMALL TIDAL STEP I T DT F ',i7,f7.1,1p,4e10.2)
102  END IF
103 *
104  RETURN
105 *
106  END