Nbody6
 All Files Functions Variables
cmfirr.f
Go to the documentation of this file.
1  SUBROUTINE cmfirr(I,IPAIR,XI,XIDOT,FIRR,FD)
2 *
3 *
4 * Irregular c.m. force & derivative.
5 * ----------------------------------
6 *
7  include 'common6.h'
8  REAL*8 xi(3),xidot(3),firr(3),fd(3),dx(3),dv(3),fp(6),fpd(6)
9 *
10 *
11 * Initialize the perturbing force & first derivative.
12  DO 1 k = 1,6
13  fp(k) = 0.0d0
14  fpd(k) = 0.0d0
15  1 CONTINUE
16 *
17 * Set individual KS components and indicators for unresolved pair.
18  i2 = 2*ipair
19  i1 = i2 - 1
20  kdum = 0
21  ifp = 0
22 *
23 * Specify perturber indices for decision-making.
24  np = list(1,i1)
25  lp = 2
26  jp = list(2,i1)
27 *
28 * Perform irregular force loop for perturbed c.m.
29  nnb0 = list(1,i)
30  DO 20 ll = 2,nnb0+1
31  k = list(ll,i)
32 * Advance lower perturber index (includes possible old neighbour).
33  2 IF (k.GT.jp.AND.lp.LE.np) THEN
34  lp = lp + 1
35  jp = list(lp,i1)
36 * Include rare case of two consecutive previous neighbours.
37  go to 2
38  END IF
39 *
40 * Distinguish between nearby perturber and more distant neighbour.
41  IF (k.NE.jp) THEN
42  IF (k.LE.n) go to 10
43 * Check c.m. approximation (point-mass assumption OK for #I).
44  a1 = x(1,k) - xi(1)
45  a2 = x(2,k) - xi(2)
46  a3 = x(3,k) - xi(3)
47  rij2 = a1*a1 + a2*a2 + a3*a3
48  IF (rij2.GT.cmsep2*r(k-n)**2) go to 10
49  kdum = 2*(k - n) - 1
50  IF (list(1,kdum).EQ.0) go to 10
51  k = kdum
52  go to 10
53 * Consider more carefully the case of identified perturber.
54  ELSE
55 * Determine next perturber index (if any).
56  IF (lp.LE.np) THEN
57  lp = lp + 1
58  jp = list(lp,i1)
59  END IF
60 * Specify first KS component if #JP is perturbed c.m.
61  IF (k.GT.n) THEN
62  kdum = 2*(k - n) - 1
63  IF (list(1,kdum).GT.0) THEN
64  k = kdum
65  END IF
66  END IF
67  END IF
68 *
69 * Evaluate perturbation on first component due to body #K.
70  ifp = 1
71  4 dr2 = 0.0
72  drdv = 0.0
73  DO 5 l = 1,3
74  dx(l) = x(l,k) - x(l,i1)
75  dv(l) = xdot(l,k) - xdot(l,i1)
76  dr2 = dr2 + dx(l)**2
77  drdv = drdv + dx(l)*dv(l)
78  5 CONTINUE
79 *
80  dr2i = 1.0/dr2
81  dr3i = body(k)*dr2i*sqrt(dr2i)
82  drdv = 3.0*drdv*dr2i
83 *
84  DO 6 l = 1,3
85  fp(l) = fp(l) + dx(l)*dr3i
86  fpd(l) = fpd(l) + (dv(l) - dx(l)*drdv)*dr3i
87  6 CONTINUE
88 *
89 * Evaluate perturbation on second component due to body #K.
90  dr2 = 0.0
91  drdv = 0.0
92  DO 7 l = 1,3
93  dx(l) = x(l,k) - x(l,i2)
94  dv(l) = xdot(l,k) - xdot(l,i2)
95  dr2 = dr2 + dx(l)**2
96  drdv = drdv + dx(l)*dv(l)
97  7 CONTINUE
98 *
99  dr2i = 1.0/dr2
100  dr3i = body(k)*dr2i*sqrt(dr2i)
101  drdv = 3.0*drdv*dr2i
102 *
103  DO 8 l = 1,3
104  fp(l+3) = fp(l+3) + dx(l)*dr3i
105  fpd(l+3) = fpd(l+3) + (dv(l) - dx(l)*drdv)*dr3i
106  8 CONTINUE
107 *
108  IF (k.GT.kdum) go to 20
109  k = k + 1
110  go to 4
111 *
112 * Sum over components of pair #J or its c.m. using c.m. approximation.
113  10 dr2 = 0.0
114  drdv = 0.0
115  DO 12 l = 1,3
116  dx(l) = x(l,k) - xi(l)
117  dv(l) = xdot(l,k) - xidot(l)
118  dr2 = dr2 + dx(l)**2
119  drdv = drdv + dx(l)*dv(l)
120  12 CONTINUE
121 *
122  dr2i = 1.0/dr2
123  dr3i = body(k)*dr2i*sqrt(dr2i)
124  drdv = 3.0*drdv*dr2i
125 *
126  DO 14 l = 1,3
127  firr(l) = firr(l) + dx(l)*dr3i
128  fd(l) = fd(l) + (dv(l) - dx(l)*drdv)*dr3i
129  14 CONTINUE
130  IF (k.EQ.kdum) THEN
131  k = k + 1
132  go to 10
133  END IF
134 
135  20 CONTINUE
136 *
137 * Add mass-weighted perturbations to force & first derivative.
138  IF (ifp.GT.0) THEN
139  bodyin = 1.0/body(i)
140  DO 30 k = 1,3
141  firr(k) = (body(i1)*fp(k) + body(i2)*fp(k+3))*bodyin +
142  & firr(k)
143  fd(k) = (body(i1)*fpd(k) + body(i2)*fpd(k+3))*bodyin +
144  & fd(k)
145  30 CONTINUE
146  END IF
147 *
148 * Check force correction due to regularized chain.
149  IF (nch.GT.0) THEN
150  CALL kcpert(i,i1,firr,fd)
151  END IF
152 *
153  RETURN
154 *
155  END