| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4DormandPrince745 implementation 26 // G4DormandPrince745 implementation
27 // 27 //
28 // DormandPrince7 - 5(4) non-FSAL 28 // DormandPrince7 - 5(4) non-FSAL
29 // definition of the stepper() method that eva 29 // definition of the stepper() method that evaluates one step in
30 // field propagation. 30 // field propagation.
31 // The coefficients and the algorithm have bee 31 // The coefficients and the algorithm have been adapted from
32 // 32 //
33 // J. R. Dormand and P. J. Prince, "A family o 33 // J. R. Dormand and P. J. Prince, "A family of embedded Runge-Kutta formulae"
34 // Journal of computational and applied Math., 34 // Journal of computational and applied Math., vol.6, no.1, pp.19-26, 1980.
35 // 35 //
36 // The Butcher table of the Dormand-Prince-7-4 36 // The Butcher table of the Dormand-Prince-7-4-5 method is as follows :
37 // 37 //
38 // 0 | 38 // 0 |
39 // 1/5 | 1/5 39 // 1/5 | 1/5
40 // 3/10| 3/40 9/40 40 // 3/10| 3/40 9/40
41 // 4/5 | 44/45 56/15 32/9 41 // 4/5 | 44/45 56/15 32/9
42 // 8/9 | 19372/6561 25360/2187 64448/6561 42 // 8/9 | 19372/6561 25360/2187 64448/6561 212/729
43 // 1 | 9017/3168 355/33 46732/5247 43 // 1 | 9017/3168 355/33 46732/5247 49/176 5103/18656
44 // 1 | 35/384 0 500/1113 44 // 1 | 35/384 0 500/1113 125/192 2187/6784 11/84
45 // ---------------------------------------- 45 // ------------------------------------------------------------------------
46 // 35/384 0 500/1113 46 // 35/384 0 500/1113 125/192 2187/6784 11/84 0
47 // 5179/57600 0 7571/16695 47 // 5179/57600 0 7571/16695 393/640 92097/339200 187/2100 1/40
48 // 48 //
49 // Created: Somnath Banerjee, Google Summer of 49 // Created: Somnath Banerjee, Google Summer of Code 2015, 25 May 2015
50 // Supervision: John Apostolakis, CERN 50 // Supervision: John Apostolakis, CERN
51 // ------------------------------------------- 51 // --------------------------------------------------------------------
52 52
53 #include "G4DormandPrince745.hh" 53 #include "G4DormandPrince745.hh"
54 #include "G4LineSection.hh" 54 #include "G4LineSection.hh"
55 55
56 #include <cstring> 56 #include <cstring>
57 57
58 using namespace field_utils; 58 using namespace field_utils;
59 59
60 // Name of this steppers <<
61 const G4String& G4DormandPrince745::StepperTyp <<
62 { <<
63 static G4String _stepperType("G4DormandPrinc <<
64 return _stepperType; <<
65 } <<
66 <<
67 // Description of this steppers - plus details <<
68 const G4String& G4DormandPrince745::StepperDes <<
69 { <<
70 static G4String _stepperDescription( <<
71 "Embedeed 5th order Runge-Kutta stepper - <<
72 return _stepperDescription; <<
73 } <<
74 <<
75 G4DormandPrince745::G4DormandPrince745(G4Equat 60 G4DormandPrince745::G4DormandPrince745(G4EquationOfMotion* equation,
76 G4int n 61 G4int noIntegrationVariables)
77 : G4MagIntegratorStepper(equation, noInteg 62 : G4MagIntegratorStepper(equation, noIntegrationVariables)
78 { 63 {
79 } 64 }
80 65
81 void G4DormandPrince745::Stepper(const G4doubl 66 void G4DormandPrince745::Stepper(const G4double yInput[],
82 const G4doubl 67 const G4double dydx[],
83 G4doubl 68 G4double hstep,
84 G4doubl 69 G4double yOutput[],
85 G4doubl 70 G4double yError[],
86 G4doubl 71 G4double dydxOutput[])
87 { 72 {
88 Stepper(yInput, dydx, hstep, yOutput, yError 73 Stepper(yInput, dydx, hstep, yOutput, yError);
89 field_utils::copy(dydxOutput, ak7); << 74 copy(dydxOutput, ak7);
90 } 75 }
91 76
92 // Stepper 77 // Stepper
93 // 78 //
94 // Passing in the value of yInput[],the first 79 // Passing in the value of yInput[],the first time dydx[] and Step length
95 // Giving back yOut and yErr arrays for output 80 // Giving back yOut and yErr arrays for output and error respectively
96 // 81 //
97 void G4DormandPrince745::Stepper(const G4doubl 82 void G4DormandPrince745::Stepper(const G4double yInput[],
98 const G4doubl 83 const G4double dydx[],
99 G4doubl 84 G4double hstep,
100 G4doubl 85 G4double yOut[],
101 G4doubl 86 G4double yErr[])
102 { 87 {
103 // The various constants defined on the ba 88 // The various constants defined on the basis of butcher tableu
104 // 89 //
105 const G4double b21 = 0.2, 90 const G4double b21 = 0.2,
106 b31 = 3.0 / 40.0, b32 = 9.0 91 b31 = 3.0 / 40.0, b32 = 9.0 / 40.0,
107 b41 = 44.0 / 45.0, b42 = -5 92 b41 = 44.0 / 45.0, b42 = -56.0 / 15.0, b43 = 32.0/9.0,
108 93
109 b51 = 19372.0 / 6561.0, b52 = -25360.0 94 b51 = 19372.0 / 6561.0, b52 = -25360.0 / 2187.0, b53 = 64448.0 / 6561.0,
110 b54 = -212.0 / 729.0, 95 b54 = -212.0 / 729.0,
111 96
112 b61 = 9017.0 / 3168.0 , b62 = -355.0 97 b61 = 9017.0 / 3168.0 , b62 = -355.0 / 33.0,
113 b63 = 46732.0 / 5247.0, b64 = 49.0 / 1 98 b63 = 46732.0 / 5247.0, b64 = 49.0 / 176.0,
114 b65 = -5103.0 / 18656.0, 99 b65 = -5103.0 / 18656.0,
115 100
116 b71 = 35.0 / 384.0, b72 = 0., 101 b71 = 35.0 / 384.0, b72 = 0.,
117 b73 = 500.0 / 1113.0, b74 = 125.0 / 19 102 b73 = 500.0 / 1113.0, b74 = 125.0 / 192.0,
118 b75 = -2187.0 / 6784.0, b76 = 11.0 / 8 103 b75 = -2187.0 / 6784.0, b76 = 11.0 / 84.0,
119 104
120 //Sum of columns, sum(bij) = ei 105 //Sum of columns, sum(bij) = ei
121 // e1 = 0. , 106 // e1 = 0. ,
122 // e2 = 1.0/5.0 , 107 // e2 = 1.0/5.0 ,
123 // e3 = 3.0/10.0 , 108 // e3 = 3.0/10.0 ,
124 // e4 = 4.0/5.0 , 109 // e4 = 4.0/5.0 ,
125 // e5 = 8.0/9.0 , 110 // e5 = 8.0/9.0 ,
126 // e6 = 1.0 , 111 // e6 = 1.0 ,
127 // e7 = 1.0 , 112 // e7 = 1.0 ,
128 113
129 // Difference between the higher and the l 114 // Difference between the higher and the lower order method coeff. :
130 // b7j are the coefficients of higher orde 115 // b7j are the coefficients of higher order
131 116
132 dc1 = -(b71 - 5179.0 / 57600.0), 117 dc1 = -(b71 - 5179.0 / 57600.0),
133 dc2 = -(b72 - .0), 118 dc2 = -(b72 - .0),
134 dc3 = -(b73 - 7571.0 / 16695.0), 119 dc3 = -(b73 - 7571.0 / 16695.0),
135 dc4 = -(b74 - 393.0 / 640.0), 120 dc4 = -(b74 - 393.0 / 640.0),
136 dc5 = -(b75 + 92097.0 / 339200.0), 121 dc5 = -(b75 + 92097.0 / 339200.0),
137 dc6 = -(b76 - 187.0 / 2100.0), 122 dc6 = -(b76 - 187.0 / 2100.0),
138 dc7 = -(- 1.0 / 40.0); 123 dc7 = -(- 1.0 / 40.0);
139 124
140 const G4int numberOfVariables = GetNumberO 125 const G4int numberOfVariables = GetNumberOfVariables();
141 State yTemp = {0., 0., 0., 0., 0., 0., 0., << 126 State yTemp;
142 127
143 // The number of variables to be integrate 128 // The number of variables to be integrated over
144 // 129 //
145 yOut[7] = yTemp[7] = yInput[7]; 130 yOut[7] = yTemp[7] = yInput[7];
146 131
147 // Saving yInput because yInput and yOut 132 // Saving yInput because yInput and yOut can be aliases for same array
148 // 133 //
149 for(G4int i = 0; i < numberOfVariables; ++ 134 for(G4int i = 0; i < numberOfVariables; ++i)
150 { 135 {
151 fyIn[i] = yInput[i]; 136 fyIn[i] = yInput[i];
152 } 137 }
153 // RightHandSide(yIn, dydx); // Not don 138 // RightHandSide(yIn, dydx); // Not done! 1st stage
154 139
155 for(G4int i = 0; i < numberOfVariables; ++ 140 for(G4int i = 0; i < numberOfVariables; ++i)
156 { 141 {
157 yTemp[i] = fyIn[i] + b21 * hstep * dyd 142 yTemp[i] = fyIn[i] + b21 * hstep * dydx[i];
158 } 143 }
159 RightHandSide(yTemp, ak2); // 144 RightHandSide(yTemp, ak2); // 2nd stage
160 145
161 for(G4int i = 0; i < numberOfVariables; ++ 146 for(G4int i = 0; i < numberOfVariables; ++i)
162 { 147 {
163 yTemp[i] = fyIn[i] + hstep * (b31 * dy 148 yTemp[i] = fyIn[i] + hstep * (b31 * dydx[i] + b32 * ak2[i]);
164 } 149 }
165 RightHandSide(yTemp, ak3); // 150 RightHandSide(yTemp, ak3); // 3rd stage
166 151
167 for(G4int i = 0; i < numberOfVariables; ++ 152 for(G4int i = 0; i < numberOfVariables; ++i)
168 { 153 {
169 yTemp[i] = fyIn[i] + hstep * ( 154 yTemp[i] = fyIn[i] + hstep * (
170 b41 * dydx[i] + b42 * ak2[i] + b43 155 b41 * dydx[i] + b42 * ak2[i] + b43 * ak3[i]);
171 } 156 }
172 RightHandSide(yTemp, ak4); // 157 RightHandSide(yTemp, ak4); // 4th stage
173 158
174 for(G4int i = 0; i < numberOfVariables; ++ 159 for(G4int i = 0; i < numberOfVariables; ++i)
175 { 160 {
176 yTemp[i] = fyIn[i] + hstep * ( 161 yTemp[i] = fyIn[i] + hstep * (
177 b51 * dydx[i] + b52 * ak2[i] + b53 162 b51 * dydx[i] + b52 * ak2[i] + b53 * ak3[i] + b54 * ak4[i]);
178 } 163 }
179 RightHandSide(yTemp, ak5); // 164 RightHandSide(yTemp, ak5); // 5th stage
180 165
181 for(G4int i = 0; i < numberOfVariables; ++ 166 for(G4int i = 0; i < numberOfVariables; ++i)
182 { 167 {
183 yTemp[i] = fyIn[i] + hstep * ( 168 yTemp[i] = fyIn[i] + hstep * (
184 b61 * dydx[i] + b62 * ak2[i] + 169 b61 * dydx[i] + b62 * ak2[i] +
185 b63 * ak3[i] + b64 * ak4[i] + b65 170 b63 * ak3[i] + b64 * ak4[i] + b65 * ak5[i]);
186 } 171 }
187 RightHandSide(yTemp, ak6); // 172 RightHandSide(yTemp, ak6); // 6th stage
188 173
189 for(G4int i = 0; i < numberOfVariables; ++ 174 for(G4int i = 0; i < numberOfVariables; ++i)
190 { 175 {
191 yOut[i] = fyIn[i] + hstep * ( 176 yOut[i] = fyIn[i] + hstep * (
192 b71 * dydx[i] + b72 * ak2[i] + b73 177 b71 * dydx[i] + b72 * ak2[i] + b73 * ak3[i] +
193 b74 * ak4[i] + b75 * ak5[i] + b76 178 b74 * ak4[i] + b75 * ak5[i] + b76 * ak6[i]);
194 } 179 }
195 RightHandSide(yOut, ak7); // 180 RightHandSide(yOut, ak7); // 7th and Final stage
196 181
197 for(G4int i = 0; i < numberOfVariables; ++ 182 for(G4int i = 0; i < numberOfVariables; ++i)
198 { 183 {
199 yErr[i] = hstep * ( 184 yErr[i] = hstep * (
200 dc1 * dydx[i] + dc2 * ak2[i] + 185 dc1 * dydx[i] + dc2 * ak2[i] +
201 dc3 * ak3[i] + dc4 * ak4[i] + 186 dc3 * ak3[i] + dc4 * ak4[i] +
202 dc5 * ak5[i] + dc6 * ak6[i] + dc7 187 dc5 * ak5[i] + dc6 * ak6[i] + dc7 * ak7[i]
203 ) + 1.5e-18; 188 ) + 1.5e-18;
204 189
205 // Store Input and Final values, for p 190 // Store Input and Final values, for possible use in calculating chord
206 // 191 //
207 fyOut[i] = yOut[i]; 192 fyOut[i] = yOut[i];
208 fdydxIn[i] = dydx[i]; 193 fdydxIn[i] = dydx[i];
209 } 194 }
210 195
211 fLastStepLength = hstep; 196 fLastStepLength = hstep;
212 } 197 }
213 198
214 G4double G4DormandPrince745::DistChord() const 199 G4double G4DormandPrince745::DistChord() const
215 { 200 {
216 // Coefficients were taken from Some Pract 201 // Coefficients were taken from Some Practical Runge-Kutta Formulas
217 // by Lawrence F. Shampine, page 149, c* 202 // by Lawrence F. Shampine, page 149, c*
218 // 203 //
219 const G4double hf1 = 6025192743.0 / 300855 204 const G4double hf1 = 6025192743.0 / 30085553152.0,
220 hf3 = 51252292925.0 / 65400 205 hf3 = 51252292925.0 / 65400821598.0,
221 hf4 = - 2691868925.0 / 4512 206 hf4 = - 2691868925.0 / 45128329728.0,
222 hf5 = 187940372067.0 / 1594 207 hf5 = 187940372067.0 / 1594534317056.0,
223 hf6 = - 1776094331.0 / 1974 208 hf6 = - 1776094331.0 / 19743644256.0,
224 hf7 = 11237099.0 / 23504338 209 hf7 = 11237099.0 / 235043384.0;
225 210
226 G4ThreeVector mid; 211 G4ThreeVector mid;
227 212
228 for(G4int i = 0; i < 3; ++i) 213 for(G4int i = 0; i < 3; ++i)
229 { 214 {
230 mid[i] = fyIn[i] + 0.5 * fLastStepLeng 215 mid[i] = fyIn[i] + 0.5 * fLastStepLength * (
231 hf1 * fdydxIn[i] + hf3 * ak3[i] + 216 hf1 * fdydxIn[i] + hf3 * ak3[i] +
232 hf4 * ak4[i] + hf5 * ak5[i] + hf6 217 hf4 * ak4[i] + hf5 * ak5[i] + hf6 * ak6[i] + hf7 * ak7[i]);
233 } 218 }
234 219
235 const G4ThreeVector begin = makeVector(fyI 220 const G4ThreeVector begin = makeVector(fyIn, Value3D::Position);
236 const G4ThreeVector end = makeVector(fyOut 221 const G4ThreeVector end = makeVector(fyOut, Value3D::Position);
237 222
238 return G4LineSection::Distline(mid, begin, 223 return G4LineSection::Distline(mid, begin, end);
239 } 224 }
240 225
241 // The lower (4th) order interpolant given by 226 // The lower (4th) order interpolant given by Dormand and Prince:
242 // J. R. Dormand and P. J. Prince, "Run 227 // J. R. Dormand and P. J. Prince, "Runge-Kutta triples"
243 // Computers & Mathematics with Applica 228 // Computers & Mathematics with Applications, vol. 12, no. 9,
244 // pp. 1007-1017, 1986. 229 // pp. 1007-1017, 1986.
245 // 230 //
246 void G4DormandPrince745:: 231 void G4DormandPrince745::
247 Interpolate4thOrder(G4double yOut[], G4double 232 Interpolate4thOrder(G4double yOut[], G4double tau) const
248 { 233 {
249 const G4int numberOfVariables = GetNumberO 234 const G4int numberOfVariables = GetNumberOfVariables();
250 235
251 const G4double tau2 = tau * tau, 236 const G4double tau2 = tau * tau,
252 tau3 = tau * tau2, 237 tau3 = tau * tau2,
253 tau4 = tau2 * tau2; 238 tau4 = tau2 * tau2;
254 239
255 const G4double bf1 = 1.0 / 11282082432.0 * 240 const G4double bf1 = 1.0 / 11282082432.0 * (
256 157015080.0 * tau4 - 13107642775.0 * t 241 157015080.0 * tau4 - 13107642775.0 * tau3 + 34969693132.0 * tau2 -
257 32272833064.0 * tau + 11282082432.0); 242 32272833064.0 * tau + 11282082432.0);
258 243
259 const G4double bf3 = - 100.0 / 32700410799 244 const G4double bf3 = - 100.0 / 32700410799.0 * tau * (
260 15701508.0 * tau3 - 914128567.0 * tau2 245 15701508.0 * tau3 - 914128567.0 * tau2 + 2074956840.0 * tau -
261 1323431896.0); 246 1323431896.0);
262 247
263 const G4double bf4 = 25.0 / 5641041216.0 * 248 const G4double bf4 = 25.0 / 5641041216.0 * tau * (
264 94209048.0 * tau3 - 1518414297.0 * tau 249 94209048.0 * tau3 - 1518414297.0 * tau2 + 2460397220.0 * tau -
265 889289856.0); 250 889289856.0);
266 251
267 const G4double bf5 = - 2187.0 / 1993167896 252 const G4double bf5 = - 2187.0 / 199316789632.0 * tau * (
268 52338360.0 * tau3 - 451824525.0 * tau2 253 52338360.0 * tau3 - 451824525.0 * tau2 + 687873124.0 * tau -
269 259006536.0); 254 259006536.0);
270 255
271 const G4double bf6 = 11.0 / 2467955532.0 * 256 const G4double bf6 = 11.0 / 2467955532.0 * tau * (
272 106151040.0 * tau3 - 661884105.0 * tau 257 106151040.0 * tau3 - 661884105.0 * tau2 +
273 946554244.0 * tau - 361440756.0); 258 946554244.0 * tau - 361440756.0);
274 259
275 const G4double bf7 = 1.0 / 29380423.0 * ta 260 const G4double bf7 = 1.0 / 29380423.0 * tau * (1.0 - tau) * (
276 8293050.0 * tau2 - 82437520.0 * tau + 261 8293050.0 * tau2 - 82437520.0 * tau + 44764047.0);
277 262
278 for(G4int i = 0; i < numberOfVariables; ++ 263 for(G4int i = 0; i < numberOfVariables; ++i)
279 { 264 {
280 yOut[i] = fyIn[i] + fLastStepLength * 265 yOut[i] = fyIn[i] + fLastStepLength * tau * (
281 bf1 * fdydxIn[i] + bf3 * ak3[i] + 266 bf1 * fdydxIn[i] + bf3 * ak3[i] + bf4 * ak4[i] +
282 bf5 * ak5[i] + bf6 * ak6[i] + bf7 267 bf5 * ak5[i] + bf6 * ak6[i] + bf7 * ak7[i]);
283 } 268 }
284 } 269 }
285 270
286 // Following interpolant of order 5 was given 271 // Following interpolant of order 5 was given by Baker,Dormand,Gilmore, Prince :
287 // T. S. Baker, J. R. Dormand, J. P. Gi 272 // T. S. Baker, J. R. Dormand, J. P. Gilmore, and P. J. Prince,
288 // "Continuous approximation with embed 273 // "Continuous approximation with embedded Runge-Kutta methods"
289 // Applied Numerical Mathematics, vol. 274 // Applied Numerical Mathematics, vol. 22, no. 1, pp. 51-62, 1996.
290 // 275 //
291 // Calculating the extra stages for the interp 276 // Calculating the extra stages for the interpolant
292 // 277 //
293 void G4DormandPrince745::SetupInterpolation5th 278 void G4DormandPrince745::SetupInterpolation5thOrder()
294 { 279 {
295 // Coefficients for the additional stages 280 // Coefficients for the additional stages
296 // 281 //
297 const G4double b81 = 6245.0 / 62208.0, 282 const G4double b81 = 6245.0 / 62208.0,
298 b82 = 0.0, 283 b82 = 0.0,
299 b83 = 8875.0 / 103032.0, 284 b83 = 8875.0 / 103032.0,
300 b84 = -125.0 / 1728.0, 285 b84 = -125.0 / 1728.0,
301 b85 = 801.0 / 13568.0, 286 b85 = 801.0 / 13568.0,
302 b86 = -13519.0 / 368064.0, 287 b86 = -13519.0 / 368064.0,
303 b87 = 11105.0 / 368064.0, 288 b87 = 11105.0 / 368064.0,
304 289
305 b91 = 632855.0 / 4478976.0, 290 b91 = 632855.0 / 4478976.0,
306 b92 = 0.0, 291 b92 = 0.0,
307 b93 = 4146875.0 / 6491016.0 292 b93 = 4146875.0 / 6491016.0,
308 b94 = 5490625.0 /14183424.0 293 b94 = 5490625.0 /14183424.0,
309 b95 = -15975.0 / 108544.0, 294 b95 = -15975.0 / 108544.0,
310 b96 = 8295925.0 / 220286304 295 b96 = 8295925.0 / 220286304.0,
311 b97 = -1779595.0 / 62938944 296 b97 = -1779595.0 / 62938944.0,
312 b98 = -805.0 / 4104.0; 297 b98 = -805.0 / 4104.0;
313 298
314 const G4int numberOfVariables = GetNumberO 299 const G4int numberOfVariables = GetNumberOfVariables();
315 State yTemp = {0., 0., 0., 0., 0., 0., 0., << 300 State yTemp;
316 301
317 // Evaluate the extra stages 302 // Evaluate the extra stages
318 // 303 //
319 for(G4int i = 0; i < numberOfVariables; ++ 304 for(G4int i = 0; i < numberOfVariables; ++i)
320 { 305 {
321 yTemp[i] = fyIn[i] + fLastStepLength * 306 yTemp[i] = fyIn[i] + fLastStepLength * (
322 b81 * fdydxIn[i] + b82 * ak2[i] + 307 b81 * fdydxIn[i] + b82 * ak2[i] + b83 * ak3[i] +
323 b84 * ak4[i] + b85 * ak5[i] + b86 308 b84 * ak4[i] + b85 * ak5[i] + b86 * ak6[i] +
324 b87 * ak7[i] 309 b87 * ak7[i]
325 ); 310 );
326 } 311 }
327 RightHandSide(yTemp, ak8); // 8th 312 RightHandSide(yTemp, ak8); // 8th Stage
328 313
329 for(G4int i = 0; i < numberOfVariables; ++ 314 for(G4int i = 0; i < numberOfVariables; ++i)
330 { 315 {
331 yTemp[i] = fyIn[i] + fLastStepLength * 316 yTemp[i] = fyIn[i] + fLastStepLength * (
332 b91 * fdydxIn[i] + b92 * ak2[i] + 317 b91 * fdydxIn[i] + b92 * ak2[i] + b93 * ak3[i] +
333 b94 * ak4[i] + b95 * ak5[i] + b96 318 b94 * ak4[i] + b95 * ak5[i] + b96 * ak6[i] +
334 b97 * ak7[i] + b98 * ak8[i] 319 b97 * ak7[i] + b98 * ak8[i]
335 ); 320 );
336 } 321 }
337 RightHandSide(yTemp, ak9); // 9th 322 RightHandSide(yTemp, ak9); // 9th Stage
338 } 323 }
339 324
340 // Calculating the interpolated result yOut wi 325 // Calculating the interpolated result yOut with the coefficients
341 // 326 //
342 void G4DormandPrince745:: 327 void G4DormandPrince745::
343 Interpolate5thOrder(G4double yOut[], G4double 328 Interpolate5thOrder(G4double yOut[], G4double tau) const
344 { 329 {
345 // Define the coefficients for the polynom 330 // Define the coefficients for the polynomials
346 // 331 //
347 G4double bi[10][5]; 332 G4double bi[10][5];
348 333
349 // COEFFICIENTS OF bi[1] 334 // COEFFICIENTS OF bi[1]
350 bi[1][0] = 1.0, 335 bi[1][0] = 1.0,
351 bi[1][1] = -38039.0 / 7040.0, 336 bi[1][1] = -38039.0 / 7040.0,
352 bi[1][2] = 125923.0 / 10560.0, 337 bi[1][2] = 125923.0 / 10560.0,
353 bi[1][3] = -19683.0 / 1760.0, 338 bi[1][3] = -19683.0 / 1760.0,
354 bi[1][4] = 3303.0 / 880.0, 339 bi[1][4] = 3303.0 / 880.0,
355 // -------------------------------------- 340 // --------------------------------------------------------
356 // 341 //
357 // COEFFICIENTS OF bi[2] 342 // COEFFICIENTS OF bi[2]
358 bi[2][0] = 0.0, 343 bi[2][0] = 0.0,
359 bi[2][1] = 0.0, 344 bi[2][1] = 0.0,
360 bi[2][2] = 0.0, 345 bi[2][2] = 0.0,
361 bi[2][3] = 0.0, 346 bi[2][3] = 0.0,
362 bi[2][4] = 0.0, 347 bi[2][4] = 0.0,
363 // -------------------------------------- 348 // --------------------------------------------------------
364 // 349 //
365 // COEFFICIENTS OF bi[3] 350 // COEFFICIENTS OF bi[3]
366 bi[3][0] = 0.0, 351 bi[3][0] = 0.0,
367 bi[3][1] = -12500.0 / 4081.0, 352 bi[3][1] = -12500.0 / 4081.0,
368 bi[3][2] = 205000.0 / 12243.0, 353 bi[3][2] = 205000.0 / 12243.0,
369 bi[3][3] = -90000.0 / 4081.0, 354 bi[3][3] = -90000.0 / 4081.0,
370 bi[3][4] = 36000.0 / 4081.0, 355 bi[3][4] = 36000.0 / 4081.0,
371 // -------------------------------------- 356 // --------------------------------------------------------
372 // 357 //
373 // COEFFICIENTS OF bi[4] 358 // COEFFICIENTS OF bi[4]
374 bi[4][0] = 0.0, 359 bi[4][0] = 0.0,
375 bi[4][1] = -3125.0 / 704.0, 360 bi[4][1] = -3125.0 / 704.0,
376 bi[4][2] = 25625.0 / 1056.0, 361 bi[4][2] = 25625.0 / 1056.0,
377 bi[4][3] = -5625.0 / 176.0, 362 bi[4][3] = -5625.0 / 176.0,
378 bi[4][4] = 1125.0 / 88.0, 363 bi[4][4] = 1125.0 / 88.0,
379 // -------------------------------------- 364 // --------------------------------------------------------
380 // 365 //
381 // COEFFICIENTS OF bi[5] 366 // COEFFICIENTS OF bi[5]
382 bi[5][0] = 0.0, 367 bi[5][0] = 0.0,
383 bi[5][1] = 164025.0 / 74624.0, 368 bi[5][1] = 164025.0 / 74624.0,
384 bi[5][2] = -448335.0 / 37312.0, 369 bi[5][2] = -448335.0 / 37312.0,
385 bi[5][3] = 295245.0 / 18656.0, 370 bi[5][3] = 295245.0 / 18656.0,
386 bi[5][4] = -59049.0 / 9328.0, 371 bi[5][4] = -59049.0 / 9328.0,
387 // -------------------------------------- 372 // --------------------------------------------------------
388 // 373 //
389 // COEFFICIENTS OF bi[6] 374 // COEFFICIENTS OF bi[6]
390 bi[6][0] = 0.0, 375 bi[6][0] = 0.0,
391 bi[6][1] = -25.0 / 28.0, 376 bi[6][1] = -25.0 / 28.0,
392 bi[6][2] = 205.0 / 42.0, 377 bi[6][2] = 205.0 / 42.0,
393 bi[6][3] = -45.0 / 7.0, 378 bi[6][3] = -45.0 / 7.0,
394 bi[6][4] = 18.0 / 7.0, 379 bi[6][4] = 18.0 / 7.0,
395 // -------------------------------------- 380 // --------------------------------------------------------
396 // 381 //
397 // COEFFICIENTS OF bi[7] 382 // COEFFICIENTS OF bi[7]
398 bi[7][0] = 0.0, 383 bi[7][0] = 0.0,
399 bi[7][1] = -2.0 / 11.0, 384 bi[7][1] = -2.0 / 11.0,
400 bi[7][2] = 73.0 / 55.0, 385 bi[7][2] = 73.0 / 55.0,
401 bi[7][3] = -171.0 / 55.0, 386 bi[7][3] = -171.0 / 55.0,
402 bi[7][4] = 108.0 / 55.0, 387 bi[7][4] = 108.0 / 55.0,
403 // -------------------------------------- 388 // --------------------------------------------------------
404 // 389 //
405 // COEFFICIENTS OF bi[8] 390 // COEFFICIENTS OF bi[8]
406 bi[8][0] = 0.0, 391 bi[8][0] = 0.0,
407 bi[8][1] = 189.0 / 22.0, 392 bi[8][1] = 189.0 / 22.0,
408 bi[8][2] = -1593.0 / 55.0, 393 bi[8][2] = -1593.0 / 55.0,
409 bi[8][3] = 3537.0 / 110.0, 394 bi[8][3] = 3537.0 / 110.0,
410 bi[8][4] = -648.0 / 55.0, 395 bi[8][4] = -648.0 / 55.0,
411 // -------------------------------------- 396 // --------------------------------------------------------
412 // 397 //
413 // COEFFICIENTS OF bi[9] 398 // COEFFICIENTS OF bi[9]
414 bi[9][0] = 0.0, 399 bi[9][0] = 0.0,
415 bi[9][1] = 351.0 / 110.0, 400 bi[9][1] = 351.0 / 110.0,
416 bi[9][2] = -999.0 / 55.0, 401 bi[9][2] = -999.0 / 55.0,
417 bi[9][3] = 2943.0 / 110.0, 402 bi[9][3] = 2943.0 / 110.0,
418 bi[9][4] = -648.0 / 55.0; 403 bi[9][4] = -648.0 / 55.0;
419 // -------------------------------------- 404 // --------------------------------------------------------
420 405
421 // Calculating the polynomials 406 // Calculating the polynomials
422 407
423 G4double b[10]; 408 G4double b[10];
424 std::memset(b, 0.0, sizeof(b)); 409 std::memset(b, 0.0, sizeof(b));
425 410
426 G4double tauPower = 1.0; 411 G4double tauPower = 1.0;
427 for(G4int j = 0; j <= 4; ++j) 412 for(G4int j = 0; j <= 4; ++j)
428 { 413 {
429 for(G4int iStage = 1; iStage <= 9; ++iS 414 for(G4int iStage = 1; iStage <= 9; ++iStage)
430 { 415 {
431 b[iStage] += bi[iStage][j] * tauPowe 416 b[iStage] += bi[iStage][j] * tauPower;
432 } 417 }
433 tauPower *= tau; 418 tauPower *= tau;
434 } 419 }
435 420
436 const G4int numberOfVariables = GetNumberO 421 const G4int numberOfVariables = GetNumberOfVariables();
437 const G4double stepLen = fLastStepLength * 422 const G4double stepLen = fLastStepLength * tau;
438 for(G4int i = 0; i < numberOfVariables; ++ 423 for(G4int i = 0; i < numberOfVariables; ++i)
439 { 424 {
440 yOut[i] = fyIn[i] + stepLen * ( 425 yOut[i] = fyIn[i] + stepLen * (
441 b[1] * fdydxIn[i] + b[2] * ak2[i] 426 b[1] * fdydxIn[i] + b[2] * ak2[i] + b[3] * ak3[i] +
442 b[4] * ak4[i] + b[5] * ak5[i] + b[ 427 b[4] * ak4[i] + b[5] * ak5[i] + b[6] * ak6[i] +
443 b[7] * ak7[i] + b[8] * ak8[i] + b[ 428 b[7] * ak7[i] + b[8] * ak8[i] + b[9] * ak9[i]
444 ); 429 );
445 } 430 }
446 } 431 }
447 432