Nbody6
xvpred.f
Go to the documentation of this file.
1  SUBROUTINE xvpred(I1,NNB)
2 *
3 *
4 * Prediction of coordinates & velocities.
5 * ---------------------------------------
6 *
7  include 'common6.h'
8 *
9 *
10 * Perform full N prediction to highest order (defined by NNB >= N).
11  IF (nnb.GE.n) go to 20
12 * Use high order for special singles (NNB < 0) or c.m. bodies (I1 > N).
13  IF (nnb.LT.0.OR.(nnb.EQ.0.AND.i1.GT.n)) go to 20
14 *
15 * Adopt low order for standard single particles (I1 <= N & NNB = 0).
16  IF (nnb.EQ.0) THEN
17  i = i1
18  s = time - t0(i)
19  IF (s.EQ.0.0d0) go to 40
20  s1 = 1.5*s
21  s2 = 2.0*s
22  x(1,i) = ((fdot(1,i)*s + f(1,i))*s + x0dot(1,i))*s + x0(1,i)
23  x(2,i) = ((fdot(2,i)*s + f(2,i))*s + x0dot(2,i))*s + x0(2,i)
24  x(3,i) = ((fdot(3,i)*s + f(3,i))*s + x0dot(3,i))*s + x0(3,i)
25  xdot(1,i) = (fdot(1,i)*s1 + f(1,i))*s2 + x0dot(1,i)
26  xdot(2,i) = (fdot(2,i)*s1 + f(2,i))*s2 + x0dot(2,i)
27  xdot(3,i) = (fdot(3,i)*s1 + f(3,i))*s2 + x0dot(3,i)
28  go to 40
29  END IF
30 *
31 * Predict all neighbours to order FDOT.
32  nnb1 = nnb + 1
33  DO 5 l = 2,nnb1
34  i = list(l,i1)
35  s = time - t0(i)
36  IF (s.EQ.0.0d0) go to 5
37  s1 = 1.5*s
38  s2 = 2.0*s
39  x(1,i) = ((fdot(1,i)*s + f(1,i))*s + x0dot(1,i))*s + x0(1,i)
40  x(2,i) = ((fdot(2,i)*s + f(2,i))*s + x0dot(2,i))*s + x0(2,i)
41  x(3,i) = ((fdot(3,i)*s + f(3,i))*s + x0dot(3,i))*s + x0(3,i)
42  xdot(1,i) = (fdot(1,i)*s1 + f(1,i))*s2 + x0dot(1,i)
43  xdot(2,i) = (fdot(2,i)*s1 + f(2,i))*s2 + x0dot(2,i)
44  xdot(3,i) = (fdot(3,i)*s1 + f(3,i))*s2 + x0dot(3,i)
45  5 CONTINUE
46 *
47 * Resolve the components of any perturbed pairs (last neighbours).
48  10 i = list(nnb1,i1)
49  IF (i.GT.n) THEN
50  jpair = i - n
51  IF (list(1,2*jpair-1).GT.0) THEN
52  CALL resolv(jpair,1)
53  END IF
54  nnb1 = nnb1 - 1
55  IF (nnb1.GT.1) go to 10
56  END IF
57  go to 40
58 *
59 * Set index of first body (#I1 if NNB = 0 or NNB >= N).
60  20 i = i1
61 * Predict coordinates & velocities of body #I to order F3DOT.
62  25 dt = time - t0(i)
63  IF ((dt.EQ.0.0d0.AND.iphase.LT.4).OR.body(i).EQ.0.0d0) go to 35
64  a1 = 0.05*dt
65  a2 = 0.25*dt
66  a3 = t0(i) - t0r(i)
67  a4 = 0.5*a3
68 *
69  DO 30 k = 1,3
70  fk = ((one6*d3r(k,i)*a3 + 0.5*d2r(k,i))*a3 + d1r(k,i))*a3
71  & + d0r(k,i) + d0(k,i)
72  f1dotk = (d3r(k,i)*a4 + d2r(k,i))*a3 + d1r(k,i) + d1(k,i)
73  f2dotk = 0.5*(d3r(k,i)*a3 + d2r(k,i) + d2(k,i))
74  f3dotk = one6*(d3r(k,i) + d3(k,i))
75  x(k,i) = ((((f3dotk*a1 + one12*f2dotk)*dt +
76  & one6*f1dotk)*dt + 0.5*fk)*dt + x0dot(k,i))*dt +
77  & x0(k,i)
78  xdot(k,i) = (((f3dotk*a2 + one3*f2dotk)*dt +
79  & 0.5*f1dotk)*dt + fk)*dt + x0dot(k,i)
80  30 CONTINUE
81 *
82 * Resolve the components of perturbed pairs.
83  35 IF (i.GT.n) THEN
84  jpair = i - n
85  IF (list(1,2*jpair-1).GT.0) THEN
86  CALL resolv(jpair,1)
87  END IF
88  END IF
89 *
90 * Check whether prediction of single particle or full N.
91  IF (nnb.GT.0) THEN
92  i = i + 1
93  IF (i.LE.nnb) go to 25
94  END IF
95 *
96  40 RETURN
97 *
98  END