Nbody6
kspert.f
Go to the documentation of this file.
1  SUBROUTINE kspert(I1,NNB0,XI,VI,FP,FD)
2 *
3 *
4 * Perturbation on KS pair.
5 * ------------------------
6 *
7  include 'common6.h'
8  common/chainc/ xc(3,ncmax),uc(3,ncmax),bodyc(ncmax),ich,
9  & listc(lmax)
10  REAL*8 xi(6),vi(6),fp(6),fd(6),tf(3),td(3)
11 *
12 *
13 * Initialize the perturbing force & first derivative.
14  DO 10 k = 1,6
15  fp(k) = 0.0d0
16  fd(k) = 0.0d0
17  10 CONTINUE
18 *
19 * Set index of the last single perturber.
20  nnb2 = nnb0 + 1
21  15 IF (list(nnb2,i1).LE.n) go to 20
22  nnb2 = nnb2 - 1
23  IF (nnb2.GT.1) go to 15
24 * Include special case of only c.m. perturbers.
25  go to 30
26 *
27 * Obtain the perturbation from single particles.
28  20 DO 25 l = 2,nnb2
29  k = list(l,i1)
30  a1 = x(1,k) - xi(1)
31  a2 = x(2,k) - xi(2)
32  a3 = x(3,k) - xi(3)
33  rij2 = a1*a1 + a2*a2 + a3*a3
34  a6 = body(k)/(rij2*sqrt(rij2))
35  fp(1) = fp(1) + a1*a6
36  fp(2) = fp(2) + a2*a6
37  fp(3) = fp(3) + a3*a6
38  v1 = xdot(1,k) - vi(1)
39  v2 = xdot(2,k) - vi(2)
40  v3 = xdot(3,k) - vi(3)
41  a9 = 3.0d0*(a1*v1 + a2*v2 + a3*v3)/rij2
42  fd(1) = fd(1) + (v1 - a1*a9)*a6
43  fd(2) = fd(2) + (v2 - a2*a9)*a6
44  fd(3) = fd(3) + (v3 - a3*a9)*a6
45 * Perturbation on first component.
46 *
47  a1 = x(1,k) - xi(4)
48  a2 = x(2,k) - xi(5)
49  a3 = x(3,k) - xi(6)
50  rij2 = a1*a1 + a2*a2 + a3*a3
51  a6 = body(k)/(rij2*sqrt(rij2))
52  fp(4) = fp(4) + a1*a6
53  fp(5) = fp(5) + a2*a6
54  fp(6) = fp(6) + a3*a6
55  v1 = xdot(1,k) - vi(4)
56  v2 = xdot(2,k) - vi(5)
57  v3 = xdot(3,k) - vi(6)
58  a9 = 3.0d0*(a1*v1 + a2*v2 + a3*v3)/rij2
59  fd(4) = fd(4) + (v1 - a1*a9)*a6
60  fd(5) = fd(5) + (v2 - a2*a9)*a6
61  fd(6) = fd(6) + (v3 - a3*a9)*a6
62 * Perturbation on second component.
63  25 CONTINUE
64 *
65 * See whether to include any remaining c.m. perturbers.
66  IF (nnb2.GT.nnb0) go to 40
67 *
68  30 kdum = 0
69 * Dummy index to enable summation of c.m. or resolved components.
70  nnb3 = nnb2 + 1
71  DO 35 l = nnb3,nnb0+1
72  k = list(l,i1)
73  a1 = x(1,k) - xi(1)
74  a2 = x(2,k) - xi(2)
75  a3 = x(3,k) - xi(3)
76  rij2 = a1*a1 + a2*a2 + a3*a3
77 * See whether c.m. approximation applies (ignore unperturbed case).
78  j = k - n
79  IF (rij2.GT.cmsep2*r(j)**2) THEN
80  v1 = xdot(1,k) - vi(1)
81  v2 = xdot(2,k) - vi(2)
82  v3 = xdot(3,k) - vi(3)
83  a9 = 3.0d0*(a1*v1 + a2*v2 + a3*v3)/rij2
84  go to 33
85  END IF
86 *
87 * Resolve pair #J and sum over individual components.
88  CALL ksres2(j,j1,j2,rij2)
89  kdum = j1
90  k = kdum
91  32 a1 = x(1,k) - xi(1)
92  a2 = x(2,k) - xi(2)
93  a3 = x(3,k) - xi(3)
94  rij2 = a1*a1 + a2*a2 + a3*a3
95  v1 = xdot(1,k) - vi(1)
96  v2 = xdot(2,k) - vi(2)
97  v3 = xdot(3,k) - vi(3)
98  a9 = 3.0d0*(a1*v1 + a2*v2 + a3*v3)/rij2
99  33 a6 = body(k)/(rij2*sqrt(rij2))
100  fp(1) = fp(1) + a1*a6
101  fp(2) = fp(2) + a2*a6
102  fp(3) = fp(3) + a3*a6
103  fd(1) = fd(1) + (v1 - a1*a9)*a6
104  fd(2) = fd(2) + (v2 - a2*a9)*a6
105  fd(3) = fd(3) + (v3 - a3*a9)*a6
106 * Perturbation on first component.
107 *
108  a1 = x(1,k) - xi(4)
109  a2 = x(2,k) - xi(5)
110  a3 = x(3,k) - xi(6)
111  rij2 = a1*a1 + a2*a2 + a3*a3
112  a6 = body(k)/(rij2*sqrt(rij2))
113  fp(4) = fp(4) + a1*a6
114  fp(5) = fp(5) + a2*a6
115  fp(6) = fp(6) + a3*a6
116  v1 = xdot(1,k) - vi(4)
117  v2 = xdot(2,k) - vi(5)
118  v3 = xdot(3,k) - vi(6)
119  a9 = 3.0d0*(a1*v1 + a2*v2 + a3*v3)/rij2
120  fd(4) = fd(4) + (v1 - a1*a9)*a6
121  fd(5) = fd(5) + (v2 - a2*a9)*a6
122  fd(6) = fd(6) + (v3 - a3*a9)*a6
123 * Perturbation on second component.
124 *
125  IF (k.EQ.kdum) THEN
126  k = k + 1
127  go to 32
128  END IF
129  35 CONTINUE
130 *
131 * Check perturbation correction due to regularized chain.
132  40 IF (nch.GT.0) THEN
133  DO 45 l = 2,nnb2
134  j = list(l,i1)
135  IF (j.GT.ich) go to 50
136  IF (j.EQ.ich) THEN
137  j1 = i1
138  CALL fchain(j1,1,xi(1),vi(1),fp(1),fd(1))
139  j1 = j1 + 1
140  CALL fchain(j1,1,xi(4),vi(4),fp(4),fd(4))
141  go to 50
142  END IF
143  45 CONTINUE
144  END IF
145 *
146 * Set the relative perturbing force and first derivative.
147  50 DO 55 k = 1,3
148  fp(k) = fp(k) - fp(k+3)
149  fd(k) = fd(k) - fd(k+3)
150  tf(k) = 0.0d0
151  td(k) = 0.0d0
152  55 CONTINUE
153 *
154 * See whether the linearized perturbation should be included.
155  IF (kz(14).GT.0.AND.kz(14).LT.3) THEN
156  q1 = xi(1) - xi(4)
157  q3 = xi(3) - xi(6)
158  CALL xtrnlp(q1,q3,tf)
159 *
160 * Use same formalism for the first derivative (omit Coriolis force).
161  vx = vi(1) - vi(4)
162  vz = vi(3) - vi(6)
163  CALL xtrnlp(vx,vz,td)
164  DO 60 k = 1,3
165  fp(k) = fp(k) + tf(k)
166  fd(k) = fd(k) + td(k)
167  60 CONTINUE
168  END IF
169 *
170 * Check optional Plummer potential.
171  IF (kz(14).EQ.4.OR.kz(14).EQ.3) THEN
172  ri2 = ap2
173  rrdot = 0.0
174 * Form one central distance and scalar product of relative motion.
175  DO 65 k = 1,3
176  ri2 = ri2 + xi(k)**2
177  rrdot = rrdot + (xi(k) - xi(k+3))*(vi(k) - vi(k+3))
178  65 CONTINUE
179  zf = 1.0/ri2
180 * Write current mass inside RI as MP*R3*ZF^{3/2} (Heggie & Hut p.73).
181  fmp = mp*zf*sqrt(zf)
182  DO 70 k = 1,3
183  xrel = xi(k) - xi(k+3)
184  vrel = vi(k) - vi(k+3)
185  fp(k) = fp(k) - xrel*fmp
186  fd(k) = fd(k) - (vrel - 3.0*rrdot*zf*xrel)*fmp
187  70 CONTINUE
188  END IF
189 *
190  RETURN
191 *
192  END