Nbody6
 All Files Functions Variables
derqp.f
Go to the documentation of this file.
1  SUBROUTINE derqp(Q,CMX,ENERGY,P,CMV,CHTIME,DQ,DX,DE,DP,DV,DT)
2 *
3 *
4 * Derivatives of chain variables.
5 * -------------------------------
6 *
7  include 'commonc.h'
8  LOGICAL kslow,kcoll
9  REAL*8 ksch
10  common/slow1/ tk2(0:nmx),ejump,ksch(nmx),kslow,kcoll
11  common/cpert/ rgrav,gpert,ipert,npert
12  common/calls/ tpr,tkpr,step,ider,icall,nfn,nreg,iter,imcirc
13  common/ccoll2/ qk(nmx4),pk(nmx4),rij(nmx,nmx),SIZE(nmx),vstar1,
14  & ecoll1,rcoll,qperi,istar(nmx),icoll,isync,ndiss1
15  common/bssave/ ep(4),dsc,facm,tfac,itfac,jc
16  common/ebsave/ ebs
17  common/ksave/ k10,k20
18  REAL*8 q(nmx4),cmx(3),p(nmx4),cmv(3)
19  REAL*8 dq(nmx4),dx(3),dp(nmx4),dv(3),fcm(3),upr(nmx4)
20  REAL*8 w(nmx4),ak(nmx4),dk(nmx),fnc(nmx3),fxtnl(nmx3)
21  REAL*8 tp(nmx4),tq(nmx4),uq(nmx4),faux(4),aux(-2:2),xaux(3)
22  SAVE
23 *
24 *
25 * Mikkola & Aarseth 1990 eqs. (77) -> (80).
26  uc=0.0
27  rsum=0.0
28 * Eq. (77).
29  DO i=1,n-1
30  ks = 4*(i-1)
31  k1 = ks + 1
32  k2 = ks + 2
33  k3 = ks + 3
34  k4 = ks + 4
35 * Obtain W = L^T(Q)P ; L^T = transpose of L times P.
36  w(k1)=(q(k1)*p(k1)-q(k2)*p(k2)-q(k3)*p(k3)+q(k4)*p(k4))
37  w(k2)=(q(k2)*p(k1)+q(k1)*p(k2)-q(k4)*p(k3)-q(k3)*p(k4))
38  w(k3)=(q(k3)*p(k1)+q(k4)*p(k2)+q(k1)*p(k3)+q(k2)*p(k4))
39  w(k4)=(q(k4)*p(k1)-q(k3)*p(k2)+q(k2)*p(k3)-q(k1)*p(k4))
40  rijl=q(k1)**2+q(k2)**2+q(k3)**2+q(k4)**2
41 * Form RSUM for decision-making.
42  rsum=rsum+rijl
43  rinv(i)=1./rijl
44  a=.5d0*rinv(i)
45  uc=uc+mkk(i)*rinv(i)
46  DO k=1,4
47  w(ks+k)=a*w(ks+k)
48  END DO
49  END DO
50  lri=n-1
51 * Eq. (78).
52  tkin=0.0
53  DO i=1,n-1
54  aux(-2)=0.5d0*tk2(i-1)
55  aux(-1)=0.5d0*tk1(i)
56  aux( 0)=tkk(i)
57  aux(+1)=0.5d0*tk1(i+1)
58  aux(+2)=0.5d0*tk2(i+1)
59  l=4*(i-1)
60  dk(i)=0.0
61  DO k=1,4
62  aa=0.0
63  DO j=-2,2
64  lj=l+4*j
65  if (lj.ge.0.and.lj.le.4*n-8) then
66  aa=aa+aux(j)*w(lj+k)
67  end if
68  END DO
69  ak(l+k)=aa
70 * Eq. (79).
71  dk(i)=dk(i)+aa*w(l+k)
72  END DO
73 * Eq. (80).
74  tkin=tkin+dk(i)
75  END DO
76 *
77 * Obtain physical coordinates.
78  DO k=1,3
79  xi(k)=0.0
80  END DO
81  DO i=1,n-1
82  l=3*(i-1)
83  ks=4*(i-1)
84  xc(l+1)=q(ks+1)**2-q(ks+2)**2-q(ks+3)**2+q(ks+4)**2
85  xc(l+2)=2.d0*(q(ks+1)*q(ks+2)-q(ks+3)*q(ks+4))
86  xc(l+3)=2.d0*(q(ks+1)*q(ks+3)+q(ks+2)*q(ks+4))
87  DO k=1,3
88  xi(l+3+k)=xi(l+k)+xc(l+k)
89  END DO
90  END DO
91 *
92 * Evaluate external force (FXTNL = W', not the 'usual' force!).
93  IF (ipert.GT.0) THEN
94  iskip=0
95  CALL xtf(fxtnl,fcm,cmx,chtime)
96 * Note that perturbation limit must be consistent with chlist.f.
97  IF (gpert.LT.1.0d-06) THEN
98  ipert=0
99  ELSE
100  ipert=1
101  END IF
102  ELSE
103  iskip=1
104  DO i=1,3*(n-1)
105  fxtnl(i)=0.0d0
106  END DO
107  DO i=1,3
108  fcm(i)=0.0d0
109  END DO
110  END IF
111 *
112 * Form non-chained contributions.
113  unc=0.0
114  DO i=1,3*(n-1)
115  fnc(i)=fxtnl(i)
116  END DO
117 *
118  DO i=1,n-2
119  li=3*(i-1)
120  DO j=i+2,n
121  lj=3*(j-1)
122  rij2=0.0
123  IF(j.GT.i+2)THEN
124  DO k=1,3
125  xaux(k)=xi(lj+k)-xi(li+k)
126  rij2=rij2+xaux(k)**2
127  END DO
128  ELSE
129  DO k=1,3
130  xaux(k)=xc(li+k)+xc(li+k+3)
131  rij2=rij2+xaux(k)**2
132  END DO
133  END IF
134  rij2in=1./rij2
135 * Introduce the inverse distances.
136  lri=lri+1
137  rinv(lri)=sqrt(rij2in)
138 *
139  fm=mij(i,j)*rinv(lri)
140  unc=unc+fm
141  fm=fm*rij2in
142 * Fij attraction.
143  DO k=1,3
144  faux(k)=-fm*xaux(k)
145  END DO
146 * Add the contribution to interactions depending on Rij.
147  DO ik=i,j-1
148  l=3*(ik-1)
149  DO k=1,3
150  fnc(l+k)=fnc(l+k)+faux(k)
151  END DO
152  END DO
153  END DO
154  END DO
155 *
156 * Evaluate UQ & TP.
157  DO i=1,n-1
158  l1=3*(i-1)+1
159  ks=4*(i-1)
160  ks1=ks+1
161  CALL qforce(q(ks1),fnc(l1),uq(ks1))
162  CALL vector(q(ks1),ak(ks1),tp(ks1))
163 * The * operation of eq. (84).
164  ak(ks+4)=-ak(ks+4)
165  CALL vector(p(ks1),ak(ks1),tq(ks1))
166 *
167  DO k=1,4
168  uq(ks+k)=uq(ks+k)-2.0d0*mkk(i)*q(ks+k)*rinv(i)**2
169  tq(ks+k)=tq(ks+k)-4.d0*dk(i)*q(ks+k)
170  END DO
171  END DO
172 * NOTE: The division by R above (in TP & TQ) is delayed.
173 *
174 * Proceed to final evaluation of derivatives (90)->(94).
175  upot=uc+unc
176  g=1./(tkin+upot)
177  h=tkin-upot
178 * Reset EJUMP after significant slow-down change (Seppo's procedure).
179 * IF (KJUMP) THEN
180 * EJUMP = H - ENERGY
181 * KJUMP = .false.
182 * END IF
183  gamma=(h-(energy+ejump))*g
184 *
185  gt= (1.-gamma)*g
186  gu=-(1.+gamma)*g
187 *
188  DO i=1,n-1
189  ks=4*(i-1)
190 * Apply the division by R here (to TP & TQ).
191  gtoverr=gt*rinv(i)
192 * NOTE: TP & TQ never get 'correct' values (thus TP = R*TPtrue).
193  DO k=1,4
194  dq(ks+k)=gtoverr*tp(ks+k)
195  dp(ks+k)=-gtoverr*tq(ks+k)-gu*uq(ks+k)
196  END DO
197  END DO
198  dt=g
199  tpr=g
200  DO k=1,3
201  dx(k)=cmv(k)*g
202  dv(k)=fcm(k)*g
203  END DO
204 * Evaluate E'.
205  de=0.0
206  IF (iskip.EQ.0) THEN
207  DO i=1,n-1
208  l1=3*(i-1)+1
209  ks=4*(i-1)
210  ks1=ks+1
211  CALL qforce(q(ks1),fxtnl(l1),faux)
212  DO k=1,4
213  de=de+dq(ks+k)*faux(k)
214  END DO
215  END DO
216  END IF
217 *
218 * Copy the time derivative for step control and reset indicator.
219 * IF (IDER.GT.0) THEN
220 * TPR = G
221 * IDER = 0
222 * END IF
223 *
224 * Check osculating pericentre of closest pair (first call only).
225  IF (icall.EQ.0) go to 50
226  IF (jc.GT.0) go to 10
227 *
228 * Perform a fast pericentre search (saves unnecessary complications).
229  rm = 0.0
230  DO 5 i = 1,n-1
231  IF (rinv(i).GT.rm) THEN
232  rm = rinv(i)
233  im = i
234  END IF
235  5 CONTINUE
236 *
237  k1 = iname(im)
238  k2 = iname(im+1)
239  k = 4*(im - 1) + 1
240  CALL peri(q(k),dq(k),tpr/ksch(im),m(k1),m(k2),qperi)
241 *
242 * Switch off indicator and exit for large pericentre.
243  IF (qperi.GT.4.0*max(SIZE(k1),SIZE(k2))) THEN
244  icall = 0
245  go to 50
246  END IF
247 *
248 * Examine all close pair indices K1 & K2 in chain vector INAME.
249  10 qpmin = 100.0
250  rpmin = 100.0
251  im0 = 1
252  DO 20 im = 1,n-1
253  k1 = iname(im)
254  k2 = iname(im+1)
255 *
256 * Set distance of nearest perturber (note extra test for N > 4).
257  IF (im.EQ.1) THEN
258  rp = 1.0/rinv(2)
259  ELSE IF (im.EQ.2.AND.n.GT.3) THEN
260  rp = min(1.0/rinv(1),1.0/rinv(3))
261  ELSE IF (im.EQ.3.AND.n.GT.4) THEN
262  rp = min(1.0/rinv(2),1.0/rinv(4))
263  ELSE
264  rp = 1.0/rinv(im-1)
265  END IF
266  rpmin = min(rp,rpmin)
267 *
268 * Determine pericentre for small perturbations by Mikkola's algorithm.
269  gi = (1.0/(rinv(im)*rp))**3
270  IF (gi.LT.0.005) THEN
271  k = 4*(im - 1) + 1
272  CALL peri(q(k),dq(k),tpr/ksch(im),m(k1),m(k2),qperi)
273  ELSE
274  qperi = 1.0/rinv(im)
275  END IF
276 *
277 * Compare pericentre with previous mutual distances (note symmetry).
278  rij(k1,k2) = min(rij(k1,k2),qperi)
279  rij(k2,k1) = min(rij(k2,k1),qperi)
280 *
281 * Save indices for smallest pericentre and switch off indicator.
282  IF (qperi.LT.qpmin.AND.im.NE.imcirc) THEN
283  g0 = gi
284  qpmin = qperi
285  rp0 = rp
286  k10 = k1
287  k20 = k2
288  im0 = im
289  icall = 0
290  END IF
291  20 CONTINUE
292 *
293 * Check smallest pericentre of current chain and copy latest value.
294  rcoll = min(rcoll,qpmin)
295  qperi = qpmin
296 *
297 * Specify KS index and closest pair indices.
298  im = im0
299  ks1 = 4*(im - 1) + 1
300  k1 = k10
301  k2 = k20
302 *
303 * Save TPR and current configuration (for EREL, TRANSK & CHAIN).
304  tkpr = tpr
305  DO 25 i = 1,n-1
306  ks = 4*(i - 1)
307  DO 24 j = 1,4
308  qk(ks+j) = q(ks+j)
309  pk(ks+j) = p(ks+j)
310  24 CONTINUE
311  25 CONTINUE
312 *
313 * Check for tidal two-body interaction or stellar collision.
314  IF (qpmin.LT.4.0*max(SIZE(k1),SIZE(k2)).OR.iter.GT.0) THEN
315 *
316 * Terminate after 25 iterations (convergence problem).
317  iter = iter + 1
318  IF (iter.GE.25) THEN
319  WRITE (6,26) iter, im, g0, qpmin, ecc
320  26 FORMAT (' WARNING! NO CONVERGENCE # IM GI QP ECC ',
321  & i5,i4,1p,3e10.2)
322  jc = 0
323  kcoll = .false.
324  icoll = -1
325  iter = 0
326  imcirc = 0
327  go to 50
328  END IF
329 *
330 * Exit if collision candidate distance is not the smallest.
331  rb = 1.0/rinv(im)
332 * IF (RB.GT.1.001*RPMIN) THEN
333 * JC = 0
334 * DSC = 1.0
335 * GO TO 50
336 * END IF
337 *
338 * Convert Q' to standard KS with T' = R and form radial velocity R'.
339  rpr = 0.0d0
340  DO 30 j = ks1,ks1+3
341  upr(j) = dq(j)*rb*ksch(im)/tpr
342  rpr = rpr + 2.0d0*q(j)*upr(j)
343  30 CONTINUE
344 *
345 * Determine small semi-major axis from non-singular expressions.
346  CALL erel(im,eb,semi)
347 *
348 * Exclude circularized binaries and consider second smallest peri.
349  ecc = 1.0 - qpmin/semi
350 * Skip exclusion in GR case (denoted by large VSTAR1).
351  IF (ecc.LT.0.002.AND.imcirc.EQ.0.AND.vstar1.LT.100.0) THEN
352  iter = 0
353  imcirc = im
354  go to 10
355  END IF
356 *
357 * Temporary exit because SEMI < 0 does not converge (tested OK 3/99).
358 * IF (SEMI.LT.0.0) THEN
359 * JC = 0
360 * GO TO 50
361 * END IF
362 *
363 * Obtain pericentre time interval from elements & Kepler's equation.
364  mb = m(k1) + m(k2)
365  CALL tperi(semi,q(ks1),upr(ks1),mb,tper)
366 *
367 * Activate collision indicator & B-S step selector (first time).
368  IF (icoll.EQ.0) THEN
369  icoll = -1
370  jc = 1
371  kcoll = .true.
372  ELSE
373 *
374 * Check convergence: radial velocity < 1.0E-09 parabolic velocity.
375  IF (abs(rpr).LT.1.0e-09*sqrt(2.0d0*mb*rb)) THEN
376 * Reset B-S step selector and copy chain index to collision indicator.
377  jc = 0
378  icoll = im
379  ebs = eb
380  iter = 0
381  imcirc = 0
382  END IF
383  END IF
384 *
385 * Set regularized step for DIFSY1 using accelerated iteration rate.
386 * IF (RB.GT.2.0*QPMIN) THEN
387 * DSC = 2.0*ABS(TPER)/TPR
388 * ELSE
389 * DSC = ABS(TPER)/TPR
390 * END IF
391 *
392 * Evaluate regularized pericentre time (Stiefel & Scheifele, p. 85).
393  hi = -0.5*mb/semi
394 * Note use of Seppo's sign convention (TPER > 0 after peri).
395  dsc = 2.0d0*(hi*tper - 0.5d0*rpr)/mb
396 * Scale from KS to chain time derivative.
397  dsc = dsc*rb/tpr
398 *
399 * Ensure negative step beyond pericentre (case RPR > 0 works OK).
400  IF (jc.GT.0.AND.rpr.LT.0.0d0) THEN
401  dsc = abs(dsc)
402  END IF
403 * Restore step to dummy value at the end (not used).
404  IF (jc.EQ.0) dsc = 1.0
405 *
406 * Switch off iteration on large perturbation after five tries.
407  IF (g0.GT.0.005.AND.iter.GT.5) THEN
408  jc = 0
409  kcoll = .false.
410  icoll = -1
411  iter = 0
412  imcirc = 0
413 * Avoid apocentre region of secondary binary (algorithmic confusion).
414  ELSE IF (rb.GT.semi.AND.imcirc.GT.0) THEN
415  jc = 0
416  iter = 0
417  imcirc = 0
418  icoll = 0
419  kcoll = .false.
420  END IF
421  ELSE
422  jc = 0
423 * Enforce restart from saved variables on failed iteration.
424  kcoll = .false.
425  icoll = -1
426  icall = 0
427  iter = 0
428  END IF
429 *
430 * Increase function call counter.
431  50 nfn = nfn + 1
432 *
433  RETURN
434 *
435  END