Nbody6
 All Files Functions Variables
fpcorr.f
Go to the documentation of this file.
1  SUBROUTINE fpcorr(I,NBLOSS,NBGAIN,XI,XIDOT)
2 *
3 *
4 * Force polynomial corrections.
5 * -----------------------------
6 *
7  include 'common6.h'
8  REAL*8 xi(3),xidot(3),save1(3),save2(3),save3(3),
9  & a(12),f2dot(3),f3dot(4)
10 *
11 *
12 * See whether there has been a change of neighbours.
13  nbsum = nbloss + nbgain
14  IF (nbsum.EQ.0) go to 100
15 *
16 * Initialize the derivative corrections.
17  DO 70 k = 1,3
18  save1(k) = 0.0
19  save2(k) = 0.0
20  save3(k) = 0.0
21  70 CONTINUE
22 *
23 * Form compact list of NBLOSS & NBGAIN.
24  nnb0 = list(1,i)
25  IF (nbgain.GT.0) THEN
26  DO 78 l = 1,nbgain
27  jlist(nbloss+l) = jlist(nnb0+l)
28  78 CONTINUE
29  END IF
30 *
31 * Accumulate derivative corrections.
32  l = 1
33  80 j = jlist(l)
34 *
35 * Use c.m. values of XDOT, F & FDOT for single KS components.
36  IF (j.LT.ifirst) THEN
37  jcm = n + kvec(j)
38 * STEPJ = STEP(JCM)
39  s = time - t0(jcm)
40 *
41  DO 82 k = 1,3
42  a(k) = x(k,j) - xi(k)
43  a(k+3) = xdot(k,j) - xidot(k)
44  a(k+6) = 2.0*(f(k,jcm) - f(k,i))
45  a(k+9) = 6.0*(fdot(k,jcm) - fdot(k,i))
46  82 CONTINUE
47 *
48  ncmder = ncmder + 1
49  ELSE
50 * STEPJ = STEP(J)
51  s = time - t0(j)
52  s3 = 3.0*s
53 *
54 * Predict F & FDOT of body #J to order FDOT.
55  DO 84 k = 1,3
56  a(k) = x(k,j) - xi(k)
57  a(k+3) = (fdot(k,j)*s3 + 2.0*f(k,j))*s + x0dot(k,j) -
58  & xidot(k)
59  a(k+6) = 2.0*(fdot(k,j)*s3 + f(k,j) - f(k,i))
60  a(k+9) = 6.0*(fdot(k,j) - fdot(k,i))
61  84 CONTINUE
62  END IF
63 *
64  a13 = 1.0/(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
65  a14 = body(j)*a13*sqrt(a13)
66  a15 = (a(1)*a(4) + a(2)*a(5) + a(3)*a(6))*a13
67  a16 = a15*a15
68  a17 = 3.0*a15
69  a18 = 6.0*a15
70  a19 = 9.0*a15
71  a20 = (a(4)*a(4) + a(5)*a(5) + a(6)*a(6) + a(1)*a(7) + a(2)*a(8)
72  & + a(3)*a(9))*a13 + a16
73  a21 = 9.0*a20
74  a20 = 3.0*a20
75  a22 = (9.0*(a(4)*a(7) + a(5)*a(8) + a(6)*a(9)) + 3.0*(a(1)*a(10)
76  & + a(2)*a(11) + a(3)*a(12)))*a13 + a17*(a20 - 4.0*a16)
77 *
78  DO 85 k = 1,3
79  f1dotk = a(k+3) - a17*a(k)
80  f2dot(k) = (a(k+6) - a18*f1dotk - a20*a(k))*a14
81  f3dot(k) = (a(k+9) - a21*f1dotk - a22*a(k))*a14 - a19*f2dot(k)
82 * Suppress F1DOT terms (already done in REGINT).
83 * F1DOT(K) = F1DOTK*A14
84  85 CONTINUE
85 *
86 * Check the third derivative in case of small neighbour steps.
87 * IF (STEPJ.LT.10.0*SMIN) THEN
88 * F3DOT(4) = ABS(F3DOT(1)) + ABS(F3DOT(2)) + ABS(F3DOT(3))
89 * A3 = ABS(D3R(1,I)) + ABS(D3R(2,I)) + ABS(D3R(3,I))
90 * IF (F3DOT(4).GT.10.0*A3) THEN
91 * Ignore large third derivative terms to avoid convergence problems.
92 * DO 86 K = 1,3
93 * F3DOT(K) = 0.0
94 * 86 CONTINUE
95 * NBDER = NBDER + 1
96 * END IF
97 * END IF
98 *
99 * Change the sign for NBLOSS contributions.
100  IF (l.LE.nbloss) THEN
101  DO 88 k = 1,3
102 * F1DOT(K) = -F1DOT(K)
103  f2dot(k) = -f2dot(k)
104  f3dot(k) = -f3dot(k)
105  88 CONTINUE
106  END IF
107 *
108 * Include derivative corrections from losses & gains.
109  DO 90 k = 1,3
110 * SAVE1(K) = SAVE1(K) + F1DOT(K)
111  save2(k) = save2(k) + f2dot(k)
112  save3(k) = save3(k) + f3dot(k)
113  90 CONTINUE
114 *
115  l = l + 1
116  IF (l.LE.nbsum) go to 80
117 *
118  nbcorr = nbcorr + 1
119 *
120 * Perform corrections to irregular and regular force derivatives.
121  DO 95 k = 1,3
122 * Note that corrected value of D1 & D1R already set in routine REGINT.
123 * D1(K,I) = D1(K,I) + SAVE1(K)
124  d2(k,i) = d2(k,i) + save2(k)
125  d3(k,i) = d3(k,i) + save3(k)
126 * D1R(K,I) = D1R(K,I) - SAVE1(K)
127  d2r(k,i) = d2r(k,i) - save2(k)
128  d3r(k,i) = d3r(k,i) - save3(k)
129  95 CONTINUE
130 *
131  100 RETURN
132 *
133  END