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