Nbody6
kspoly.f
Go to the documentation of this file.
1  SUBROUTINE kspoly(IPAIR,IMOD)
2 *
3 *
4 * Initialization of KS polynomials.
5 * ---------------------------------
6 *
7  include 'common6.h'
8  common/slow0/ range,islow(10)
9  REAL*8 a1(3,4),a(8),fp(4),fpdot(3),ui(4),uidot(4)
10  REAL*8 b2(3,4),u2(4),q2(4)
11 *
12 *
13 * Specify indices of c.m. & components and perturber membership (+ 1).
14  i = n + ipair
15  i2 = 2*ipair
16  i1 = i2 - 1
17  nnb1 = list(1,i1) + 1
18  semi = -0.5*body(i)/h(ipair)
19 * Save new slow-down index (note direct CALL KSPOLY).
20  kslow(ipair) = imod
21 *
22 * Predict current coordinates & velocities of the perturbers.
23  DO 5 l = 2,nnb1
24  j = list(l,i1)
25  CALL xvpred(j,0)
26  5 CONTINUE
27 *
28 * Ensure current values of existing KS solution.
29  IF (iphase.EQ.0) THEN
30  CALL resolv(ipair,1)
31  END IF
32 *
33 * Initialize variables for accumulating contributions.
34  DO 10 k = 1,3
35  fp(k) = 0.0d0
36  fpdot(k) = 0.0d0
37  10 CONTINUE
38 * Set dummy index for summation of c.m. or resolved components.
39  jdum = 0
40 *
41 * Obtain the perturbation & first derivative.
42  DO 40 l = 2,nnb1
43  j = list(l,i1)
44  IF (j.GT.n) THEN
45 * See whether c.m. approximation applies.
46  rij2 = (x(1,j) - x(1,i))**2 + (x(2,j) - x(2,i))**2 +
47  & (x(3,j) - x(3,i))**2
48  k = j - n
49  IF (rij2.LT.cmsep2*r(k)**2.AND.list(1,2*k-1).GT.0) THEN
50  jdum = 2*k - 1
51  j = jdum
52  END IF
53  END IF
54 *
55 * Sum over both components (reverse sign for second component).
56  20 ii = i1
57  DO 30 kcomp = 1,2
58  DO 25 k = 1,3
59  a(k) = x(k,j) - x(k,ii)
60  a(k+3) = xdot(k,j) - xdot(k,ii)
61  25 CONTINUE
62 * Current velocities are predicted in routines INTGRT, KSMOD, etc.
63  rij2 = a(1)*a(1) + a(2)*a(2) + a(3)*a(3)
64  a8 = body(j)/(rij2*sqrt(rij2))
65  a9 = 3.0d0*(a(1)*a(4) + a(2)*a(5) + a(3)*a(6))/rij2
66  IF (kcomp.EQ.2) a8 = -a8
67 *
68  DO 28 k = 1,3
69  fp(k) = fp(k) + a(k)*a8
70  fpdot(k) = fpdot(k) + (a(k+3) - a(k)*a9)*a8
71  28 CONTINUE
72  ii = i2
73  30 CONTINUE
74 *
75  IF (j.EQ.jdum) THEN
76  j = j + 1
77  go to 20
78  END IF
79  40 CONTINUE
80 *
81 * Check for external perturbations.
82  IF (kz(14).GT.0.AND.kz(14).LT.3) THEN
83  q1 = x(1,i1) - x(1,i2)
84  q3 = x(3,i1) - x(3,i2)
85  CALL xtrnlp(q1,q3,fp)
86 *
87 * Use same formalism for the first derivative (omit Coriolis force).
88  vx = xdot(1,i1) - xdot(1,i2)
89  vz = xdot(3,i1) - xdot(3,i2)
90  CALL xtrnlp(vx,vz,fpdot)
91  END IF
92 *
93 * Transform to regularized force derivative using T' = R.
94  DO 45 k = 1,3
95  fpdot(k) = r(ipair)*fpdot(k)
96  45 CONTINUE
97 *
98 * Save the relative perturbation.
99  fp(4) = sqrt(fp(1)**2 + fp(2)**2 + fp(3)**2)
100  gamma(ipair) = fp(4)*r(ipair)**2/body(i)
101 *
102 * Copy U & UDOT to scalars and set current transformation matrix.
103  DO 48 k = 1,4
104  ui(k) = u(k,ipair)
105  uidot(k) = udot(k,ipair)
106  48 CONTINUE
107  CALL matrix(ui,a1)
108 *
109 * Construct regularized polynomials from explicit derivatives.
110  tdot2(ipair) = 0.0d0
111  tdot3(ipair) = 0.0d0
112  tdot4 = 0.0d0
113  tdot5 = 0.0d0
114  tdot6 = 0.0d0
115  hdot(ipair) = 0.0d0
116  hdot2(ipair) = 0.0d0
117  hdot3(ipair) = 0.0d0
118  hdot4(ipair) = 0.0d0
119 *
120 * Scale perturbing force & first derivative by modification factor.
121  IF (imod.GT.1) THEN
122  zmod = float(islow(imod))
123  DO 50 k = 1,3
124  fp(k) = zmod*fp(k)
125  fpdot(k) = zmod*fpdot(k)
126  50 CONTINUE
127  END IF
128 *
129 * Form regularized force & two derivatives of time & binding energy.
130  DO 60 k = 1,4
131  a(k) = a1(1,k)*fp(1) + a1(2,k)*fp(2) + a1(3,k)*fp(3)
132  a(k+4) = a1(1,k)*fpdot(1) + a1(2,k)*fpdot(2) +
133  & a1(3,k)*fpdot(3)
134  fu(k,ipair) = 0.5d0*h(ipair)*u(k,ipair) + 0.5d0*r(ipair)*a(k)
135  tdot2(ipair) = tdot2(ipair) + 2.0d0*u(k,ipair)*udot(k,ipair)
136  tdot3(ipair) = tdot3(ipair) + 2.0d0*udot(k,ipair)**2 +
137  & 2.0d0*u(k,ipair)*fu(k,ipair)
138  hdot(ipair) = hdot(ipair) + 2.0d0*udot(k,ipair)*a(k)
139  hdot2(ipair) = hdot2(ipair) + 2.0d0*fu(k,ipair)*a(k)
140  fp0(k,ipair) = 0.5d0*r(ipair)*a(k)
141  u2(k) = fu(k,ipair)
142  60 CONTINUE
143 *
144 * Set regularized velocity matrix (Levi-Civita matrix not required).
145  CALL matrix(uidot,a1)
146 *
147 * Include the whole (L*F)' term in explicit derivatives of FU & H.
148  DO 65 k = 1,4
149  a(k+4) = a(k+4) + a1(1,k)*fp(1) + a1(2,k)*fp(2) +
150  & a1(3,k)*fp(3)
151  hdot2(ipair) = hdot2(ipair) + 2.0d0*udot(k,ipair)*a(k+4)
152  fudot(k,ipair) = 0.5d0*(h(ipair)*udot(k,ipair) +
153  & hdot(ipair)*u(k,ipair) +
154  & tdot2(ipair)*a(k) + r(ipair)*a(k+4))
155  hdot3(ipair) = hdot3(ipair) + 2.0*fudot(k,ipair)*a(k) +
156  & 4.0*fu(k,ipair)*a(k+4)
157  tdot4 = tdot4 + 2.0*fudot(k,ipair)*u(k,ipair) +
158  & 6.0*fu(k,ipair)*udot(k,ipair)
159  fd0(k,ipair) = 0.5d0*(tdot2(ipair)*a(k) + r(ipair)*a(k+4) +
160  & hdot(ipair)*u(k,ipair))
161  65 CONTINUE
162 *
163 * Form L(U'') and add two of the three terms in (L*F)'' (11/1/98).
164  CALL matrix(u2,b2)
165  DO 66 k = 1,4
166  q2(k) = b2(1,k)*fp(1) + b2(2,k)*fp(2) + b2(3,k)*fp(3) +
167  & 2.0*(a1(1,k)*fpdot(1) + a1(2,k)*fpdot(2) +
168  & a1(3,k)*fpdot(3))
169  hdot3(ipair) = hdot3(ipair) + 2.0*udot(k,ipair)*q2(k)
170  66 CONTINUE
171 *
172 * Form higher derivatives by boot-strapping (Q2 added 11/1/98).
173  DO 70 k = 1,4
174  fudot2(k,ipair) = 0.5*(h(ipair)*fu(k,ipair) +
175  & hdot2(ipair)*u(k,ipair) + tdot3(ipair)*a(k)) +
176  & tdot2(ipair)*a(k+4) + hdot(ipair)*udot(k,ipair)
177  fudot3(k,ipair) = 0.5*(h(ipair)*fudot(k,ipair) + tdot4*a(k) +
178  & hdot3(ipair)*u(k,ipair)) +
179  & 1.5*(hdot(ipair)*fu(k,ipair) +
180  & hdot2(ipair)*udot(k,ipair) +
181  & tdot3(ipair)*a(k+4))
182  fudot2(k,ipair) = fudot2(k,ipair) + 0.5*r(ipair)*q2(k)
183  fudot3(k,ipair) = fudot3(k,ipair) + 1.5*tdot2(ipair)*q2(k)
184  hdot3(ipair) = hdot3(ipair) + 2.0*udot(k,ipair)*q2(k)
185  hdot4(ipair) = hdot4(ipair) + 2.0*fudot2(k,ipair)*a(k) +
186  & 6.0*(fudot(k,ipair)*a(k+4) + fu(k,ipair)*q2(k))
187  tdot5 = tdot5 + 2.0*fudot2(k,ipair)*u(k,ipair) +
188  & 8.0*fudot(k,ipair)*udot(k,ipair) + 6.0*fu(k,ipair)**2
189  tdot6 = tdot6 + 2.0*fudot3(k,ipair)*u(k,ipair) +
190  & 10.0*fudot2(k,ipair)*udot(k,ipair) +
191  & 20.0*fudot(k,ipair)*fu(k,ipair)
192  70 CONTINUE
193 *
194 * Check maximum square step (soft binaries & weak hyperbolic pairs).
195  IF (abs(h(ipair)).GT.eclose) THEN
196  a2 = 0.5/abs(h(ipair))
197  ELSE
198  a2 = min(r(ipair)/body(i),0.5/abs(h(ipair)))
199  END IF
200 *
201 * Assign a conservative value of the initial step (except IMOD > 1).
202  fac = 0.8
203  IF (imod.GT.1.OR.h(ipair).GT.0.0) fac = 1.0
204  dtu = fac*etau*sqrt(a2)/(1.0 + 1000.0*gamma(ipair))**0.333
205 *
206 * Include convergence criterion DH = H'*DTU + H''*DTU**2/2 = 0.001*|H|.
207  IF (imod.EQ.1.AND.hdot2(ipair).NE.0.0d0) THEN
208  dh = 1.0e-03*max(abs(h(ipair)),eclose)
209  xf = 2.0*dh/abs(hdot2(ipair))
210  yf = hdot(ipair)/hdot2(ipair)
211  dtu1 = sqrt(xf + yf**2) - abs(yf)
212  dtu = min(dtu1,dtu)
213  END IF
214 *
215 * Initialize reference energy.
216  h0(ipair) = h(ipair)
217 *
218 * Skip Stumpff part and use ZMOD*TK for unperturbed binary.
219  IF (list(1,i1).EQ.0.AND.semi.GT.0.0) THEN
220  tk = twopi*semi*sqrt(semi/body(i))
221  step(i1) = min(tk,step(i))
222  go to 80
223  END IF
224 *
225 * Obtain Stumpff coefficients.
226  75 z = -0.5d0*h(ipair)*dtu**2
227  CALL stumpf(ipair,z)
228  z5 = sf(6,ipair)
229  z6 = sf(7,ipair)
230  s6 = z6*tdot6/720.0
231 *
232 * Convert regularized step to physical units and initialize time T0.
233  step(i1) = ((((s6 + tdot5*dtu*z5/120.0d0 + tdot4/24.0d0)*dtu +
234  & tdot3(ipair)*one6)*dtu + 0.5d0*tdot2(ipair))*dtu +
235  & r(ipair))*dtu
236 *
237 * Ensure that c.m. step is not exceeded (H > 0 is OK).
238  IF (step(i1).GT.step(i).AND.h(ipair).LT.0.0) THEN
239  dtu = 0.5*dtu
240  go to 75
241  END IF
242 *
243 * Initialize step & time and include possible slow-down factor.
244  80 dtau(ipair) = dtu
245  t0(i1) = time
246  IF (imod.GT.1) THEN
247  step(i1) = zmod*step(i1)
248  END IF
249 *
250 * Include factorials in force and first derivative.
251  DO 85 k = 1,4
252  fu(k,ipair) = 0.5d0*fu(k,ipair)
253  fudot(k,ipair) = one6*fudot(k,ipair)
254  85 CONTINUE
255 *
256  RETURN
257 *
258  END