Nbody6
 All Files Functions Variables
kcpert.f
Go to the documentation of this file.
1  SUBROUTINE kcpert(II,I1,FIRR,FD)
2 *
3 *
4 * Differential force correction on active KS from chain.
5 * ------------------------------------------------------
6 *
7  include 'common6.h'
8  common/chainc/ xc(3,ncmax),uc(3,ncmax),bodyc(ncmax),ich,
9  & listc(lmax)
10  REAL*8 firr(3),fd(3),dx(3),dv(3),fp(3),fpd(3),fcm(3),fcmd(3)
11 *
12 *
13 * See whether perturber list contains chain c.m. body #ICH.
14  j = 0
15  nnb1 = list(1,i1) + 1
16  DO 1 l = 2,nnb1
17  jj = list(l,i1)
18  IF (jj.EQ.ich) j = jj
19  1 CONTINUE
20 * Exit if chain c.m. not identified.
21  IF (j.EQ.0) go to 30
22 *
23 * Obtain chain c.m. force on each KS component for diff correction.
24  bodyin = 1.0/body(ii)
25  i = i1
26  5 a1 = x(1,j) - x(1,i)
27  a2 = x(2,j) - x(2,i)
28  a3 = x(3,j) - x(3,i)
29  rij2 = a1*a1 + a2*a2 + a3*a3
30 *
31  dv(1) = xdot(1,j) - xdot(1,i)
32  dv(2) = xdot(2,j) - xdot(2,i)
33  dv(3) = xdot(3,j) - xdot(3,i)
34  dr2i = 1.0/rij2
35  dr3i = body(j)*dr2i*sqrt(dr2i)
36  drdv = 3.0*(a1*dv(1) + a2*dv(2) + a3*dv(3))*dr2i
37 *
38 * Save c.m. contribution (to be subtracted).
39  fcm(1) = a1*dr3i
40  fcm(2) = a2*dr3i
41  fcm(3) = a3*dr3i
42  fcmd(1) = (dv(1) - a1*drdv)*dr3i
43  fcmd(2) = (dv(2) - a2*drdv)*dr3i
44  fcmd(3) = (dv(3) - a3*drdv)*dr3i
45 *
46  DO 10 k = 1,3
47  fp(k) = 0.0
48  fpd(k) = 0.0
49  10 CONTINUE
50 *
51 * Form contributions from each chain member for each component.
52  DO 20 jj = 1,nch
53  dr2 = 0.0
54  drdv = 0.0
55  DO 15 l = 1,3
56  dx(l) = xc(l,jj) - x(l,i)
57  dv(l) = uc(l,jj) - xdot(l,i)
58  dr2 = dr2 + dx(l)**2
59  drdv = drdv + dx(l)*dv(l)
60  15 CONTINUE
61  dr2i = 1.0/dr2
62  dr3i = bodyc(jj)*dr2i*sqrt(dr2i)
63  drdv = 3.0*drdv*dr2i
64 *
65 * Sum force & first derivative.
66  DO 18 l = 1,3
67  fp(l) = fp(l) + dx(l)*dr3i
68  fpd(l) = fpd(l) + (dv(l) - dx(l)*drdv)*dr3i
69  18 CONTINUE
70  20 CONTINUE
71 *
72 * Add differential corrections (chain contribution minus c.m. term).
73  DO 25 k = 1,3
74  firr(k) = firr(k) + body(i)*(fp(k) - fcm(k))*bodyin
75  fd(k) = fd(k) + body(i)*(fpd(k) - fcmd(k))*bodyin
76  25 CONTINUE
77 *
78 * Treat the second component in a similar way.
79  IF (i.EQ.i1) THEN
80  i = i1 + 1
81  go to 5
82  END IF
83 *
84  30 RETURN
85 *
86  END