Nbody6
 All Files Functions Variables
chfirr.f
Go to the documentation of this file.
1  SUBROUTINE chfirr(I,IR,XI,XIDOT,FIRR,FD)
2 *
3 *
4 * Irregular force & derivative on chain c.m.
5 * ------------------------------------------
6 *
7  include 'common6.h'
8  REAL*8 m,mass,mc,mij,mkk
9  parameter(nmx=10,nmx3=3*nmx,nmxm=nmx*(nmx-1)/2)
10  common/chain1/ xch(nmx3),vch(nmx3),m(nmx),
11  & zz(nmx3),wc(nmx3),mc(nmx),
12  & xj(nmx3),pi(nmx3),mass,rinv(nmxm),rsum,mkk(nmx),
13  & mij(nmx,nmx),tkk(nmx),tk1(nmx),iname(nmx),nn
14  common/chainc/ xc(3,ncmax),uc(3,ncmax),bodyc(ncmax),ich,
15  & listc(lmax)
16  REAL*8 xi(3),xidot(3),dx(3),dv(3),firr(3),fd(3),fp(3*ncmax),
17  & fpsum(3),fpd(3*ncmax),fpdsum(3)
18  SAVE nskip
19  DATA nskip /0/
20 *
21 *
22 * Update perturber list at frequent intervals.
23  IF (ir.GT.0) THEN
24  CALL chlist(ich)
25 * Note XC & UC for chain members predicted at new block-time.
26  ELSE
27  nskip = nskip + 1
28  IF (nskip.GE.10) THEN
29  CALL chlist(ich)
30  nskip = 0
31  END IF
32  END IF
33 *
34 * Subtract all close perturber forces acting on c.m. body #ICH.
35  rpert2 = cmsep2*(0.5*rsum)**2
36  npc = listc(1) + 1
37  jdum = 0
38  DO 30 l = 2,npc
39  j = listc(l)
40  dx(1) = x(1,j) - xi(1)
41  dx(2) = x(2,j) - xi(2)
42  dx(3) = x(3,j) - xi(3)
43  rij2 = dx(1)**2 + dx(2)**2 + dx(3)**2
44 * Skip correcting all interactions outside RPERT2.
45  IF (rij2.GT.rpert2) go to 30
46  IF (j.LE.n) go to 5
47  jp = j - n
48  j1 = 2*jp - 1
49 * See whether c.m. approximation applies.
50  IF (rij2.GT.cmsep2*r(jp)**2.OR.list(1,j1).EQ.0) go to 5
51 * Resolve components of pair #JP.
52  CALL ksres2(jp,j1,j2,rij2)
53 * Note danger of low-order prediction (high-order KSRES2 in INTGRT).
54  j = j1
55  jdum = j1
56 * Sum over single body or individual KS components.
57  5 dr2 = 0.0
58  drdv = 0.0
59  DO 10 k = 1,3
60  dx(k) = x(k,j) - xi(k)
61  dv(k) = xdot(k,j) - xidot(k)
62  dr2 = dr2 + dx(k)**2
63  drdv = drdv + dx(k)*dv(k)
64  10 CONTINUE
65  dr2i = 1.0/dr2
66  dr3i = body(j)*dr2i*sqrt(dr2i)
67  drdv = 3.0*drdv*dr2i
68 *
69 * Subtract force and first derivative.
70  DO 20 k = 1,3
71  firr(k) = firr(k) - dx(k)*dr3i
72  fd(k) = fd(k) - (dv(k) - dx(k)*drdv)*dr3i
73  20 CONTINUE
74  IF (j.EQ.jdum) THEN
75  j = j + 1
76  go to 5
77  END IF
78  30 CONTINUE
79 *
80 * Initialize the perturbing force components and vectorial sum.
81  DO 35 k = 1,3*nch
82  fp(k) = 0.0d0
83  fpd(k) = 0.0d0
84  35 CONTINUE
85  DO 40 k = 1,3
86  fpsum(k) = 0.0d0
87  fpdsum(k) = 0.0d0
88  40 CONTINUE
89 *
90 * Sum each perturber contribution over the chain components.
91  kdum = 0
92  im1 = 0
93  DO 70 ll = 2,npc
94  k = listc(ll)
95  a1 = x(1,k) - xi(1)
96  a2 = x(2,k) - xi(2)
97  a3 = x(3,k) - xi(3)
98  rij2 = a1*a1 + a2*a2 + a3*a3
99 *
100 * Skip perturbers consistent with c.m. approximation test above.
101  IF (rij2.GT.rpert2) go to 70
102 * Decide appropriate summation (components or c.m. approximation).
103  IF (k.GT.n) THEN
104 * Resolve the binary inside 20*R (OK for tidal effect on chain).
105  IF (rij2.LT.400.0*r(k-n)**2) THEN
106  kdum = 2*(k - n) - 1
107  k = kdum
108  go to 55
109  END IF
110  END IF
111 *
112 * Obtain perturbation on each component.
113  DO 50 im = 1,nch
114  im1 = 3*(im - 1)
115  46 dr2 = 0.0
116  drdv = 0.0
117  DO 48 l = 1,3
118  dx(l) = x(l,k) - xc(l,im)
119  dv(l) = xdot(l,k) - uc(l,im)
120  dr2 = dr2 + dx(l)**2
121  drdv = drdv + dx(l)*dv(l)
122  48 CONTINUE
123  dr2i = 1.0/dr2
124  dr3i = body(k)*dr2i*sqrt(dr2i)
125  drdv = 3.0*drdv*dr2i
126 *
127 * Accumulate perturbing force & derivative.
128  DO 49 l = 1,3
129  fp(im1+l) = fp(im1+l) + dx(l)*dr3i
130  fpd(im1+l) = fpd(im1+l) + (dv(l) - dx(l)*drdv)*dr3i
131  49 CONTINUE
132 *
133  IF (k.EQ.kdum) THEN
134  k = k + 1
135  go to 46
136  ELSE
137  IF (kdum.GT.0) k = kdum
138  END IF
139  50 CONTINUE
140 *
141 * Reset dummy index after use (otherwise bug with two KS pairs).
142  kdum = 0
143  go to 70
144 *
145 * Sum over individual components of pair #J using c.m. approximation.
146  55 dr2 = 0.0
147  drdv = 0.0
148  DO 60 l = 1,3
149  dx(l) = x(l,k) - xi(l)
150  dv(l) = xdot(l,k) - xidot(l)
151  dr2 = dr2 + dx(l)**2
152  drdv = drdv + dx(l)*dv(l)
153  60 CONTINUE
154  dr2i = 1.0/dr2
155  dr3i = body(k)*dr2i*sqrt(dr2i)
156  drdv = 3.0*drdv*dr2i
157 *
158 * Add force & first derivative.
159  DO 65 l = 1,3
160  firr(l) = firr(l) + dx(l)*dr3i
161  fd(l) = fd(l) + (dv(l) - dx(l)*drdv)*dr3i
162  65 CONTINUE
163  IF (k.EQ.kdum) THEN
164  k = k + 1
165  go to 55
166  END IF
167  70 CONTINUE
168 *
169 * Add any perturbations on components to the c.m. force & derivative.
170  IF (im1.GT.0) THEN
171 * Sum individual perturbations.
172  DO 80 l = 1,nch
173  l1 = 3*(l - 1)
174  DO 75 k = 1,3
175  fpsum(k) = fpsum(k) + bodyc(l)*fp(l1+k)
176  fpdsum(k) = fpdsum(k) + bodyc(l)*fpd(l1+k)
177  75 CONTINUE
178  80 CONTINUE
179 *
180 * Include mass-weighted contributions.
181  bodyin = 1.0/body(i)
182  DO 90 k = 1,3
183  firr(k) = fpsum(k)*bodyin + firr(k)
184  fd(k) = fpdsum(k)*bodyin + fd(k)
185  90 CONTINUE
186  END IF
187 *
188  RETURN
189 *
190  END