ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DormandPrinceRK56.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DormandPrinceRK56.cc
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 //
37 G4DormandPrinceRK56::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 //
84 G4DormandPrinceRK56::~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 //
119 void 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 //
357 G4double 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 //
401 void 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 
432 void 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 //
489 void 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 //
569 void 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 }