Nbody6
kscorr.f
Go to the documentation of this file.
1  SUBROUTINE kscorr(IPAIR,UI,UIDOT,FP,FD,TD2,TDOT4,TDOT5,TDOT6)
2 *
3 *
4 * Corrector for KS motion.
5 * ------------------------
6 *
7  include 'common6.h'
8  common/slow0/ range,islow(10)
9  REAL*8 ui(4),uidot(4),fp(6),fd(6),freg(4),frd(4),a1(3,4),a(8),
10  & pr(4),prd(4),u2(4),u3(4),u4(4),u5(4),fd1(3)
11  parameter(itlow=2)
12 *
13 *
14 * Convert from physical to regularized derivative using T' = R.
15  ri = r(ipair)
16  DO 1 k = 1,3
17  fd1(k) = fd(k)
18  fd(k) = ri*fd(k)
19  1 CONTINUE
20 *
21 * Include KS slow-down factor in the perturbation if ZMOD > 1.
22  IF (kz(26).GT.0) THEN
23  imod = kslow(ipair)
24  IF (imod.GT.1) THEN
25  zmod = float(islow(imod))
26  DO 5 k = 1,3
27  fp(k) = zmod*fp(k)
28  fd(k) = zmod*fd(k)
29  fd1(k) = zmod*fd1(k)
30  5 CONTINUE
31  END IF
32  END IF
33 *
34 * Predict DH from Taylor series derivatives and save temporary value.
35  dtu = dtau(ipair)
36  dt3 = one6*dtu
37  dt4 = 0.25d0*dtu
38  dh = (((hdot4(ipair)*dt4 + hdot3(ipair))*dt3 +
39  & 0.5d0*hdot2(ipair))*dtu + hdot(ipair))*dtu
40 *
41 * Set time-step factors and copy Stumpff coefficients to scalars.
42  dtsq = dtu**2
43  dt6 = 6.0/(dtu*dtsq)
44  dt2 = 2.0/dtsq
45  z3 = sf(3,ipair)
46  z4 = sf(4,ipair)
47  z5 = sf(5,ipair)
48 *
49 * Define optimized scalars outside iteration loop.
50  dt5 = 0.2d0*dtu*z5
51  dtz = dt4*z4
52  dh0 = (0.5d0*hdot2(ipair)*dtu + hdot(ipair))*dtu
53 *
54 * Use extra iteration and new KSPERT for large perturbation.
55 * IF (GAMMA(IPAIR).LT.1.0D-20) THEN
56  itmax = itlow
57 * ELSE
58 * ITMAX = ITLOW + 1
59 * I1 = 2*IPAIR - 1
60 * NNB0 = LIST(1,I1)
61 * I = N + IPAIR
62 * BODYIN = 1.0/BODY(I)
63 * END IF
64 *
65 * Perform iteration with or without re-calculating perturbation.
66  DO 40 iter = 1,itmax
67 *
68 * Obtain new transformation matrix.
69  CALL matrix(ui,a1)
70 *
71 * Form twice regularized force and half first derivative of H.
72  hd = 0.0d0
73  td2 = 0.0d0
74  DO 10 k = 1,4
75  a(k) = a1(1,k)*fp(1) + a1(2,k)*fp(2) + a1(3,k)*fp(3)
76  a(k+4) = a1(1,k)*fd(1) + a1(2,k)*fd(2) + a1(3,k)*fd(3)
77  pr(k) = ri*a(k)
78  freg(k) = dh*ui(k) + pr(k)
79  hd = hd + uidot(k)*a(k)
80  td2 = td2 + ui(k)*uidot(k)
81  10 CONTINUE
82 *
83 * Set regularized velocity matrix (Levi-Civita matrix not required).
84  CALL matrix(uidot,a1)
85 *
86 * Include the whole (L*F)' term in explicit derivatives of FU & H'.
87  hd2 = 0.0d0
88  DO 15 k = 1,4
89  ak4 = a(k+4) + a1(1,k)*fp(1) + a1(2,k)*fp(2) +
90  & a1(3,k)*fp(3)
91  hd2 = hd2 + (h(ipair)*ui(k) + freg(k))*a(k) +
92  & 2.0d0*uidot(k)*ak4
93  prd(k) = 0.5d0*ri*ak4 + td2*a(k)
94  frd(k) = hd*ui(k) + 0.5d0*dh*uidot(k) + prd(k)
95 * Form the regularized perturbation and modified force.
96  pr(k) = 0.5d0*pr(k)
97  freg(k) = 0.5d0*freg(k)
98  15 CONTINUE
99 *
100 * Determine new derivatives of U evaluated at beginning of the step.
101  DO 20 k = 1,4
102  df = fp0(k,ipair) - freg(k)
103  sum = fd0(k,ipair) + frd(k)
104  bt2 = -3.0d0*df - (sum + fd0(k,ipair))*dtu
105  at3 = 2.0d0*df + sum*dtu
106 *
107  u2(k) = 0.5d0*h0(ipair)*u0(k,ipair) + fp0(k,ipair)
108  u3(k) = 0.5d0*h0(ipair)*udot(k,ipair) + fd0(k,ipair)
109  u4(k) = 0.5d0*h0(ipair)*u2(k) + bt2*dt2
110  u5(k) = 0.5d0*h0(ipair)*u3(k) + at3*dt6
111  20 CONTINUE
112 *
113 * Improve the solution of U and UDOT.
114  DO 25 k = 1,4
115  ui(k) = ((((u5(k)*dt5 + u4(k)*z4)*dt4 + u3(k))*dt3 +
116  & 0.5*u2(k))*dtu + udot(k,ipair))*dtu + u0(k,ipair)
117  uidot(k) = (((u5(k)*dtz + u4(k)*z3)*dt3 +
118  & 0.5*u3(k))*dtu + u2(k))*dtu + udot(k,ipair)
119  25 CONTINUE
120 *
121 * Update the physical distance for next iteration or final solution.
122  ri = ui(1)**2 + ui(2)**2 + ui(3)**2 + ui(4)**2
123 *
124 * Choose between old and new perturbation (skip last time).
125 * IF (ITER.GT.ITMAX) THEN
126 * IF (GAMMA(IPAIR).LT.1.0D-04) THEN
127 * DO 30 K = 1,3
128 * FD(K) = RI*FD1(K)
129 * 30 CONTINUE
130 * ELSE
131 * Transform to improved coordinates & velocities.
132 * CALL KSTRAN(IPAIR,I1,I,BODYIN,RI,UI,UIDOT,XI,VI)
133 *
134 * Re-calculate the perturbing force & derivative.
135 * CALL KSPERT(I1,NNB0,XI,VI,FP,FD)
136 *
137 * Convert to regularized derivative (note slow-down not active).
138 * DO 35 K = 1,3
139 * FD(K) = RI*FD(K)
140 * 35 CONTINUE
141 * END IF
142 * END IF
143 *
144 * Re-evaluate DH by adding Hermite corrector.
145  hd = 2.0d0*hd
146  dhd = hdot(ipair) - hd
147  sum = hdot2(ipair) + hd2
148  bt2 = -3.0d0*dhd - (sum + hdot2(ipair))*dtu
149  at3 = 2.0d0*dhd + sum*dtu
150  dh = dh0 + (0.25d0*at3 + one3*bt2)*dtu
151  40 CONTINUE
152 *
153 * Copy final values and set higher derivatives.
154  h(ipair) = h(ipair) + dh
155  DO 50 k = 1,4
156  u(k,ipair) = ui(k)
157  u0(k,ipair) = ui(k)
158  udot(k,ipair) = uidot(k)
159  fudot2(k,ipair) = u4(k) + u5(k)*dtu
160  fudot3(k,ipair) = u5(k)
161  ak4 = a(k+4) + a1(1,k)*fp(1) + a1(2,k)*fp(2) + a1(3,k)*fp(3)
162  pr(k) = 0.5d0*ri*a(k)
163  prd(k) = 0.5d0*ri*ak4 + td2*a(k) + 0.5d0*hd*ui(k)
164  fp0(k,ipair) = pr(k)
165  fd0(k,ipair) = prd(k)
166  50 CONTINUE
167 *
168 * Save new derivatives of H.
169  hdot(ipair) = hd
170  hdot2(ipair) = hd2
171  hdot3(ipair) = (3.0d0*at3 + bt2)*dt2
172  hdot4(ipair) = at3*dt6
173 *
174 * Set new FU/2 & FUDOT/6 and form scalar terms for time derivatives.
175  td2 = 0.0d0
176  td3 = 0.0d0
177  tdot4 = 0.0d0
178  tdot5 = 0.0d0
179  tdot6 = 0.0d0
180  h2 = 0.5d0*h(ipair)
181  DO 60 k = 1,4
182  u2k = pr(k) + h2*ui(k)
183  u3k = prd(k) + h2*uidot(k)
184  fu(k,ipair) = 0.5d0*u2k
185  fudot(k,ipair) = one6*u3k
186  td2 = td2 + ui(k)*uidot(k)
187  td3 = td3 + uidot(k)**2 + ui(k)*u2k
188  tdot4 = tdot4 + ui(k)*u3k + 3.0d0*uidot(k)*u2k
189  tdot5 = tdot5 + 0.5d0*fudot2(k,ipair)*ui(k) +
190  & 2.0d0*u3k*uidot(k) + 1.5d0*u2k**2
191  tdot6 = tdot6 + u5(k)*ui(k) + 5.0d0*u4(k)*uidot(k) +
192  & 10.0d0*u3k*u2k
193 * Note that TDOT4/2, TDOT5/4 and TDOT6/2 are calculated.
194  60 CONTINUE
195 *
196 * Save distance and second & third time derivatives.
197  r(ipair) = ri
198  tdot2(ipair) = 2.0d0*td2
199  tdot3(ipair) = 2.0d0*td3
200 *
201  RETURN
202 *
203  END