Nbody6
 All Files Functions Variables
fpoly2.f
Go to the documentation of this file.
1  SUBROUTINE fpoly2(I1,I2,KCASE)
2 *
3 *
4 * Second & third force derivative.
5 * --------------------------------
6 *
7  include 'common6.h'
8  REAL*8 a(12),f2dot(3),f3dot(3)
9 *
10 *
11 * Include standard case, new c.m. or KS termination (KCASE = 0, 1, 2).
12  jlast = ntot
13 * Reduce loop size for new c.m. polynomial.
14  IF (kcase.EQ.1) jlast = ntot - 1
15 *
16 * Loop over all bodies, pair #ICOMP & JCOMP, c.m. or one single body.
17  DO 70 i = i1,i2
18 *
19 * Initialize the higher differences for body #I.
20  DO 10 k = 1,3
21  d2(k,i) = 0.0d0
22  d3(k,i) = 0.0d0
23  d2r(k,i) = 0.0d0
24  d3r(k,i) = 0.0d0
25  10 CONTINUE
26 *
27  nnb = list(1,i)
28 * Neglect F2DOT & F3DOT outside 5*RS unless high accuracy is needed.
29  rcrit2 = 25.0*rs(i)**2*(1.0 + 1.0/dfloat(nnb+1))
30 * Specify index of first neighbour to be identified.
31  l = 2
32  namej = list(l,i)
33 *
34 * Sum over c.m. instead of regularized neighbours since F not known.
35  DO 60 j = ifirst,jlast
36 * Note IFIRST = 2*NPAIRS + 1, JCOMP + 1, ICOMP for KCASE = 0, 1, 2.
37  DO 15 k = 1,3
38  a(k) = x(k,j) - x(k,i)
39  15 CONTINUE
40  rij2 = a(1)*a(1) + a(2)*a(2) + a(3)*a(3)
41 *
42 * Ensure that all neighbours are considered.
43  IF (rij2.GT.rcrit2.AND.j.NE.namej) go to 60
44 * Distant bodies do not contribute significantly to F2DOT & F3DOT.
45 *
46  IF (kcase.GT.0) THEN
47  IF (j.GT.jcomp.OR.kcase.EQ.1) go to 30
48  END IF
49  IF (j.EQ.i) go to 60
50 *
51 * See whether F & FDOT extrapolation is required.
52  IF (kcase.EQ.0.AND.time.GT.0.0d0) THEN
53  IF (t0(j).LT.time) go to 30
54 * Note that routine FCLOSE sets T0(J) = TIME for dominant bodies.
55  END IF
56 *
57 * Copy F & FDOT (all J at TIME = 0, otherwise dominant bodies only).
58  DO 20 k = 1,3
59  a(k+6) = f(k,j)
60  a(k+9) = fdot(k,j)
61  20 CONTINUE
62  go to 40
63 *
64 * Obtain current force and first derivative to second order.
65  30 dt = time - t0(j)
66  dt1 = 0.5*dt
67  dtr = time - t0r(j)
68  dt1r = 0.5*dtr
69 *
70  DO 35 k = 1,3
71  a(k+6) = (d2r(k,j)*dt1r + d1r(k,j))*dtr + fr(k,j) +
72  & (d2(k,j)*dt1 + d1(k,j))*dt + fi(k,j)
73  a(k+9) = d2r(k,j)*dtr + d1r(k,j) + d2(k,j)*dt + d1(k,j)
74  35 CONTINUE
75 *
76  40 DO 45 k = 1,3
77  a(k+3) = xdot(k,j) - xdot(k,i)
78  a(k+6) = a(k+6) - f(k,i)
79  a(k+9) = a(k+9) - fdot(k,i)
80  45 CONTINUE
81 *
82  a13 = 1.0/rij2
83  a14 = body(j)*a13*sqrt(a13)
84  a15 = (a(1)*a(4) + a(2)*a(5) + a(3)*a(6))*a13
85  a16 = a15*a15
86  a17 = 3.0*a15
87  a18 = 6.0*a15
88  a19 = 9.0*a15
89  a20 = (a(4)*a(4) + a(5)*a(5) + a(6)*a(6) + a(1)*a(7) +
90  & a(2)*a(8) + a(3)*a(9))*a13 + a16
91  a21 = 9.0*a20
92  a20 = 3.0*a20
93  a22 = (9.0*(a(4)*a(7) + a(5)*a(8) + a(6)*a(9)) +
94  & 3.0*(a(1)*a(10) + a(2)*a(11) + a(3)*a(12)))*a13 +
95  & a17*(a20 - 4.0*a16)
96 *
97  DO 50 k = 1,3
98  f1dotk = a(k+3) - a17*a(k)
99  f2dot(k) = (a(k+6) - a18*f1dotk - a20*a(k))*a14
100  f3dot(k) = (a(k+9) - a21*f1dotk - a22*a(k))*a14 -
101  & a19*f2dot(k)
102  50 CONTINUE
103 *
104 * See whether body #J is a neighbour of body #I.
105  IF (j.NE.namej) THEN
106  DO 52 k = 1,3
107  d2r(k,i) = d2r(k,i) + f2dot(k)
108  d3r(k,i) = d3r(k,i) + f3dot(k)
109  52 CONTINUE
110  ELSE
111  DO 55 k = 1,3
112  d2(k,i) = d2(k,i) + f2dot(k)
113  d3(k,i) = d3(k,i) + f3dot(k)
114  55 CONTINUE
115 *
116 * Advance the neighbour list until last member is identified.
117  IF (l.LE.nnb) THEN
118  l = l + 1
119  namej = list(l,i)
120  END IF
121  END IF
122  60 CONTINUE
123  70 CONTINUE
124 *
125 * Check option for external force.
126  IF (kz(14).GT.0) THEN
127  CALL xtrnld(i1,i2,2)
128  END IF
129 *
130 * Set new time-steps and initialize prediction variables.
131  CALL steps(i1,i2)
132 *
133  RETURN
134 *
135  END