176{
177
178
179 const G4double b21 = 0.2, b31 = 3.0 / 40.0, b32 = 9.0 / 40.0, b41 = 0.3,
180 b42 = -0.9, b43 = 1.2,
181
182 b51 = -11.0 / 54.0, b52 = 2.5, b53 = -70.0 / 27.0,
183 b54 = 35.0 / 27.0,
184
185 b61 = 1631.0 / 55296.0, b62 = 175.0 / 512.0,
186 b63 = 575.0 / 13824.0, b64 = 44275.0 / 110592.0,
187 b65 = 253.0 / 4096.0,
188
189 c1 = 37.0 / 378.0, c3 = 250.0 / 621.0, c4 = 125.0 / 594.0,
190 c6 = 512.0 / 1771.0, dc5 = -277.0 / 14336.0;
191
192 const G4double dc1 = c1 - 2825.0 / 27648.0, dc3 = c3 - 18575.0 / 48384.0,
193 dc4 = c4 - 13525.0 / 55296.0, dc6 = c6 - 0.25;
194
195
196
197
198
199
200
201 for(
unsigned int i = 0; i <
N; ++i)
202 {
203 yIn[i] = yInput[i];
204 }
205
206
207 for(
unsigned int i = 0; i <
N; ++i)
208 {
209 yTemp[i] = yIn[i] + b21 * Step * dydx[i];
210 }
212
213 for(
unsigned int i = 0; i <
N; ++i)
214 {
215 yTemp[i] = yIn[i] + Step * (b31 * dydx[i] + b32 * ak2[i]);
216 }
218
219 for(
unsigned int i = 0; i <
N; ++i)
220 {
221 yTemp[i] = yIn[i] + Step * (b41 * dydx[i] + b42 * ak2[i] + b43 * ak3[i]);
222 }
224
225 for(
unsigned int i = 0; i <
N; ++i)
226 {
227 yTemp[i] = yIn[i] + Step * (b51 * dydx[i] + b52 * ak2[i] + b53 * ak3[i] +
228 b54 * ak4[i]);
229 }
231
232 for(
unsigned int i = 0; i <
N; ++i)
233 {
234 yTemp[i] = yIn[i] + Step * (b61 * dydx[i] + b62 * ak2[i] + b63 * ak3[i] +
235 b64 * ak4[i] + b65 * ak5[i]);
236 }
238
239 for(
unsigned int i = 0; i <
N; ++i)
240 {
241
242
243 yOut[i] = yIn[i] +
244 Step * (c1 * dydx[i] + c3 * ak3[i] + c4 * ak4[i] + c6 * ak6[i]);
245 }
246 for(
unsigned int i = 0; i <
N; ++i)
247 {
248
249
250
251 yErr[i] = Step * (dc1 * dydx[i] + dc3 * ak3[i] + dc4 * ak4[i] +
252 dc5 * ak5[i] + dc6 * ak6[i]);
253 }
254 for(
unsigned int i = 0; i <
N; ++i)
255 {
256
257 fLastInitialVector[i] = yIn[i];
258 fLastFinalVector[i] = yOut[i];
259 fLastDyDx[i] = dydx[i];
260 }
261
262
263 fLastStepLength = Step;
264
265 return;
266}
void RightHandSideInl(const G4double y[], G4double dydx[])