Nbody6
derqp4.f
Go to the documentation of this file.
1  SUBROUTINE derqp4(P,Q,DP,DQ)
2 *
3 *
4 * Equations of motion for chain regularization.
5 * ---------------------------------------------
6 *
7  IMPLICIT REAL*8 (a-h,o-z)
8  REAL*8 m,mij,xnr(9),fnr(9),wi(9),at(3),aq(3),ap(9),fq(9),wp(9),
9  & p(12),q(12),dp(12),dq(13),upr(12),mb
10  LOGICAL switch,gtype,gtype0
11  common/creg/ m(4),x(12),xd(12),pp(12),qq(12),time4,energy,epsr2,
12  & xr(9),w(9),r(6),ta(6),mij(6),cm(10),rmax4,tmax,
13  & ds,tstep,eps,nstep4,name4(4),kz15,kz27,nreg,nfn
14  common/tpr/ switch,gtype,gtype0
15  common/iconf/ i1,i2,i3,i4
16  common/close/ rij(4,4),rcoll,qperi,SIZE(4),ecoll3,ip(4)
17  common/ccoll/ qk(12),pk(12),icall,icoll,ndiss4
18  common/bssave/ ep(4),dsc,facm,tfac,itfac,jc
19  common/ksave/ k10,k20
20 * EQUIVALENCE (T11,TA(1)),(T22,TA(2)),(T33,TA(3)),(T12,TA(4)),
21 * & (T23,TA(5))
22 *
23 * Copy old equivalence variables.
24  t11 = ta(1)
25  t22 = ta(2)
26  t33 = ta(3)
27  t12 = ta(4)
28  t23 = ta(5)
29 *
30  nfn = nfn + 1
31  k = 0
32  kp = 0
33  DO 1 l = 1,3
34  k1 = k + 1
35  k2 = k + 2
36  k3 = k + 3
37  k4 = k + 4
38  kp1 = kp + 1
39  kp2 = kp + 2
40  kp3 = kp + 3
41 *
42 * Form physical vector & scalar distance.
43  xr(kp1) = q(k1)**2 - q(k2)**2 - q(k3)**2 + q(k4)**2
44  xr(kp2) = 2.d0*(q(k1)*q(k2) - q(k3)*q(k4))
45  xr(kp3) = 2.d0*(q(k1)*q(k3) + q(k2)*q(k4))
46  r(l) = q(k1)**2 + q(k2)**2 + q(k3)**2 + q(k4)**2
47 * Set physical momentum times distance.
48  wi(kp1) = q(k1)*p(k1) - q(k2)*p(k2) - q(k3)*p(k3) +
49  & q(k4)*p(k4)
50  wi(kp2) = q(k2)*p(k1) + q(k1)*p(k2) - q(k4)*p(k3) -
51  & q(k3)*p(k4)
52  wi(kp3) = q(k3)*p(k1) + q(k4)*p(k2) + q(k1)*p(k3) +
53  & q(k2)*p(k4)
54  at(l) = ta(l)*(p(k1)**2 + p(k2)**2 + p(k3)**2 +
55  & p(k4)**2) - mij(l)
56  k = k4
57  kp = kp3
58  1 CONTINUE
59 *
60 * Obtain irregular vectors & distances.
61  r(4) = 0.0d0
62  r(5) = 0.0d0
63  r(6) = 0.0d0
64  DO 2 k = 1,3
65  xnr(k ) = xr(k ) + xr(k+3)
66  xnr(k+3) = xr(k+3) + xr(k+6)
67  xnr(k+6) = xnr(k ) + xr(k+6)
68  r(4) = r(4) + xnr(k )**2
69  r(5) = r(5) + xnr(k+3)**2
70  r(6) = r(6) + xnr(k+6)**2
71  2 CONTINUE
72 *
73 * Evaluate irregular physical forces.
74  unr = 0.0d0
75  DO 5 i = 4,6
76  dist = sqrt(r(i))
77  fc = mij(i)/dist
78  unr = unr + fc
79  fc = fc/r(i)
80  r(i) = dist
81  ik = 3*(i - 4)
82  DO 3 k = 1,3
83  fnr(ik+k) = fc*xnr(ik+k)
84  3 CONTINUE
85  5 CONTINUE
86  w1w2 = wi(1)*wi(4) + wi(2)*wi(5) + wi(3)*wi(6)
87  w2w3 = wi(4)*wi(7) + wi(5)*wi(8) + wi(6)*wi(9)
88  unre = unr + energy
89  r12 = r(1)*r(2)
90  r13 = r(1)*r(3)
91  r23 = r(2)*r(3)
92  aq(1) = at(2)*r(3) + at(3)*r(2) + t23*w2w3-r23*unre
93  aq(2) = at(1)*r(3) + at(3)*r(1) - r13*unre
94  aq(3) = at(2)*r(1) + at(1)*r(2) + t12*w1w2-r12*unre
95  r123 = r12*r(3)
96  ap(1) = 2.0d0*t11*r23
97  ap(2) = 2.0d0*t22*r13
98  ap(3) = 2.0d0*t33*r12
99  wk12 = t12*r(3)
100  wk23 = t23*r(1)
101  DO 6 k = 1,3
102  wp(k ) = wk12*wi(k+3)
103  wp(k+3) = wk12*wi(k ) + wk23*wi(k+6)
104  wp(k+6) = wk23*wi(k+3)
105  fq(k ) = fnr(k ) + fnr(k+6)
106  fq(k+3) = fnr(k ) + fnr(k+3) + fnr(k+6)
107  fq(k+6) = fnr(k+3) + fnr(k+6)
108  6 CONTINUE
109 *
110 * Form regularized derivatives.
111  kq = 0
112  l = 0
113  DO 10 i = 1,3
114  k1 = kq + 1
115  k2 = kq + 2
116  k3 = kq + 3
117  k4 = kq + 4
118  l1 = l + 1
119  l2 = l + 2
120  l3 = l + 3
121 *
122  f1 = +fq(l1)*q(k1) + fq(l2)*q(k2) + fq(l3)*q(k3)
123  f2 = -fq(l1)*q(k2) + fq(l2)*q(k1) + fq(l3)*q(k4)
124  f3 = -fq(l1)*q(k3) - fq(l2)*q(k4) + fq(l3)*q(k1)
125  f4 = +fq(l1)*q(k4) - fq(l2)*q(k3) + fq(l3)*q(k2)
126 *
127  g1 = +wp(l1)*p(k1) + wp(l2)*p(k2) + wp(l3)*p(k3)
128  g2 = -wp(l1)*p(k2) + wp(l2)*p(k1) + wp(l3)*p(k4)
129  g3 = -wp(l1)*p(k3) - wp(l2)*p(k4) + wp(l3)*p(k1)
130  g4 = +wp(l1)*p(k4) - wp(l2)*p(k3) + wp(l3)*p(k2)
131 *
132  rrr = r123 + r123
133  dp(k1) = -(2.d0*aq(i)*q(k1) + g1+rrr*f1)
134  dp(k2) = -(2.d0*aq(i)*q(k2) + g2+rrr*f2)
135  dp(k3) = -(2.d0*aq(i)*q(k3) + g3+rrr*f3)
136  dp(k4) = -(2.d0*aq(i)*q(k4) + g4+rrr*f4)
137 *
138  g1 = +wp(l1)*q(k1) + wp(l2)*q(k2) + wp(l3)*q(k3)
139  g2 = -wp(l1)*q(k2) + wp(l2)*q(k1) + wp(l3)*q(k4)
140  g3 = -wp(l1)*q(k3) - wp(l2)*q(k4) + wp(l3)*q(k1)
141  g4 = +wp(l1)*q(k4) - wp(l2)*q(k3) + wp(l3)*q(k2)
142 *
143  dq(k1) = ap(i)*p(k1) + g1
144  dq(k2) = ap(i)*p(k2) + g2
145  dq(k3) = ap(i)*p(k3) + g3
146  dq(k4) = ap(i)*p(k4) + g4
147 *
148  kq = kq + 4
149  l = l + 3
150  10 CONTINUE
151 *
152 * Set the time derivative and check tolerance scaling (ITFAC > 1).
153  dq(13) = r123
154  IF (itfac.GT.0) THEN
155  tfac = facm*(r(1)*r(2) + r(1)*r(3) + r(2)*r(3))
156  itfac = -1
157  END IF
158 *
159 * Procedure for obtaining GAMMA (only after first call).
160 * IF (IDIAG.EQ.1) THEN
161 * GAMMA = AT(1)*R23 + AT(2)*R13 + AT(3)*R12 +
162 * & (R(3)*T12*W1W2 + T23*R(1)*W2W3) - R123*UNRE
163 * IDIAG = 0
164 * NB! IDIAG must be in COMMON and set > 0 by RCHAIN.
165 * END IF
166 *
167  IF (gtype) THEN
168 * Use modified time transformation for critical case.
169  gamma = at(1)*r23 + at(2)*r13 + at(3)*r12 +
170  & r(3)*t12*w1w2 + t23*r(1)*w2w3 - r123*unre
171  sigin = 1.d0/(r12 + r23 + r13)
172  gsigin = 2.d0*gamma*sigin
173  sumr = r(1) + r(2) + r(3)
174  DO 15 i = 1,3
175  si = (sumr - r(i))*gsigin
176  kq = 4*(i - 1)
177  DO 12 k = kq+1,kq+4
178  dq(k) = sigin*dq(k)
179  dp(k) = sigin*(dp(k) + si*q(k))
180  12 CONTINUE
181  15 CONTINUE
182  dq(13) = dq(13)*sigin
183  IF (itfac.NE.0) THEN
184  tfac = tfac*sigin
185  itfac = 0
186  END IF
187  END IF
188 *
189 * Check osculating pericentre of closest pair (first call only).
190  IF (icall.EQ.0) go to 50
191 *
192 * Examine all close pair indices K1 & K2 in chain vector /ICONF/.
193  qpmin = 100.0
194  im0 = 1
195  DO 20 im = 1,3
196  IF (im.EQ.1) THEN
197  k1 = i1
198  k2 = i2
199  ELSE IF (im.EQ.2) THEN
200  k1 = i2
201  k2 = i3
202  ELSE
203  k1 = i3
204  k2 = i4
205  END IF
206 *
207 * Set distance of nearest perturber.
208  IF (im.EQ.1.OR.im.EQ.3) THEN
209  rp = r(2)
210  ELSE
211  rp = min(r(1),r(3))
212  END IF
213 *
214 * Obtain pericentre for small perturbations (ignore mass effect).
215  gi = (r(im)/rp)**3
216  IF (gi.LT.0.005) THEN
217  k = 4*(im - 1) + 1
218  CALL peri(q(k),dq(k),dq(13),m(k1),m(k2),qperi)
219  ELSE
220  qperi = r(im)
221  END IF
222 *
223 * Compare pericentre with previous mutual distances.
224  rij(k1,k2) = min(rij(k1,k2),qperi)
225  rij(k2,k1) = min(rij(k2,k1),qperi)
226 *
227 * Save indices for smallest pericentre and switch off indicator.
228  IF (im.EQ.1.OR.r(im).LE.r(im0)) THEN
229  qpmin = qperi
230  rp0 = rp
231  k10 = k1
232  k20 = k2
233  im0 = im
234  icall = 0
235  END IF
236  20 CONTINUE
237 *
238 * Check smallest pericentre and reset corresponding indices.
239  rcoll = min(rcoll,qpmin)
240  qperi = qpmin
241  im = im0
242  k = 4*(im - 1) + 1
243  k1 = k10
244  k2 = k20
245 *
246 * Save minimum configuration (for EREL4 & TRANS4).
247  DO 25 j = 1,12
248  qk(j) = q(j)
249  pk(j) = p(j)
250  25 CONTINUE
251 *
252 * Check for tidal two-body interaction or stellar collision.
253  IF (qpmin.LT.4.0*max(SIZE(k1),SIZE(k2))) THEN
254 *
255 * Exit if collision distance is not the smallest.
256  IF (r(im).GT.rp0) go to 50
257 *
258 * Convert Q' to standard KS with T' = R and form radial velocity R'.
259  rpr = 0.0d0
260  DO 30 j = k,k+3
261  upr(j) = dq(j)*r(im)/dq(13)
262  rpr = rpr + 2.0d0*q(j)*upr(j)
263  30 CONTINUE
264 *
265 * Determine small semi-major axis from non-singular expressions.
266  CALL erel4(im,eb,semi)
267 *
268 * Obtain pericentre time interval from elements & Kepler's equation.
269  mb = m(k1) + m(k2)
270  CALL tperi(semi,q(k),upr(k),mb,tp)
271 *
272 * Activate collision indicator & B-S step selector (first time).
273  IF (icoll.EQ.0) THEN
274  icoll = -1
275  jc = 1
276  ELSE
277 *
278 * Check convergence: radial velocity < 1.0E-06 parabolic velocity.
279  IF (abs(rpr).LT.1.0e-06*sqrt(2.0d0*mb*r(im)).OR.
280  & (abs(dsc).LT.1.0e-06)) THEN
281 * Reset B-S step selector and set collision indicator to CHAIN index.
282  jc = 0
283  icoll = im
284  END IF
285  END IF
286 *
287 * Set regularized step for DIFSY4 using accelerated iteration rate.
288  IF (r(im).GT.2.0*qpmin) THEN
289  dsc = 2.0*abs(tp)/dq(13)
290  ELSE
291  dsc = abs(tp)/dq(13)
292  END IF
293 *
294 * Ensure negative step beyond pericentre (restore step at end).
295  IF (jc.GT.0.AND.rpr.GT.0.0d0) THEN
296  dsc = -dsc
297  END IF
298  IF (jc.EQ.0) dsc = 1.0
299 *
300 * Begin iteration if R(IM) is smallest or wide pair is past pericentre.
301  IF (r(im).LT.rp0.OR.rpr.GT.0.0d0) THEN
302  dsfac = 0.2*abs(ds/dsc)
303  IF (r(im).LT.rp0) dsfac = 1.0
304  ds = dsfac*dsc
305  ELSE
306 * Switch off indicators and carry on normally until next check.
307  icoll = 0
308  jc = 0
309  END IF
310  END IF
311 *
312  50 RETURN
313 *
314  END