Nbody6
 All Files Functions Variables
fpoly1.f
Go to the documentation of this file.
1  SUBROUTINE fpoly1(I1,I2,KCASE)
2 *
3 *
4 * Force & first derivative.
5 * -------------------------
6 *
7  include 'common6.h'
8  REAL*8 a(9),f1(3),f1dot(3)
9 *
10 *
11 * Standard case, new c.m. or KS termination (KCASE = 0, 1, 2).
12  jlast = ntot
13 * Reduce loop size for new c.m. polynomial.
14  IF (kcase.EQ.1) jlast = ntot - 1
15 *
16 * Loop over all bodies, pair #ICOMP & JCOMP or one single body.
17  DO 40 i = i1,i2
18 *
19 * Initialize forces & first differences for body #I.
20  DO 10 k = 1,3
21  fi(k,i) = 0.0d0
22  fr(k,i) = 0.0d0
23  d1(k,i) = 0.0d0
24  d1r(k,i) = 0.0d0
25  10 CONTINUE
26 *
27 * Include exceptional case of no neighbours (bug fix 2/3/12).
28  IF (list(1,i).EQ.0) THEN
29  rs0 = rs(i)
30  CALL nblist(i,rs0)
31  END IF
32 *
33 * Obtain force & first derivative by summing over all bodies.
34  kdum = 0
35  nnb = list(1,i)
36 * Set index of first neighbour to be identified in force loop.
37  l = 2
38  namej = list(l,i)
39 *
40 * Sum over all other bodies.
41  DO 30 jdum = ifirst,jlast
42  IF (jdum.EQ.i) go to 30
43  j = jdum
44  IF (j.GT.n.AND.j.EQ.namej) THEN
45  jpair = j - n
46 * Use c.m. approximation for unperturbed binary.
47  IF (list(1,2*jpair-1).GT.0) THEN
48  kdum = 2*jpair - 1
49  j = kdum
50  END IF
51  END IF
52 *
53  12 DO 15 k = 1,3
54  a(k) = x(k,j) - x(k,i)
55  a(k+3) = xdot(k,j) - xdot(k,i)
56  15 CONTINUE
57 *
58  a(7) = 1.0/(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
59  a(8) = body(j)*a(7)*sqrt(a(7))
60  a(9) = 3.0*(a(1)*a(4) + a(2)*a(5) + a(3)*a(6))*a(7)
61 *
62  DO 20 k = 1,3
63  f1(k) = a(k)*a(8)
64  f1dot(k) = (a(k+3) - a(k)*a(9))*a(8)
65  20 CONTINUE
66 *
67 * See whether summation index is equal to either component.
68  IF (j.EQ.icomp.OR.j.EQ.jcomp) THEN
69  IF (kcase.EQ.1) go to 30
70 * Note that dominant terms cancel analytically in c.m. polynomial.
71  END IF
72 *
73  IF (jdum.NE.namej) THEN
74  DO 25 k = 1,3
75  fr(k,i) = fr(k,i) + f1(k)
76  d1r(k,i) = d1r(k,i) + f1dot(k)
77  25 CONTINUE
78  ELSE
79  DO 28 k = 1,3
80  fi(k,i) = fi(k,i) + f1(k)
81  d1(k,i) = d1(k,i) + f1dot(k)
82  28 CONTINUE
83 *
84  IF (j.EQ.kdum) THEN
85  j = j + 1
86  go to 12
87  END IF
88 *
89 * Advance the neighbour list until last member has been considered.
90  IF (l.LE.nnb) THEN
91  l = l + 1
92  namej = list(l,i)
93  END IF
94  END IF
95  30 CONTINUE
96 *
97 * Check option for interstellar clouds (force & first derivative).
98  IF (kz(13).GT.0.AND.time.GT.0.0d0) THEN
99  CALL fcloud(i,f1,f1dot,2)
100  END IF
101  40 CONTINUE
102 *
103 * Check option for external force.
104  IF (kz(14).GT.0) THEN
105  CALL xtrnld(i1,i2,1)
106  END IF
107 *
108 * Set total force & first derivative.
109  DO 50 i = i1,i2
110 *
111  DO 45 k = 1,3
112  f(k,i) = fi(k,i) + fr(k,i)
113  fdot(k,i) = d1(k,i) + d1r(k,i)
114  45 CONTINUE
115  50 CONTINUE
116 *
117 * Check case of new c.m. (KCASE = 1 with I1 = ICOMP, I2 = JCOMP).
118  IF (kcase.EQ.1) THEN
119 * Form total c.m. force & first derivative and also AC components.
120  a1 = body(icomp)
121  a2 = body(jcomp)
122  DO 60 k = 1,3
123  f(k,ntot) = (a1*f(k,icomp) + a2*f(k,jcomp))/body(ntot)
124  fdot(k,ntot) = (a1*fdot(k,icomp) + a2*fdot(k,jcomp))/
125  & body(ntot)
126  fi(k,ntot) = (a1*fi(k,icomp) + a2*fi(k,jcomp))/body(ntot)
127  d1(k,ntot) = (a1*d1(k,icomp) + a2*d1(k,jcomp))/body(ntot)
128  fr(k,ntot) = (a1*fr(k,icomp) + a2*fr(k,jcomp))/body(ntot)
129  d1r(k,ntot) = (a1*d1r(k,icomp) + a2*d1r(k,jcomp))/
130  & body(ntot)
131  60 CONTINUE
132  j1 = ntot
133  j2 = ntot
134  ELSE
135  j1 = i1
136  j2 = i2
137  END IF
138 *
139  DO 80 i = j1,j2
140  DO 70 k = 1,3
141  d0(k,i) = fi(k,i)
142  d0r(k,i) = fr(k,i)
143  fidot(k,i) = d1(k,i)
144  frdot(k,i) = d1r(k,i)
145  70 CONTINUE
146  80 CONTINUE
147 *
148  RETURN
149 *
150  END