Nbody6
 All Files Functions Variables
derqp3.f
Go to the documentation of this file.
1  SUBROUTINE derqp3(Q,P,QPR,PPR,TPR)
2 *
3 *
4 * Equations of motion for AZ regularization.
5 * ------------------------------------------
6 *
7  IMPLICIT REAL*8 (a-h,m,o-z)
8  common/azreg/ time3,tmax,v(8),w(8),r1,r2,r3,energy,m(3),x3(3,3),
9  & xdot3(3,3),cm(10),c11,c12,c19,c20,c24,c25,
10  & nstep3,name3(3),kz15,kz27
11 * Note COMMON locations of Q & P are replaced by dummy variables.
12  common/close/ rij(4,4),rcoll,qperi,SIZE(4),ecoll3,ip(4)
13  common/azcoll/ rk(3),qk(8),pk(8),icall,icoll,ndiss3
14  common/bssave/ ep(4),dsc,facm,tfac,itfac,jc
15  common/azconf/ iconf(3)
16  REAL*8 q(8),p(8),qpr(8),ppr(8),s2(8),s5(8),s8(8),upr(8)
17 *
18 *
19 * Form scalar distances & coefficients.
20  r1 = q(1)*q(1) + q(2)*q(2) + q(3)*q(3) + q(4)*q(4)
21  r2 = q(5)*q(5) + q(6)*q(6) + q(7)*q(7) + q(8)*q(8)
22  c3 = q(1)*p(1) - q(2)*p(2) - q(3)*p(3) + q(4)*p(4)
23  c4 = q(5)*p(5) - q(6)*p(6) - q(7)*p(7) + q(8)*p(8)
24  c5 = q(2)*p(1) + q(1)*p(2) - q(4)*p(3) - q(3)*p(4)
25  c6 = q(6)*p(5) + q(5)*p(6) - q(8)*p(7) - q(7)*p(8)
26  c7 = q(3)*p(1) + q(4)*p(2) + q(1)*p(3) + q(2)*p(4)
27  c8 = q(7)*p(5) + q(8)*p(6) + q(5)*p(7) + q(6)*p(8)
28  c9 = p(1)*p(1) + p(2)*p(2) + p(3)*p(3) + p(4)*p(4)
29  c10 = p(5)*p(5) + p(6)*p(6) + p(7)*p(7) + p(8)*p(8)
30  c13 = c11*r2
31  c14 = c12*r1
32  c15 = c12*c10
33  c16 = c11*c9
34  c17 = r2*energy
35  c18 = r1*energy
36 * Note that twice the energy is saved in COMMON.
37  c21 = q(1)*q(1) - q(2)*q(2) - q(3)*q(3) + q(4)*q(4)
38  & - q(5)*q(5) + q(6)*q(6) + q(7)*q(7) - q(8)*q(8)
39  c22 = q(1)*q(2) - q(3)*q(4) - q(5)*q(6) + q(7)*q(8)
40  c23 = q(1)*q(3) + q(2)*q(4) - q(5)*q(7) - q(6)*q(8)
41  c22 = c22 + c22
42  c23 = c23 + c23
43  rr = c21*c21 + c22*c22 + c23*c23
44  r3 = sqrt(rr)
45  a = c25/r3
46 *
47 * Set first derivative of the physical time (standard transformation).
48  tpr = r1*r2
49  b = a*tpr/rr
50 *
51 * Combine vectorial components with relevant coefficients.
52  s2(1) = q(1)*c4 + q(2)*c6 + q(3)*c8
53  s2(2) =-q(2)*c4 + q(1)*c6 + q(4)*c8
54  s2(3) =-q(3)*c4 - q(4)*c6 + q(1)*c8
55  s2(4) = q(4)*c4 - q(3)*c6 + q(2)*c8
56  s2(5) = q(5)*c3 + q(6)*c5 + q(7)*c7
57  s2(6) =-q(6)*c3 + q(5)*c5 + q(8)*c7
58  s2(7) =-q(7)*c3 - q(8)*c5 + q(5)*c7
59  s2(8) = q(8)*c3 - q(7)*c5 + q(6)*c7
60  s5(1) = p(1)*c4 + p(2)*c6 + p(3)*c8
61  s5(2) =-p(2)*c4 + p(1)*c6 + p(4)*c8
62  s5(3) =-p(3)*c4 - p(4)*c6 + p(1)*c8
63  s5(4) = p(4)*c4 - p(3)*c6 + p(2)*c8
64  s5(5) = p(5)*c3 + p(6)*c5 + p(7)*c7
65  s5(6) =-p(6)*c3 + p(5)*c5 + p(8)*c7
66  s5(7) =-p(7)*c3 - p(8)*c5 + p(5)*c7
67  s5(8) = p(8)*c3 - p(7)*c5 + p(6)*c7
68  s8(1) = q(1)*c21 + q(2)*c22 + q(3)*c23
69  s8(2) =-q(2)*c21 + q(1)*c22 + q(4)*c23
70  s8(3) =-q(3)*c21 - q(4)*c22 + q(1)*c23
71  s8(4) = q(4)*c21 - q(3)*c22 + q(2)*c23
72  s8(5) =-q(5)*c21 - q(6)*c22 - q(7)*c23
73  s8(6) = q(6)*c21 - q(5)*c22 - q(8)*c23
74  s8(7) = q(7)*c21 + q(8)*c22 - q(5)*c23
75  s8(8) =-q(8)*c21 + q(7)*c22 - q(6)*c23
76  c1 = c17 - c15 + c19 + a*r2
77  c2 = c18 - c16 + c20 + a*r1
78 *
79 * Form derivatives for standard equations of motion.
80  DO 10 i = 1,4
81  k = i + 4
82  qpr(i) = c13*p(i) + c24*s2(i)
83  qpr(k) = c14*p(k) + c24*s2(k)
84  ppr(i) = c1*q(i) - c24*s5(i) - b*s8(i)
85  ppr(k) = c2*q(k) - c24*s5(k) - b*s8(k)
86  10 CONTINUE
87 *
88 * Check the Bulirsch-Stoer tolerance scaling TFAC = TPR*U (ITFAC > 1).
89  IF (itfac.GT.0) THEN
90  tfac = facm*(r1 + r2 + tpr/r3)
91  itfac = -1
92  END IF
93 *
94 * Evaluate the regularized Hamiltonian (denoted GAMMA in 1974 paper).
95  gamma = 0.5d0*(-r1*c17 + r2*c11*c9 + r1*c12*c10 - c20*r2 - c19*r1
96  & - a*tpr) + c24*(c3*c4 + c5*c6 + c7*c8)
97  g2 = 1.0/sqrt(r1 + r2)
98 * Time transformation of 1976 paper DT/DTAU = R1*R2/(R1 + R2)**0.5.
99 *
100 * ------------------------------------------------------------------
101 * U = M(1)*M(3)*R2 + M(2)*M(3)*R1 + 0.5D0*A*TPR
102 * G2 = 1.0/U
103 * Alternative transformation DT/DTAU = R1*R2/U allows explicit time.
104 * ------------------------------------------------------------------
105 *
106 * Set the explicit time derivative and check tolerance scaling.
107  tpr = tpr*g2
108  IF (itfac.NE.0) THEN
109  tfac = tfac*g2
110  itfac = 0
111  END IF
112 *
113 * Include stabilization terms due to modified time transformation.
114  gamma = gamma*g2
115 * S6 = C19 + A*R2
116  DO 20 i = 1,8
117  qpr(i) = g2*qpr(i)
118  ppr(i) = g2*(ppr(i) + gamma*q(i))
119 * PPR(I) = G2*(PPR(I) + GAMMA*(S6*Q(I) - B*S8(I)))
120 * IF (I.EQ.4) S6 = C20 + A*R1
121  20 CONTINUE
122 *
123 * Check osculating pericentre of closest pair (first call only).
124  IF (icall.EQ.0) go to 50
125 *
126 * Find index & distance of closest pair and set corresponding KS index.
127  IF (r1.LT.r2) THEN
128  im = 1
129  rm = r1
130  ELSE
131  im = 2
132  rm = r2
133  END IF
134  k = 1 + 4*(im - 1)
135 *
136 * Obtain pericentre by Mikkola's algorithm (perturbation < 0.005).
137  gi = 2.0*m(3-im)*(rm/r3)**3/(m(im) + m(3))
138  IF (gi.LT.0.005) THEN
139  CALL peri(q(k),qpr(k),tpr,m(im),m(3),qperi)
140  ELSE
141  qperi = rm
142  END IF
143 *
144 * Identify the close bodies and check mutual distances.
145  i = iconf(im)
146  j = iconf(3)
147  rij(i,j) = min(rij(i,j),qperi)
148  rij(j,i) = min(rij(j,i),qperi)
149 *
150 * Ckeck minimum two-body distance and switch off indicator.
151  rcoll = min(rcoll,qperi)
152  icall = 0
153 *
154 * Check for tidal two-body interaction or stellar collision.
155  IF (qperi.LT.4.0*max(SIZE(im),SIZE(3))) THEN
156 *
157 * Convert QPR to standard KS with T' = R and form radial velocity R'.
158  rpr = 0.0d0
159  DO 25 j = k,k+3
160  upr(j) = qpr(j)*rm/tpr
161  rpr = rpr + 2.0d0*q(j)*upr(j)
162  25 CONTINUE
163 *
164 * Save minimum configuration.
165  DO 30 j = 1,8
166  qk(j) = q(j)
167  pk(j) = p(j)
168  30 CONTINUE
169 *
170 * Determine small semi-major axis from non-singular expressions.
171  CALL erel3(im,eb,semi)
172 *
173 * Obtain pericentre time interval from elements & Kepler's equation.
174  mb = m(im) + m(3)
175  CALL tperi(semi,q(k),upr(k),mb,tp)
176 *
177 * Activate collision indicator & B-S step selector (first time).
178  IF (icoll.EQ.0) THEN
179  icoll = -1
180  jc = 1
181  ELSE
182 *
183 * Check convergence: radial velocity < 1.0E-06 parabolic velocity.
184  IF (abs(rpr).LT.1.0e-08*sqrt(2.0d0*mb*rm).OR.
185  & (abs(dsc).LT.1.0e-08)) THEN
186 * Reset B-S step selector and set collision indicator to AZ pair index.
187  jc = 0
188  icoll = im
189  END IF
190  END IF
191 *
192 * Set regularized step for DIFSY3 using accelerated iteration rate.
193  IF (rm.GT.2.0*qperi) THEN
194  dsc = 2.0*abs(tp)/tpr
195  ELSE
196  dsc = abs(tp)/tpr
197  END IF
198 *
199 * Ensure negative step beyond pericentre (restore step at end).
200  IF (jc.GT.0.AND.rpr.GT.0.0d0) THEN
201  dsc = -dsc
202  END IF
203  IF (jc.EQ.0) dsc = 1.0
204  END IF
205 *
206  50 RETURN
207 *
208  END