58 using namespace field_utils;
61 G4int noIntegrationVariables)
73 Stepper(yInput, dydx, hstep, yOutput, yError);
91 b31 = 3.0 / 40.0, b32 = 9.0 / 40.0,
92 b41 = 44.0 / 45.0, b42 = -56.0 / 15.0, b43 = 32.0/9.0,
94 b51 = 19372.0 / 6561.0, b52 = -25360.0 / 2187.0, b53 = 64448.0 / 6561.0,
97 b61 = 9017.0 / 3168.0 , b62 = -355.0 / 33.0,
98 b63 = 46732.0 / 5247.0, b64 = 49.0 / 176.0,
99 b65 = -5103.0 / 18656.0,
101 b71 = 35.0 / 384.0, b72 = 0.,
102 b73 = 500.0 / 1113.0, b74 = 125.0 / 192.0,
103 b75 = -2187.0 / 6784.0, b76 = 11.0 / 84.0,
117 dc1 = -(b71 - 5179.0 / 57600.0),
119 dc3 = -(b73 - 7571.0 / 16695.0),
120 dc4 = -(b74 - 393.0 / 640.0),
121 dc5 = -(b75 + 92097.0 / 339200.0),
122 dc6 = -(b76 - 187.0 / 2100.0),
123 dc7 = -(- 1.0 / 40.0);
130 yOut[7] = yTemp[7] = yInput[7];
134 for(
G4int i = 0; i < numberOfVariables; ++i)
140 for(
G4int i = 0; i < numberOfVariables; ++i)
142 yTemp[i] =
fyIn[i] + b21 * hstep * dydx[i];
146 for(
G4int i = 0; i < numberOfVariables; ++i)
148 yTemp[i] =
fyIn[i] + hstep * (b31 * dydx[i] + b32 *
ak2[i]);
152 for(
G4int i = 0; i < numberOfVariables; ++i)
154 yTemp[i] =
fyIn[i] + hstep * (
155 b41 * dydx[i] + b42 *
ak2[i] + b43 *
ak3[i]);
159 for(
G4int i = 0; i < numberOfVariables; ++i)
161 yTemp[i] =
fyIn[i] + hstep * (
162 b51 * dydx[i] + b52 *
ak2[i] + b53 *
ak3[i] + b54 *
ak4[i]);
166 for(
G4int i = 0; i < numberOfVariables; ++i)
168 yTemp[i] =
fyIn[i] + hstep * (
169 b61 * dydx[i] + b62 *
ak2[i] +
170 b63 *
ak3[i] + b64 *
ak4[i] + b65 *
ak5[i]);
174 for(
G4int i = 0; i < numberOfVariables; ++i)
176 yOut[i] =
fyIn[i] + hstep * (
177 b71 * dydx[i] + b72 *
ak2[i] + b73 *
ak3[i] +
178 b74 *
ak4[i] + b75 *
ak5[i] + b76 *
ak6[i]);
182 for(
G4int i = 0; i < numberOfVariables; ++i)
185 dc1 * dydx[i] + dc2 *
ak2[i] +
186 dc3 *
ak3[i] + dc4 *
ak4[i] +
187 dc5 *
ak5[i] + dc6 *
ak6[i] + dc7 *
ak7[i]
204 const G4double hf1 = 6025192743.0 / 30085553152.0,
205 hf3 = 51252292925.0 / 65400821598.0,
206 hf4 = - 2691868925.0 / 45128329728.0,
207 hf5 = 187940372067.0 / 1594534317056.0,
208 hf6 = - 1776094331.0 / 19743644256.0,
209 hf7 = 11237099.0 / 235043384.0;
213 for(
G4int i = 0; i < 3; ++i)
217 hf4 *
ak4[i] + hf5 *
ak5[i] + hf6 *
ak6[i] + hf7 *
ak7[i]);
240 const G4double bf1 = 1.0 / 11282082432.0 * (
241 157015080.0 * tau4 - 13107642775.0 * tau3 + 34969693132.0 * tau2 -
242 32272833064.0 * tau + 11282082432.0);
244 const G4double bf3 = - 100.0 / 32700410799.0 * tau * (
245 15701508.0 * tau3 - 914128567.0 * tau2 + 2074956840.0 * tau -
248 const G4double bf4 = 25.0 / 5641041216.0 * tau * (
249 94209048.0 * tau3 - 1518414297.0 * tau2 + 2460397220.0 * tau -
252 const G4double bf5 = - 2187.0 / 199316789632.0 * tau * (
253 52338360.0 * tau3 - 451824525.0 * tau2 + 687873124.0 * tau -
256 const G4double bf6 = 11.0 / 2467955532.0 * tau * (
257 106151040.0 * tau3 - 661884105.0 * tau2 +
258 946554244.0 * tau - 361440756.0);
260 const G4double bf7 = 1.0 / 29380423.0 * tau * (1.0 - tau) * (
261 8293050.0 * tau2 - 82437520.0 * tau + 44764047.0);
263 for(
G4int i = 0; i < numberOfVariables; ++i)
267 bf5 *
ak5[i] + bf6 *
ak6[i] + bf7 *
ak7[i]);
282 const G4double b81 = 6245.0 / 62208.0,
284 b83 = 8875.0 / 103032.0,
285 b84 = -125.0 / 1728.0,
286 b85 = 801.0 / 13568.0,
287 b86 = -13519.0 / 368064.0,
288 b87 = 11105.0 / 368064.0,
290 b91 = 632855.0 / 4478976.0,
292 b93 = 4146875.0 / 6491016.0,
293 b94 = 5490625.0 /14183424.0,
294 b95 = -15975.0 / 108544.0,
295 b96 = 8295925.0 / 220286304.0,
296 b97 = -1779595.0 / 62938944.0,
297 b98 = -805.0 / 4104.0;
304 for(
G4int i = 0; i < numberOfVariables; ++i)
308 b84 *
ak4[i] + b85 *
ak5[i] + b86 *
ak6[i] +
314 for(
G4int i = 0; i < numberOfVariables; ++i)
318 b94 *
ak4[i] + b95 *
ak5[i] + b96 *
ak6[i] +
319 b97 *
ak7[i] + b98 *
ak8[i]
336 bi[1][1] = -38039.0 / 7040.0,
337 bi[1][2] = 125923.0 / 10560.0,
338 bi[1][3] = -19683.0 / 1760.0,
339 bi[1][4] = 3303.0 / 880.0,
352 bi[3][1] = -12500.0 / 4081.0,
353 bi[3][2] = 205000.0 / 12243.0,
354 bi[3][3] = -90000.0 / 4081.0,
355 bi[3][4] = 36000.0 / 4081.0,
360 bi[4][1] = -3125.0 / 704.0,
361 bi[4][2] = 25625.0 / 1056.0,
362 bi[4][3] = -5625.0 / 176.0,
363 bi[4][4] = 1125.0 / 88.0,
368 bi[5][1] = 164025.0 / 74624.0,
369 bi[5][2] = -448335.0 / 37312.0,
370 bi[5][3] = 295245.0 / 18656.0,
371 bi[5][4] = -59049.0 / 9328.0,
376 bi[6][1] = -25.0 / 28.0,
377 bi[6][2] = 205.0 / 42.0,
378 bi[6][3] = -45.0 / 7.0,
379 bi[6][4] = 18.0 / 7.0,
384 bi[7][1] = -2.0 / 11.0,
385 bi[7][2] = 73.0 / 55.0,
386 bi[7][3] = -171.0 / 55.0,
387 bi[7][4] = 108.0 / 55.0,
392 bi[8][1] = 189.0 / 22.0,
393 bi[8][2] = -1593.0 / 55.0,
394 bi[8][3] = 3537.0 / 110.0,
395 bi[8][4] = -648.0 / 55.0,
400 bi[9][1] = 351.0 / 110.0,
401 bi[9][2] = -999.0 / 55.0,
402 bi[9][3] = 2943.0 / 110.0,
403 bi[9][4] = -648.0 / 55.0;
409 std::memset(b, 0.0,
sizeof(b));
412 for(
G4int j = 0; j <= 4; ++j)
414 for(
G4int iStage = 1; iStage <= 9; ++iStage)
416 b[iStage] += bi[iStage][j] * tauPower;
423 for(
G4int i = 0; i < numberOfVariables; ++i)
425 yOut[i] =
fyIn[i] + stepLen * (
427 b[4] *
ak4[i] + b[5] *
ak5[i] + b[6] *
ak6[i] +
428 b[7] *
ak7[i] + b[8] *
ak8[i] + b[9] *
ak9[i]