Nbody6
 All Files Functions Variables
cmfreg.f
Go to the documentation of this file.
1  SUBROUTINE cmfreg(I,RS2,RCRIT2,VRFAC,NNB,XI,XID,FIRR,FREG,FD,FDR)
2 *
3 *
4 * Regular & irregular c.m. force & first derivative.
5 * --------------------------------------------------
6 *
7  include 'common6.h'
8  REAL*8 xi(3),xid(3),firr(3),freg(3),dv(3),fd(3),fdr(3)
9 *
10 *
11 * Set non-zero indicator for perturbed c.m.
12  np = 0
13  IF (i.GT.n) THEN
14  ipair = i - n
15  IF (list(1,2*ipair-1).GT.0) np = 1
16  END IF
17 *
18 * Prepare case of single particle or unperturbed c.m. (second call).
19  IF (i.LE.n.OR.np.EQ.0) THEN
20 * Copy all KS pairs to JLIST and find MAX(R) for joint treatment.
21  nnb1 = npairs
22  rmax1 = 0.0
23  DO 10 lj = 1,nnb1
24  jlist(lj) = n + lj
25  rmax1 = max(rmax1,r(lj))
26  10 CONTINUE
27 *
28 * Adopt adequate square distance for c.m. approximation.
29  rcm2 = max(rcrit2,cmsep2*rmax1**2)
30 * Define dummy indices for skipping perturber test.
31  jp = 0
32  lp = 1
33  go to 25
34  END IF
35 *
36 * Specify variables for treatment of perturbed c.m. particle.
37  i2 = 2*ipair
38  i1 = i2 - 1
39  rpert2 = cmsep2*r(ipair)**2
40  bodyin = 1.0/body(i)
41 * Initialize perturber list for decision-making.
42  np = list(1,i1)
43  lp = 2
44  jp = list(2,i1)
45 *
46 * Use fast force loop for particles satisfying c.m. approximation.
47  rcm2 = max(rcrit2,rpert2)
48  nnb1 = 0
49  DO 20 j = ifirst,ntot
50  a1 = x(1,j) - xi(1)
51  a2 = x(2,j) - xi(2)
52  a3 = x(3,j) - xi(3)
53  dv(1) = xdot(1,j) - xid(1)
54  dv(2) = xdot(2,j) - xid(2)
55  dv(3) = xdot(3,j) - xid(3)
56  rij2 = a1*a1 + a2*a2 + a3*a3
57 *
58 * Form a list of particles for more careful consideration.
59  IF (rij2.LT.rcm2) THEN
60  nnb1 = nnb1 + 1
61  jlist(nnb1) = j
62  go to 20
63  END IF
64 *
65  dr2i = 1.0/rij2
66  dr3i = body(j)*dr2i*sqrt(dr2i)
67  drdv = 3.0*(a1*dv(1) + a2*dv(2) + a3*dv(3))*dr2i
68 *
69  freg(1) = freg(1) + a1*dr3i
70  freg(2) = freg(2) + a2*dr3i
71  freg(3) = freg(3) + a3*dr3i
72  fdr(1) = fdr(1) + (dv(1) - a1*drdv)*dr3i
73  fdr(2) = fdr(2) + (dv(2) - a2*drdv)*dr3i
74  fdr(3) = fdr(3) + (dv(3) - a3*drdv)*dr3i
75  20 CONTINUE
76 *
77 * Begin dual purpose force loop (all RIJ2 < RCM2, J > N or I <= N).
78  25 DO 60 lj = 1,nnb1
79  jdum = jlist(lj)
80  a1 = x(1,jdum) - xi(1)
81  a2 = x(2,jdum) - xi(2)
82  a3 = x(3,jdum) - xi(3)
83  dv(1) = xdot(1,jdum) - xid(1)
84  dv(2) = xdot(2,jdum) - xid(2)
85  dv(3) = xdot(3,jdum) - xid(3)
86  rij2 = a1*a1 + a2*a2 + a3*a3
87 * First see if the distance exceeds c.m. approximation limit.
88  IF (rij2.GT.rcm2) go to 56
89 *
90  j = jdum
91 * Check whether particle #J satisfies neighbour criteria.
92  IF (rij2.GT.rcrit2) go to 54
93 *
94 * Consider small step particle (may give large correction terms).
95  IF (rij2.GT.rs2.AND.step(j).GT.smin) THEN
96  a7 = a1*dv(1) + a2*dv(2) + a3*dv(3)
97 * Accept member if maximum penetration factor exceeds 8 per cent.
98  IF (a7.GT.vrfac) go to 54
99  END IF
100  IF (jdum.EQ.i) go to 60
101 *
102 * Obtain force due to current neighbours.
103  nnb = nnb + 1
104  ilist(nnb) = j
105  kcm = 1
106 *
107 * Advance lower perturber index (includes possible old neighbour).
108  26 IF (lp.LE.np.AND.j.GT.jp) THEN
109  lp = lp + 1
110  jp = list(lp,i1)
111 * Include rare case of two consecutive previous neighbours.
112  go to 26
113  END IF
114 *
115 * Decide appropriate expressions from perturber comparison.
116  IF (j.NE.jp) THEN
117  IF (j.LE.n) go to 30
118  IF (rij2.GT.cmsep2*r(j-n)**2) go to 30
119  kdum = 2*(j - n) - 1
120  IF (list(1,kdum).GT.0) THEN
121  k = kdum
122  j2 = k + 1
123  go to 50
124  END IF
125  ELSE
126 * Treat perturbers more carefully.
127  IF (lp.LE.np) THEN
128  lp = lp + 1
129  jp = list(lp,i1)
130  END IF
131  j2 = 0
132  IF (j.GT.n) THEN
133  kdum = 2*(j - n) - 1
134  IF (list(1,kdum).GT.0) THEN
135  j = kdum
136  j2 = j + 1
137  END IF
138  END IF
139  go to 40
140  END IF
141 *
142  30 dr2i = 1.0/rij2
143  dr3i = body(j)*dr2i*sqrt(dr2i)
144  drdv = 3.0*(a1*dv(1) + a2*dv(2) + a3*dv(3))*dr2i
145 *
146  firr(1) = firr(1) + a1*dr3i
147  firr(2) = firr(2) + a2*dr3i
148  firr(3) = firr(3) + a3*dr3i
149  fd(1) = fd(1) + (dv(1) - a1*drdv)*dr3i
150  fd(2) = fd(2) + (dv(2) - a2*drdv)*dr3i
151  fd(3) = fd(3) + (dv(3) - a3*drdv)*dr3i
152  go to 60
153 *
154 * Obtain relevant force on c.m (KCM = 0 denotes regular force).
155  40 k = j
156  42 l = i1
157 * Individual components I1 & I2 are resolved in routine INTGRT.
158  45 a1 = x(1,k) - x(1,l)
159  a2 = x(2,k) - x(2,l)
160  a3 = x(3,k) - x(3,l)
161  rij2 = a1*a1 + a2*a2 + a3*a3
162  dv(1) = xdot(1,k) - xdot(1,l)
163  dv(2) = xdot(2,k) - xdot(2,l)
164  dv(3) = xdot(3,k) - xdot(3,l)
165 *
166  dr2i = 1.0/rij2
167  dr3i = body(k)*body(l)*dr2i*sqrt(dr2i)*bodyin
168  drdv = 3.0*(a1*dv(1) + a2*dv(2) + a3*dv(3))*dr2i
169 *
170  IF (kcm.NE.0) THEN
171  firr(1) = firr(1) + a1*dr3i
172  firr(2) = firr(2) + a2*dr3i
173  firr(3) = firr(3) + a3*dr3i
174  fd(1) = fd(1) + (dv(1) - a1*drdv)*dr3i
175  fd(2) = fd(2) + (dv(2) - a2*drdv)*dr3i
176  fd(3) = fd(3) + (dv(3) - a3*drdv)*dr3i
177  ELSE
178  freg(1) = freg(1) + a1*dr3i
179  freg(2) = freg(2) + a2*dr3i
180  freg(3) = freg(3) + a3*dr3i
181  fdr(1) = fdr(1) + (dv(1) - a1*drdv)*dr3i
182  fdr(2) = fdr(2) + (dv(2) - a2*drdv)*dr3i
183  fdr(3) = fdr(3) + (dv(3) - a3*drdv)*dr3i
184  END IF
185 *
186  l = l + 1
187  IF (l.EQ.i2) go to 45
188  k = k + 1
189  IF (k.EQ.j2) go to 42
190  go to 60
191 *
192 * Treat c.m. approximation for #I and #K as single or composite.
193  50 a1 = x(1,k) - xi(1)
194  a2 = x(2,k) - xi(2)
195  a3 = x(3,k) - xi(3)
196  dv(1) = xdot(1,k) - xid(1)
197  dv(2) = xdot(2,k) - xid(2)
198  dv(3) = xdot(3,k) - xid(3)
199  rij2 = a1*a1 + a2*a2 + a3*a3
200 *
201  dr2i = 1.0/rij2
202  dr3i = body(k)*dr2i*sqrt(dr2i)
203  drdv = 3.0*(a1*dv(1) + a2*dv(2) + a3*dv(3))*dr2i
204 *
205  IF (kcm.NE.0) THEN
206  firr(1) = firr(1) + a1*dr3i
207  firr(2) = firr(2) + a2*dr3i
208  firr(3) = firr(3) + a3*dr3i
209  fd(1) = fd(1) + (dv(1) - a1*drdv)*dr3i
210  fd(2) = fd(2) + (dv(2) - a2*drdv)*dr3i
211  fd(3) = fd(3) + (dv(3) - a3*drdv)*dr3i
212  ELSE
213  freg(1) = freg(1) + a1*dr3i
214  freg(2) = freg(2) + a2*dr3i
215  freg(3) = freg(3) + a3*dr3i
216  fdr(1) = fdr(1) + (dv(1) - a1*drdv)*dr3i
217  fdr(2) = fdr(2) + (dv(2) - a2*drdv)*dr3i
218  fdr(3) = fdr(3) + (dv(3) - a3*drdv)*dr3i
219  END IF
220 *
221  k = k + 1
222  IF (k.EQ.j2) go to 50
223  go to 60
224 *
225 * Define regular force indicator.
226  54 kcm = 0
227 * Distinguish between second and first call (I > N & I <= N, J > N)
228  IF (jp.EQ.0) THEN
229 * Note that first case is for J > N and #I single or unperturbed c.m.
230  IF (rij2.LT.cmsep2*r(j-n)**2) THEN
231  j2 = 2*(j - n)
232  k = j2 - 1
233  go to 50
234  END IF
235  ELSE IF (j.LE.n) THEN
236 * Consider case of single #J and perturbed c.m.
237  IF (rij2.LT.rpert2) THEN
238  j2 = 0
239  go to 40
240  END IF
241  ELSE
242 * Split final case I > N & J > N into two parts according to RPERT2.
243  IF (rij2.GT.rpert2) THEN
244  IF (rij2.GT.cmsep2*r(j-n)**2) THEN
245  k = j
246  j2 = 0
247  ELSE
248  j2 = 2*(j - n)
249  k = j2 - 1
250  END IF
251 * Adopt c.m. approximation for #I.
252  go to 50
253  ELSE
254 * See whether both c.m. bodies should be resolved.
255  IF (rij2.GT.cmsep2*r(j-n)**2) THEN
256  j2 = 0
257  go to 40
258  ELSE
259  j2 = 2*(j - n)
260  k = j2 - 1
261  go to 42
262  END IF
263  END IF
264  END IF
265 *
266 * Obtain the regular force due to single body or c.m. particle.
267  56 dr2i = 1.0/rij2
268  dr3i = body(jdum)*dr2i*sqrt(dr2i)
269  drdv = 3.0*(a1*dv(1) + a2*dv(2) + a3*dv(3))*dr2i
270 *
271  freg(1) = freg(1) + a1*dr3i
272  freg(2) = freg(2) + a2*dr3i
273  freg(3) = freg(3) + a3*dr3i
274  fdr(1) = fdr(1) + (dv(1) - a1*drdv)*dr3i
275  fdr(2) = fdr(2) + (dv(2) - a2*drdv)*dr3i
276  fdr(3) = fdr(3) + (dv(3) - a3*drdv)*dr3i
277  60 CONTINUE
278 *
279 * Check force correction due to regularized chain (same as CMFIRR).
280  IF (i.GT.n.AND.nch.GT.0) THEN
281  IF (jp.GT.0) THEN
282  CALL kcpert(i,i1,firr,fd)
283  END IF
284  END IF
285 *
286  RETURN
287 *
288  END