Nbody6
cmcorr.f
Go to the documentation of this file.
1  SUBROUTINE cmcorr
2 *
3 *
4 * Center of mass & total force corrections.
5 * -----------------------------------------
6 *
7  include 'common6.h'
8 *
9 *
10 * Initialize centre of mass variables.
11  DO 10 k = 1,3
12  cmr(k) = 0.0d0
13  cmrdot(k) = 0.0d0
14  10 CONTINUE
15 *
16 * Form c.m. coordinate & velocity displacements.
17  DO 20 i = ifirst,ntot
18  DO 15 k = 1,3
19  cmr(k) = cmr(k) + body(i)*x(k,i)
20  cmrdot(k) = cmrdot(k) + body(i)*xdot(k,i)
21  15 CONTINUE
22  20 CONTINUE
23 *
24  DO 30 k = 1,3
25  cmr(k) = cmr(k)/zmass
26  cmrdot(k) = cmrdot(k)/zmass
27  30 CONTINUE
28 *
29 * Include effect of c.m. motion in Plummer potential.
30  IF (kz(14).EQ.4) THEN
31  CALL plpot1(phi1)
32  END IF
33 *
34 * Apply c.m. corrections to X & XDOT and accumulate energy changes.
35  errx = 0.0d0
36  errv = 0.0d0
37  DO 40 i = ifirst,ntot
38  DO 35 k = 1,3
39  xi2 = x(k,i)**2
40  vi2 = xdot(k,i)**2
41  x(k,i) = x(k,i) - cmr(k)
42  xdot(k,i) = xdot(k,i) - cmrdot(k)
43 * Note TIDAL(K) = 0 for KZ(14) = 3, hence no skip on ERRX is needed.
44  errx = errx - tidal(k)*body(i)*(x(k,i)**2 - xi2)
45  errv = errv + body(i)*(xdot(k,i)**2 - vi2)
46  35 CONTINUE
47  40 CONTINUE
48 *
49 * Adjust the total energy to new kinetic energy & tidal potential.
50  be(3) = be(3) + 0.5*(errx + errv)
51  e(11) = e(11) - 0.5*(errx + errv)
52 *
53 * Perform a consistent shift of the density centre.
54  DO 50 k = 1,3
55  rdens(k) = rdens(k) - cmr(k)
56  50 CONTINUE
57 *
58 * Subtract tidal corrections from total force & first derivative.
59  IF (kz(14).GT.0.AND.kz(14).LE.2) THEN
60  DO 60 i = ifirst,ntot
61 * Skip ghosts to avoid spurious prediction inside 1.0E+10.
62  IF (body(i).EQ.0.0d0) go to 60
63  DO 55 k = 1,3
64  df = tidal(k)*cmr(k)
65  dd = tidal(k)*cmrdot(k)
66  fr(k,i) = fr(k,i) - df
67  f(k,i) = f(k,i) - 0.5*df
68  d1r(k,i) = d1r(k,i) - dd
69  frdot(k,i) = frdot(k,i) - dd
70  fdot(k,i) = fdot(k,i) - one6*dd
71  IF (k.EQ.1) THEN
72  fi(1,i) = fi(1,i) - tidal(4)*cmrdot(2)
73  ELSE IF (k.EQ.2) THEN
74  fi(2,i) = fi(2,i) + tidal(4)*cmrdot(1)
75  END IF
76  55 CONTINUE
77  60 CONTINUE
78  END IF
79 *
80 * Re-determine X0 & X0DOT consistently with current corrected X & XDOT.
81  DO 70 i = ifirst,ntot
82  IF (body(i).EQ.0.0d0) go to 70
83  dt = time - t0(i)
84  dtr = time - t0r(i)
85  a1 = 0.2*dt
86  a2 = dt/24.0
87 *
88  DO 65 k = 1,3
89  f2dotk = d3r(k,i)*dtr + d2r(k,i) + d3(k,i)*dt + d2(k,i)
90  f3dotk = d3r(k,i) + d3(k,i)
91  dv0 = (((f3dotk*a2 + one6*f2dotk)*dt +
92  & 3.0*fdot(k,i))*dt + 2.0*f(k,i))*dt
93  x0dot(k,i) = xdot(k,i) - dv0
94  dx0 = ((((f3dotk*a1 + f2dotk)*a2 + fdot(k,i))*dt +
95  & f(k,i))*dt + x0dot(k,i))*dt
96  x0(k,i) = x(k,i) - dx0
97  x(k,i) = x0(k,i)
98  xdot(k,i) = x0dot(k,i)
99  65 CONTINUE
100  70 CONTINUE
101 *
102 * Check differential correction for Plummer potential.
103  IF (kz(14).EQ.4) THEN
104  CALL plpot1(phi2)
105  emdot = emdot + (phi1 - phi2)
106  END IF
107 *
108 * Ensure consistent coordinates & velocities for binary components.
109  DO 80 ipair = 1,npairs
110  IF (body(n+ipair).GT.0.0d0) THEN
111  CALL resolv(ipair,1)
112  END IF
113  80 CONTINUE
114 *
115  RETURN
116 *
117  END