Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4DormandPrinceRK56.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// G4DormandPrinceRK56 implementation
27//
28// Created: Somnath Banerjee, Google Summer of Code 2015, 26 June 2015
29// Supervision: John Apostolakis, CERN
30// --------------------------------------------------------------------
31
32#include "G4DormandPrinceRK56.hh"
33#include "G4LineSection.hh"
34
35// Constructor
36//
37G4DormandPrinceRK56::G4DormandPrinceRK56(G4EquationOfMotion* EqRhs,
38 G4int noIntegrationVariables,
39 G4bool primary)
40 : G4MagIntegratorStepper(EqRhs, noIntegrationVariables)
41{
42 const G4int numberOfVariables = noIntegrationVariables;
43
44 // New Chunk of memory being created for use by the stepper
45
46 // aki - for storing intermediate RHS
47 //
48 ak2 = new G4double[numberOfVariables];
49 ak3 = new G4double[numberOfVariables];
50 ak4 = new G4double[numberOfVariables];
51 ak5 = new G4double[numberOfVariables];
52 ak6 = new G4double[numberOfVariables];
53 ak7 = new G4double[numberOfVariables];
54 ak8 = new G4double[numberOfVariables];
55 ak9 = new G4double[numberOfVariables];
56
57 // Memory for Additional stages
58 //
59 ak10 = new G4double[numberOfVariables];
60 ak11 = new G4double[numberOfVariables];
61 ak12 = new G4double[numberOfVariables];
62 ak10_low = new G4double[numberOfVariables];
63
64 const G4int numStateVars = std::max(noIntegrationVariables, 8);
65 yTemp = new G4double[numStateVars];
66 yIn = new G4double[numStateVars] ;
67
68 fLastInitialVector = new G4double[numStateVars] ;
69 fLastFinalVector = new G4double[numStateVars] ;
70
71 fLastDyDx = new G4double[numStateVars];
72
73 fMidVector = new G4double[numStateVars];
74 fMidError = new G4double[numStateVars];
75
76 if( primary )
77 {
78 fAuxStepper = new G4DormandPrinceRK56(EqRhs, numberOfVariables, !primary);
79 }
80}
81
82// Destructor
83//
84G4DormandPrinceRK56::~G4DormandPrinceRK56()
85{
86 // clear all previously allocated memory for stepper and DistChord
87
88 delete [] ak2;
89 delete [] ak3;
90 delete [] ak4;
91 delete [] ak5;
92 delete [] ak6;
93 delete [] ak7;
94 delete [] ak8;
95 delete [] ak9;
96
97 delete [] ak10;
98 delete [] ak10_low;
99 delete [] ak11;
100 delete [] ak12;
101
102 delete [] yTemp;
103 delete [] yIn;
104
105 delete [] fLastInitialVector;
106 delete [] fLastFinalVector;
107 delete [] fLastDyDx;
108 delete [] fMidVector;
109 delete [] fMidError;
110
111 delete fAuxStepper;
112}
113
114// Stepper
115//
116// Passing in the value of yInput[],the first time dydx[] and Step length
117// Giving back yOut and yErr arrays for output and error respectively
118//
119void G4DormandPrinceRK56::Stepper(const G4double yInput[],
120 const G4double dydx[],
121 G4double Step,
122 G4double yOut[],
123 G4double yErr[] )
124 // G4double nextDydx[] ) -- Output:
125 // endpoint DyDx ( for future FSAL version )
126{
127 G4int i;
128
129 // The various constants defined on the basis of butcher tableu
130 // Old Coefficients from
131 // P.J.Prince and J.R.Dormand, "High order embedded Runge-Kutta formulae"
132 // Journal of Computational and Applied Math., vol.7, no.1, pp.67-75, 1980.
133 //
134 const G4double b21 = 1.0/10.0 ,
135 b31 = -2.0/81.0 ,
136 b32 = 20.0/81.0 ,
137
138 b41 = 615.0/1372.0 ,
139 b42 = -270.0/343.0 ,
140 b43 = 1053.0/1372.0 ,
141
142 b51 = 3243.0/5500.0 ,
143 b52 = -54.0/55.0 ,
144 b53 = 50949.0/71500.0 ,
145 b54 = 4998.0/17875.0 ,
146
147 b61 = -26492.0/37125.0 ,
148 b62 = 72.0/55.0 ,
149 b63 = 2808.0/23375.0 ,
150 b64 = -24206.0/37125.0 ,
151 b65 = 338.0/459.0 ,
152
153 b71 = 5561.0/2376.0 ,
154 b72 = -35.0/11.0 ,
155 b73 = -24117.0/31603.0 ,
156 b74 = 899983.0/200772.0 ,
157 b75 = -5225.0/1836.0 ,
158 b76 = 3925.0/4056.0 ,
159
160 b81 = 465467.0/266112.0 ,
161 b82 = -2945.0/1232.0 ,
162 b83 = -5610201.0/14158144.0 ,
163 b84 = 10513573.0/3212352.0 ,
164 b85 = -424325.0/205632.0 ,
165 b86 = 376225.0/454272.0 ,
166 b87 = 0.0 ,
167
168 c1 = 61.0/864.0 ,
169 c2 = 0.0 ,
170 c3 = 98415.0/321776.0 ,
171 c4 = 16807.0/146016.0 ,
172 c5 = 1375.0/7344.0 ,
173 c6 = 1375.0/5408.0 ,
174 c7 = -37.0/1120.0 ,
175 c8 = 1.0/10.0 ,
176
177 b91 = 61.0/864.0 ,
178 b92 = 0.0 ,
179 b93 = 98415.0/321776.0 ,
180 b94 = 16807.0/146016.0 ,
181 b95 = 1375.0/7344.0 ,
182 b96 = 1375.0/5408.0 ,
183 b97 = -37.0/1120.0 ,
184 b98 = 1.0/10.0 ,
185
186 dc1 = c1 - 821.0/10800.0 ,
187 dc2 = c2 - 0.0 ,
188 dc3 = c3 - 19683.0/71825,
189 dc4 = c4 - 175273.0/912600.0 ,
190 dc5 = c5 - 395.0/3672.0 ,
191 dc6 = c6 - 785.0/2704.0 ,
192 dc7 = c7 - 3.0/50.0 ,
193 dc8 = c8 - 0.0 ,
194 dc9 = 0.0;
195
196
197// New Coefficients obtained from
198// Table 3 RK6(5)9FM with corrected coefficients
199//
200// T. S. Baker, J. R. Dormand, J. P. Gilmore, and P. J. Prince,
201// "Continuous approximation with embedded Runge-Kutta methods"
202// Applied Numerical Mathematics, vol. 22, no. 1, pp. 51-62, 1996.
203//
204// b21 = 1.0/9.0 ,
205//
206// b31 = 1.0/24.0 ,
207// b32 = 1.0/8.0 ,
208//
209// b41 = 1.0/16.0 ,
210// b42 = 0.0 ,
211// b43 = 3.0/16.0 ,
212//
213// b51 = 280.0/729.0 ,
214// b52 = 0.0 ,
215// b53 = -325.0/243.0 ,
216// b54 = 1100.0/729.0 ,
217//
218// b61 = 6127.0/14680.0 ,
219// b62 = 0.0 ,
220// b63 = -1077.0/734.0 ,
221// b64 = 6494.0/4037.0 ,
222// b65 = -9477.0/161480.0 ,
223//
224// b71 = -13426273320.0/14809773769.0 ,
225// b72 = 0.0 ,
226// b73 = 4192558704.0/2115681967.0 ,
227// b74 = 14334750144.0/14809773769.0 ,
228// b75 = 117092732328.0/14809773769.0 ,
229// b76 = -361966176.0/40353607.0 ,
230//
231// b81 = -2340689.0/1901060.0 ,
232// b82 = 0.0 ,
233// b83 = 31647.0/13579.0 ,
234// b84 = 253549596.0/149518369.0 ,
235// b85 = 10559024082.0/977620105.0 ,
236// b86 = -152952.0/12173.0 ,
237// b87 = -5764801.0/186010396.0 ,
238//
239// b91 = 203.0/2880.0 ,
240// b92 = 0.0 ,
241// b93 = 0.0 ,
242// b94 = 30208.0/70785.0 ,
243// b95 = 177147.0/164560.0 ,
244// b96 = -536.0/705.0 ,
245// b97 = 1977326743.0/3619661760.0 ,
246// b98 = -259.0/720.0 ,
247//
248//
249// dc1 = 36567.0/458800.0 - b91,
250// dc2 = 0.0 - b92,
251// dc3 = 0.0 - b93,
252// dc4 = 9925984.0/27063465.0 - b94,
253// dc5 = 85382667.0/117968950.0 - b95,
254// dc6 = - 310378.0/808635.0 - b96 ,
255// dc7 = 262119736669.0/345979336560.0 - b97,
256// dc8 = - 1.0/2.0 - b98 ,
257// dc9 = -101.0/2294.0 ;
258
259 // end of declaration
260
261 const G4int numberOfVariables = GetNumberOfVariables();
262
263 // The number of variables to be integrated over
264 //
265 yOut[7] = yTemp[7] = yIn[7] = yInput[7];
266
267 // Saving yInput because yInput and yOut can be aliases for same array
268 //
269 for(i=0; i<numberOfVariables; ++i)
270 {
271 yIn[i]=yInput[i];
272 }
273 // RightHandSide(yIn, dydx) ; // 1st Stage - Not doing, getting passed
274
275 for(i=0; i<numberOfVariables; ++i)
276 {
277 yTemp[i] = yIn[i] + b21*Step*dydx[i] ;
278 }
279 RightHandSide(yTemp, ak2) ; // 2nd Stage
280
281 for(i=0; i<numberOfVariables; ++i)
282 {
283 yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b32*ak2[i]) ;
284 }
285 RightHandSide(yTemp, ak3) ; // 3rd Stage
286
287 for(i=0; i<numberOfVariables; ++i)
288 {
289 yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]) ;
290 }
291 RightHandSide(yTemp, ak4) ; // 4th Stage
292
293 for(i=0; i<numberOfVariables; ++i)
294 {
295 yTemp[i] = yIn[i] + Step*(b51*dydx[i] + b52*ak2[i] + b53*ak3[i] +
296 b54*ak4[i]) ;
297 }
298 RightHandSide(yTemp, ak5) ; // 5th Stage
299
300 for(i=0; i<numberOfVariables; ++i)
301 {
302 yTemp[i] = yIn[i] + Step*(b61*dydx[i] + b62*ak2[i] + b63*ak3[i] +
303 b64*ak4[i] + b65*ak5[i]) ;
304 }
305 RightHandSide(yTemp, ak6) ; // 6th Stage
306
307 for(i=0; i<numberOfVariables; ++i)
308 {
309 yTemp[i] = yIn[i] + Step*(b71*dydx[i] + b72*ak2[i] + b73*ak3[i] +
310 b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
311 }
312 RightHandSide(yTemp, ak7); // 7th Stage
313
314 for(i=0; i<numberOfVariables; ++i)
315 {
316 yTemp[i] = yIn[i] + Step*(b81*dydx[i] + b82*ak2[i] + b83*ak3[i] +
317 b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
318 b87*ak7[i]);
319 }
320 RightHandSide(yTemp, ak8); // 8th Stage
321
322 for(i=0; i<numberOfVariables; ++i)
323 {
324 yOut[i] = yIn[i] + Step*(b91*dydx[i] + b92*ak2[i] + b93*ak3[i] +
325 b94*ak4[i] + b95*ak5[i] + b96*ak6[i] +
326 b97*ak7[i] + b98*ak8[i] );
327 }
328 RightHandSide(yOut, ak9); // 9th Stage
329
330
331 for(i=0; i<numberOfVariables; ++i)
332 {
333 // Estimate error as difference between 5th and
334 // 6th order methods
335 //
336 yErr[i] = Step*( dc1*dydx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i]
337 + dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] + dc8*ak8[i]
338 + dc9*ak9[i] ) ;
339
340 // Saving 'estimated' derivative at end-point
341 // nextDydx[i] = ak9[i];
342
343 // Store Input and Final values, for possible use in calculating chord
344 //
345 fLastInitialVector[i] = yIn[i] ;
346 fLastFinalVector[i] = yOut[i];
347 fLastDyDx[i] = dydx[i];
348 }
349
350 fLastStepLength = Step;
351
352 return ;
353}
354
355// DistChord
356//
357G4double G4DormandPrinceRK56::DistChord() const
358{
359 G4double distLine, distChord;
360 G4ThreeVector initialPoint, finalPoint, midPoint;
361
362 // Store last initial and final points
363 // (they will be overwritten in self-Stepper call!)
364 //
365 initialPoint = G4ThreeVector( fLastInitialVector[0],
366 fLastInitialVector[1], fLastInitialVector[2]);
367 finalPoint = G4ThreeVector( fLastFinalVector[0],
368 fLastFinalVector[1], fLastFinalVector[2]);
369
370 // Do half a step using StepNoErr
371
372 fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
373 fMidVector, fMidError );
374
375 midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
376
377 // Use stored values of Initial and Endpoint + new Midpoint to evaluate
378 // distance of Chord
379 //
380 if (initialPoint != finalPoint)
381 {
382 distLine = G4LineSection::Distline( midPoint,initialPoint,finalPoint );
383 distChord = distLine;
384 }
385 else
386 {
387 distChord = (midPoint-initialPoint).mag();
388 }
389 return distChord;
390}
391
392// The following interpolation scheme has been obtained from
393// Table 5. The RK6(5)9FM process and associated dense formula
394//
395// J. R. Dormand, M. A. Lockyer, N. E. McGorrigan, and P. J. Prince,
396// "Global error estimation with runge-kutta triples"
397// Computers & Mathematics with Applications, vol.18, no.9, pp.835-846, 1989.
398//
399// Fifth order interpolant with one extra function evaluation per step
400//
401void G4DormandPrinceRK56::SetupInterpolate_low( const G4double yInput[],
402 const G4double dydx[],
403 const G4double Step )
404{
405 const G4int numberOfVariables= this->GetNumberOfVariables();
406
407 G4double b_101 = 33797.0/460800.0 ,
408 b_102 = 0. ,
409 b_103 = 0. ,
410 b_104 = 27757.0/70785.0 ,
411 b_105 = 7923501.0/26329600.0 ,
412 b_106 = -927.0/3760.0 ,
413 b_107 = -3314760575.0/23165835264.0 ,
414 b_108 = 2479.0/23040.0 ,
415 b_109 = 1.0/64.0 ;
416
417 for(G4int i=0; i<numberOfVariables; ++i)
418 {
419 yIn[i]=yInput[i];
420 }
421
422
423 for(G4int i=0; i<numberOfVariables; ++i)
424 {
425 yTemp[i] = yIn[i] + Step*(b_101*dydx[i] + b_102*ak2[i] + b_103*ak3[i] +
426 b_104*ak4[i] + b_105*ak5[i] + b_106*ak6[i] +
427 b_107*ak7[i] + b_108*ak8[i] + b_109*ak9[i]);
428 }
429 RightHandSide(yTemp, ak10_low); // 10th Stage
430}
431
432void G4DormandPrinceRK56::Interpolate_low( const G4double yInput[],
433 const G4double dydx[],
434 const G4double Step,
435 G4double yOut[],
436 G4double tau )
437{
438 G4double bf1, bf4, bf5, bf6, bf7, bf8, bf9, bf10;
439
440 G4double tau0 = tau;
441 const G4int numberOfVariables= this->GetNumberOfVariables();
442
443 for(G4int i=0; i<numberOfVariables; ++i)
444 {
445 yIn[i]=yInput[i];
446 }
447
448 G4double tau_2 = tau0*tau0 ,
449 tau_3 = tau0*tau_2,
450 tau_4 = tau_2*tau_2;
451
452 // bf2 = bf3 = 0.0
453 bf1 = (66480.0*tau_4-206243.0*tau_3+237786.0*tau_2-124793.0*tau+28800.0)
454 / 28800.0 ;
455 bf4 = -16.0*tau*(45312.0*tau_3 - 125933.0*tau_2 + 119706.0*tau -40973.0)
456 / 70785.0 ;
457 bf5 = -2187.0*tau*(19440.0*tau_3 - 45743.0*tau_2 + 34786.0*tau - 9293.0)
458 / 1645600.0 ;
459 bf6 = tau*(12864.0*tau_3 - 30653.0*tau_2 + 23786.0*tau - 6533.0)
460 / 705.0 ;
461 bf7 = -5764801.0*tau*(16464.0*tau_3 - 32797.0*tau_2 + 17574.0*tau - 1927.0)
462 / 7239323520.0 ;
463 bf8 = 37.0*tau*(336.0*tau_3 - 661.0*tau_2 + 342.0*tau -31.0)
464 / 1440.0 ;
465 bf9 = tau*(tau-1.0)*(16.0*tau_2 - 15.0*tau +3.0)
466 / 4.0 ;
467 bf10 = 8.0*tau*(tau - 1.0)*(tau - 1.0)*(2.0*tau - 1.0) ;
468
469 for( G4int i=0; i<numberOfVariables; ++i)
470 {
471 yOut[i] = yIn[i] + Step*tau*( bf1*dydx[i] + bf4*ak4[i] + bf5*ak5[i] +
472 bf6*ak6[i] + bf7*ak7[i] + bf8*ak8[i] +
473 bf9*ak9[i] + bf10*ak10_low[i] ) ;
474 }
475}
476
477// The following scheme and set of coefficients have been obtained from
478// Table 2. Sixth order dense formula based on linear optimisation for
479// RK6(5)9FM with extra stages C1O= 1/2, C11 =1/6, c12= 5/12
480//
481// T. S. Baker, J. R. Dormand, J. P. Gilmore, and P. J. Prince,
482// "Continuous approximation with embedded Runge-Kutta methods"
483// Applied Numerical Mathematics, vol. 22, no. 1, pp. 51-62, 1996.
484//
485// --- Sixth order interpolant with 3 additional stages per step ---
486//
487// Function for calculating the additional stages
488//
489void G4DormandPrinceRK56::SetupInterpolate_high( const G4double yInput[],
490 const G4double dydx[],
491 const G4double Step )
492{
493 // Coefficients for the additional stages
494 //
495 G4double b101 = 33797.0/460800.0 ,
496 b102 = 0.0 ,
497 b103 = 0.0 ,
498 b104 = 27757.0/70785.0 ,
499 b105 = 7923501.0/26329600.0 ,
500 b106 = -927.0/3760.0 ,
501 b107 = -3314760575.0/23165835264.0 ,
502 b108 = 2479.0/23040.0 ,
503 b109 = 1.0/64.0 ,
504
505 b111 = 5843.0/76800.0 ,
506 b112 = 0.0 ,
507 b113 = 0.0 ,
508 b114 = 464.0/2673.0 ,
509 b115 = 353997.0/1196800.0 ,
510 b116 = -15068.0/57105.0 ,
511 b117 = -282475249.0/3644974080.0 ,
512 b118 = 8678831.0/156245760.0 ,
513 b119 = 116113.0/11718432.0 ,
514 b1110 = -25.0/243.0 ,
515
516 b121 = 15088049.0/199065600.0 ,
517 b122 = 0.0 ,
518 b123 = 0.0 ,
519 b124 = 2.0/5.0 ,
520 b125 = 92222037.0/268083200.0 ,
521 b126 = -433420501.0/1528586640.0 ,
522 b127 = -11549242677007.0/83630285291520.0 ,
523 b128 = 2725085557.0/26167173120.0 ,
524 b129 = 235429367.0/16354483200.0 ,
525 b1210 = -90924917.0/1040739840.0 ,
526 b1211 = -271149.0/21414400.0 ;
527
528 const G4int numberOfVariables = GetNumberOfVariables();
529
530 // Saving yInput because yInput and yOut can be aliases for same array
531 //
532 for(G4int i=0; i<numberOfVariables; ++i)
533 {
534 yIn[i]=yInput[i];
535 }
536 yTemp[7] = yIn[7];
537
538 // Evaluate the extra stages
539 //
540 for(G4int i=0; i<numberOfVariables; ++i)
541 {
542 yTemp[i] = yIn[i] + Step*(b101*dydx[i] + b102*ak2[i] + b103*ak3[i] +
543 b104*ak4[i] + b105*ak5[i] + b106*ak6[i] +
544 b107*ak7[i] + b108*ak8[i] + b109*ak9[i]);
545 }
546 RightHandSide(yTemp, ak10); // 10th Stage
547
548 for(G4int i=0; i<numberOfVariables; ++i)
549 {
550 yTemp[i] = yIn[i] + Step*(b111*dydx[i] + b112*ak2[i] + b113*ak3[i] +
551 b114*ak4[i] + b115*ak5[i] + b116*ak6[i] +
552 b117*ak7[i] + b118*ak8[i] + b119*ak9[i] +
553 b1110*ak10[i]);
554 }
555 RightHandSide(yTemp, ak11); // 11th Stage
556
557 for(G4int i=0; i<numberOfVariables; ++i)
558 {
559 yTemp[i] = yIn[i] + Step*(b121*dydx[i] + b122*ak2[i] + b123*ak3[i] +
560 b124*ak4[i] + b125*ak5[i] + b126*ak6[i] +
561 b127*ak7[i] + b128*ak8[i] + b129*ak9[i] +
562 b1210*ak10[i] + b1211*ak11[i]);
563 }
564 RightHandSide(yTemp, ak12); // 12th Stage
565}
566
567// Function to interpolate to tau(passed in) fraction of the step
568//
569void G4DormandPrinceRK56::Interpolate_high( const G4double yInput[],
570 const G4double dydx[],
571 const G4double Step,
572 G4double yOut[],
573 G4double tau )
574{
575 // Define the coefficients for the polynomials
576 //
577 G4double bi[13][6], b[13];
578 G4int numberOfVariables = GetNumberOfVariables();
579
580
581 // COEFFICIENTS OF bi[ 1]
582 bi[1][0] = 1.0 ,
583 bi[1][1] = -18487.0/2880.0 ,
584 bi[1][2] = 139189.0/7200.0 ,
585 bi[1][3] = -53923.0/1800.0 ,
586 bi[1][4] = 13811.0/600,
587 bi[1][5] = -2071.0/300,
588 // --------------------------------------------------------
589 //
590 // COEFFICIENTS OF bi[2]
591 bi[2][0] = 0.0 ,
592 bi[2][1] = 0.0 ,
593 bi[2][2] = 0.0 ,
594 bi[2][3] = 0.0 ,
595 bi[2][4] = 0.0 ,
596 bi[2][5] = 0.0 ,
597 // --------------------------------------------------------
598 //
599 // COEFFICIENTS OF bi[3]
600 bi[3][0] = 0.0 ,
601 bi[3][1] = 0.0 ,
602 bi[3][2] = 0.0 ,
603 bi[3][3] = 0.0 ,
604 bi[3][4] = 0.0 ,
605 bi[3][5] = 0.0 ,
606 // --------------------------------------------------------
607 //
608 // COEFFICIENTS OF bi[4]
609 bi[4][0] = 0.0 ,
610 bi[4][1] = -30208.0/14157.0 ,
611 bi[4][2] = 1147904.0/70785.0 ,
612 bi[4][3] = -241664.0/5445.0 ,
613 bi[4][4] = 241664.0/4719.0 ,
614 bi[4][5] = -483328.0/23595.0 ,
615 // --------------------------------------------------------
616 //
617 // COEFFICIENTS OF bi[5]
618 bi[5][0] = 0.0 ,
619 bi[5][1] = -177147.0/32912.0 ,
620 bi[5][2] = 3365793.0/82280.0 ,
621 bi[5][3] = -2302911.0/20570.0 ,
622 bi[5][4] = 531441.0/4114.0 ,
623 bi[5][5] = -531441.0/10285.0 ,
624 // --------------------------------------------------------
625 //
626 // COEFFICIENTS OF bi[6]
627 bi[6][0] = 0.0 ,
628 bi[6][1] = 536.0/141.0 ,
629 bi[6][2] = -20368.0/705.0 ,
630 bi[6][3] = 55744.0/705.0 ,
631 bi[6][4] = -4288.0/47.0 ,
632 bi[6][5] = 8576.0/235,
633 // --------------------------------------------------------
634 //
635 // COEFFICIENTS OF bi[7]
636 bi[7][0] = 0.0 ,
637 bi[7][1] = -1977326743.0/723932352.0 ,
638 bi[7][2] = 37569208117.0/1809830880.0 ,
639 bi[7][3] = -1977326743.0/34804440.0 ,
640 bi[7][4] = 1977326743.0/30163848.0 ,
641 bi[7][5] = -1977326743.0/75409620.0 ,
642 // --------------------------------------------------------
643 //
644 // COEFFICIENTS OF bi[8]
645 bi[8][0] = 0.0 ,
646 bi[8][1] = 259.0/144.0 ,
647 bi[8][2] = -4921.0/360.0 ,
648 bi[8][3] = 3367.0/90.0 ,
649 bi[8][4] = -259.0/6.0 ,
650 bi[8][5] = 259.0/15.0 ,
651 // --------------------------------------------------------
652 //
653 // COEFFICIENTS OF bi[9]
654 bi[9][0] = 0.0 ,
655 bi[9][1] = 62.0/105.0 ,
656 bi[9][2] = -2381.0/525.0 ,
657 bi[9][3] = 949.0/75.0 ,
658 bi[9][4] = -2636.0/175.0 ,
659 bi[9][5] = 1112.0/175.0 ,
660 // --------------------------------------------------------
661 //
662 // COEFFICIENTS OF bi[10]
663 bi[10][0] = 0.0 ,
664 bi[10][1] = 43.0/3.0 ,
665 bi[10][2] = -1534.0/15.0 ,
666 bi[10][3] = 3767.0/15.0 ,
667 bi[10][4] = -1264.0/5.0 ,
668 bi[10][5] = 448.0/5.0 ,
669 // --------------------------------------------------------
670 //
671 // COEFFICIENTS OF bi[11]
672 bi[11][0] = 0.0 ,
673 bi[11][1] = 63.0/5.0 ,
674 bi[11][2] = -1494.0/25.0 ,
675 bi[11][3] = 2907.0/25.0 ,
676 bi[11][4] = -2592.0/25.0 ,
677 bi[11][5] = 864.0/25.0 ,
678 // --------------------------------------------------------
679 //
680 // COEFFICIENTS OF bi[12]
681 bi[12][0] = 0.0 ,
682 bi[12][1] = -576.0/35.0 ,
683 bi[12][2] = 19584.0/175.0 ,
684 bi[12][3] = -6336.0/25.0 ,
685 bi[12][4] = 41472.0/175.0 ,
686 bi[12][5] = -13824.0/175.0 ;
687 // --------------------------------------------------------
688
689 for(G4int i = 0; i< numberOfVariables; ++i)
690 {
691 yIn[i] = yInput[i];
692 }
693
694 G4double tau0 = tau;
695
696 // Calculating the polynomials (coefficents for the respective stages) :
697 //
698 for(auto i=1; i<=12; ++i) // i is NOT the coordinate no., it's stage no.
699 {
700 b[i] = 0;
701 tau = 1.0;
702 for(auto j=0; j<=5; ++j)
703 {
704 b[i] += bi[i][j]*tau;
705 tau*=tau0;
706 }
707 }
708
709 // Calculating the interpolation at the fraction tau of the step using
710 // the polynomial coefficients and the respective stages
711 //
712 for(G4int i=0; i<numberOfVariables; ++i) // Here i IS the coordinate no.
713 {
714 yOut[i] = yIn[i] + Step*tau0*(b[1]*dydx[i] + b[2]*ak2[i] + b[3]*ak3[i] +
715 b[4]*ak4[i] + b[5]*ak5[i] + b[6]*ak6[i] +
716 b[7]*ak7[i] + b[8]*ak8[i] + b[9]*ak9[i] +
717 b[10]*ak10[i] + b[11]*ak11[i] + b[12]*ak12[i]);
718 }
719}
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
G4double
double G4double
Definition G4Types.hh:83
G4bool
bool G4bool
Definition G4Types.hh:86
G4int
int G4int
Definition G4Types.hh:85
G4DormandPrinceRK56::SetupInterpolate_low
void SetupInterpolate_low(const G4double yInput[], const G4double dydx[], const G4double Step)
Definition G4DormandPrinceRK56.cc:401
G4DormandPrinceRK56::G4DormandPrinceRK56
G4DormandPrinceRK56(G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
Definition G4DormandPrinceRK56.cc:37
G4DormandPrinceRK56::Interpolate_low
void Interpolate_low(const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
Definition G4DormandPrinceRK56.cc:432
G4DormandPrinceRK56::Stepper
void Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[]) override
Definition G4DormandPrinceRK56.cc:119
G4DormandPrinceRK56::Interpolate_high
void Interpolate_high(const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
Definition G4DormandPrinceRK56.cc:569
G4DormandPrinceRK56::DistChord
G4double DistChord() const override
Definition G4DormandPrinceRK56.cc:357
G4DormandPrinceRK56::SetupInterpolate_high
void SetupInterpolate_high(const G4double yInput[], const G4double dydx[], const G4double Step)
Definition G4DormandPrinceRK56.cc:489
G4LineSection::Distline
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
Definition G4LineSection.hh:72
G4MagIntegratorStepper::GetNumberOfVariables
G4int GetNumberOfVariables() const
G4MagIntegratorStepper::RightHandSide
void RightHandSide(const G4double y[], G4double dydx[]) const