Nbody6
fcorr.f
Go to the documentation of this file.
1  SUBROUTINE fcorr(I,DM,KW)
2 *
3 *
4 * Total force corrections due to masss loss.
5 * ------------------------------------------
6 *
7  include 'common6.h'
8  REAL*8 a(9)
9  LOGICAL ikick
10 *
11 *
12 * Save the velocity components and square velocity.
13  vi2 = 0.0
14  DO 1 k = 1,3
15  a(k+6) = xdot(k,i)
16  vi2 = vi2 + xdot(k,i)**2
17  1 CONTINUE
18 *
19 * Include velocity kick in case of new WD, NS, BH or massless SN.
20  IF (kw.GE.10.AND.kw.LE.15) THEN
21  ikick = .true.
22 * Ensure WD kicks are treated consistently (cf. define.f and #25).
23  IF ((kw.EQ.10.OR.kw.EQ.11).AND.kz(25).NE.1) ikick = .false.
24  IF (kw.EQ.12.AND.kz(25).NE.2) ikick = .false.
25 * Distinguish between single star (first time only) and binary.
26  IF (i.LE.n.AND.kw.NE.kstar(i).AND.ikick) THEN
27  CALL kick(i,1,kw,dm)
28  ELSE IF (i.GT.n.AND.ikick) THEN
29  ipair = i - n
30  CALL kick(ipair,0,kw,dm)
31  END IF
32  ELSE
33  ikick = .false.
34  END IF
35 *
36 * Define consistent c.m. variables for KS mass loss (exclude Roche).
37  IF (i.GT.n.AND.kstar(i).LE.10) THEN
38  i2 = 2*(i - n)
39  i1 = i2 - 1
40  vf2 = 0.0
41  dv2 = 0.0
42  bodyi = body(i)
43  IF (body(i).EQ.0.0d0) bodyi = body(i1) + body(i2)
44  DO 8 k = 1,3
45  x(k,i) = (body(i1)*x(k,i1) + body(i2)*x(k,i2))/bodyi
46  xdot(k,i) = (body(i1)*xdot(k,i1) + body(i2)*xdot(k,i2))/
47  & bodyi
48  x0(k,i) = x(k,i)
49  x0dot(k,i) = xdot(k,i)
50  vf2 = vf2 + xdot(k,i)**2
51  dv2 = dv2 + (xdot(k,i) - a(k+6))**2
52  8 CONTINUE
53  vfac = sqrt(vf2/vi2)
54  END IF
55 *
56 * Correct potential energy, all forces & first derivatives.
57  potj = 0.0d0
58  DO 40 j = ifirst,ntot
59  IF (j.EQ.i) go to 40
60  rij2 = 0.0d0
61  rijdot = 0.0d0
62  rdvdot = 0.0d0
63 *
64  DO 10 k = 1,3
65  a(k) = x(k,i) - x(k,j)
66  a(k+3) = a(k+6) - xdot(k,j)
67  rij2 = rij2 + a(k)**2
68  rijdot = rijdot + a(k)*a(k+3)
69  rdvdot = rdvdot + a(k)*(xdot(k,i) - a(k+6))
70  10 CONTINUE
71 *
72  rij = sqrt(rij2)
73  potj = potj + body(j)/rij
74  a3 = 1.0/(rij2*rij)
75  a4 = body(i)*a3
76  a5 = dm*a3
77  a6 = 3.0*rijdot/rij2
78  a7 = 3.0*rdvdot/rij2
79 *
80  DO 15 k = 1,3
81  a(k+3) = (a(k+3) - a6*a(k))*a5
82  IF (ikick) THEN
83 * Include FDOT corrections due to increased velocity.
84  a(k+3) = a(k+3) + (xdot(k,i) - a(k+6))*a4
85  a(k+3) = a(k+3) - a7*a(k)*a4
86  END IF
87  15 CONTINUE
88 *
89 * Use neighbour list of #J to distinguish irregular & regular terms.
90  nnb1 = list(1,j) + 1
91  DO 25 l = 2,nnb1
92  IF (list(l,j).EQ.i) THEN
93  DO 20 k = 1,3
94  f(k,j) = f(k,j) - 0.5*a(k)*a5
95  fi(k,j) = fi(k,j) - a(k)*a5
96  fdot(k,j) = fdot(k,j) - one6*a(k+3)
97  d1(k,j) = d1(k,j) - a(k+3)
98  fidot(k,j) = fidot(k,j) - a(k+3)
99  20 CONTINUE
100  go to 40
101  END IF
102  25 CONTINUE
103  DO 30 k = 1,3
104  f(k,j) = f(k,j) - 0.5*a(k)*a5
105  fr(k,j) = fr(k,j) - a(k)*a5
106  fdot(k,j) = fdot(k,j) - one6*a(k+3)
107  d1r(k,j) = d1r(k,j) - a(k+3)
108  frdot(k,j) = frdot(k,j) - a(k+3)
109  30 CONTINUE
110  40 CONTINUE
111 *
112 * Update the potential and kinetic energy loss.
113  emdot = emdot - dm*potj + 0.5*dm*vi2
114 *
115 * Modify energy loss further for c.m. body (exclude Roche cases).
116  IF (i.GT.n.AND.kstar(i).LE.10) THEN
117  ecdot = ecdot - 0.5*body(i)*vi2*(vfac**2 - 1.0)
118  END IF
119 *
120 * See whether linearized tidal terms should be included.
121  IF (kz(14).GT.0.AND.kz(14).LT.3) THEN
122  emdot = emdot - 0.5*dm*(tidal(1)*x(1,i)**2 +
123  & tidal(3)*x(3,i)**2)
124  END IF
125 *
126 * Check optional Plummer potential.
127  IF (kz(14).EQ.4.OR.kz(14).EQ.3) THEN
128  ri2 = ap2
129  DO 50 k = 1,3
130  ri2 = ri2 + x(k,i)**2
131  50 CONTINUE
132  emdot = emdot - dm*mp/sqrt(ri2)
133  END IF
134 *
135 * Accumulate energy loss for conservation check (not used).
136  e(12) = emdot
137 *
138  RETURN
139 *
140  END