Nbody6
ficorr.f
Go to the documentation of this file.
1  SUBROUTINE ficorr(I,DM)
2 *
3 *
4 * Local force corrections due to mass loss.
5 * -----------------------------------------
6 *
7  include 'common6.h'
8  REAL*8 a(6)
9 *
10 *
11 * Correct force & first derivative for neighbours of body #I.
12  nnb1 = list(1,i) + 1
13  DO 20 l = 2,nnb1
14  j = list(l,i)
15  nnb2 = list(1,j) + 1
16 *
17 * Ensure that body #I is a neighbour of body #J.
18  DO 1 k = 2,nnb2
19  IF (list(k,j).EQ.i) go to 5
20  1 CONTINUE
21  go to 20
22 *
23  5 rij2 = 0.0d0
24  a7 = 0.0d0
25 *
26  DO 10 k = 1,3
27  a(k) = x(k,i) - x(k,j)
28  a(k+3) = xdot(k,i) - xdot(k,j)
29  rij2 = rij2 + a(k)**2
30  a7 = a7 + a(k)*a(k+3)
31  10 CONTINUE
32 *
33  a7 = 3.0*a7/rij2
34  a8 = dm/(rij2*sqrt(rij2))
35 *
36  DO 15 k = 1,3
37  a(k+3) = (a(k+3) - a7*a(k))*a8
38  f(k,j) = f(k,j) - 0.5*a(k)*a8
39  fi(k,j) = fi(k,j) - a(k)*a8
40  fdot(k,j) = fdot(k,j) - one6*a(k+3)
41  d1(k,j) = d1(k,j) - a(k+3)
42  fidot(k,j) = fidot(k,j) - a(k+3)
43  15 CONTINUE
44  20 CONTINUE
45 *
46 * Obtain the potential energy due to all particles.
47  potj = 0.0d0
48  DO 30 j = ifirst,ntot
49  IF (j.EQ.i) go to 30
50  rij2 = 0.0d0
51  DO 25 k = 1,3
52  rij2 = rij2 + (x(k,i) - x(k,j))**2
53  25 CONTINUE
54  potj = potj + body(j)/sqrt(rij2)
55  30 CONTINUE
56 *
57 * Update the potential energy loss.
58  emdot = emdot - dm*potj
59 *
60 * Include kinetic energy correction.
61  vi2 = xdot(1,i)**2 + xdot(2,i)**2 + xdot(3,i)**2
62  emdot = emdot + 0.5*dm*vi2
63 *
64 * See whether linearized tidal terms should be included.
65  IF (kz(14).GT.0.AND.kz(14).LT.3) THEN
66  ecdot = ecdot - 0.5*dm*(tidal(1)*x(1,i)**2 +
67  & tidal(3)*x(3,i)**2)
68  END IF
69 *
70 * Check optional Plummer potential.
71  IF (kz(14).EQ.4.OR.kz(14).EQ.3) THEN
72  ri2 = ap2
73  DO 50 k = 1,3
74  ri2 = ri2 + x(k,i)**2
75  50 CONTINUE
76  ecdot = ecdot - dm*mp/sqrt(ri2)
77  END IF
78 *
79  RETURN
80 *
81  END