Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/magneticfield/src/G4DormandPrinceRK56.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /geometry/magneticfield/src/G4DormandPrinceRK56.cc (Version 11.3.0) and /geometry/magneticfield/src/G4DormandPrinceRK56.cc (Version 9.4)


  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 24 // *******************************************    
 25 //                                                
 26 // G4DormandPrinceRK56 implementation             
 27 //                                                
 28 // Created: Somnath Banerjee, Google Summer of    
 29 // Supervision: John Apostolakis, CERN            
 30 // -------------------------------------------    
 31                                                   
 32 #include "G4DormandPrinceRK56.hh"                 
 33 #include "G4LineSection.hh"                       
 34                                                   
 35 // Constructor                                    
 36 //                                                
 37 G4DormandPrinceRK56::G4DormandPrinceRK56(G4Equ    
 38                                          G4int    
 39                                          G4boo    
 40   : G4MagIntegratorStepper(EqRhs, noIntegratio    
 41 {                                                 
 42     const G4int numberOfVariables = noIntegrat    
 43                                                   
 44     // New Chunk of memory being created for u    
 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(noInte    
 65     yTemp = new G4double[numStateVars];           
 66     yIn = new G4double[numStateVars] ;            
 67                                                   
 68     fLastInitialVector = new G4double[numState    
 69     fLastFinalVector = new G4double[numStateVa    
 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(Eq    
 79     }                                             
 80 }                                                 
 81                                                   
 82 // Destructor                                     
 83 //                                                
 84 G4DormandPrinceRK56::~G4DormandPrinceRK56()       
 85 {                                                 
 86     // clear all previously allocated memory f    
 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     
117 // Giving back yOut and yErr arrays for output    
118 //                                                
119 void G4DormandPrinceRK56::Stepper(const G4doub    
120                                   const G4doub    
121                                         G4doub    
122                                         G4doub    
123                                         G4doub    
124                                  //   G4double    
125                                  //   endpoint    
126 {                                                 
127     G4int i;                                      
128                                                   
129     // The various constants defined on the ba    
130     // Old Coefficients from                      
131     // P.J.Prince and J.R.Dormand, "High order    
132     // Journal of Computational and Applied Ma    
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    
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 coefficien    
199 //                                                
200 //    T. S. Baker, J. R. Dormand, J. P. Gilmor    
201 //    "Continuous approximation with embedded     
202 //    Applied Numerical Mathematics, vol. 22,     
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 - b    
256 //    dc8 = - 1.0/2.0 - b98 ,                     
257 //    dc9 = -101.0/2294.0 ;                       
258                                                   
259     // end of declaration                         
260                                                   
261     const G4int numberOfVariables = GetNumberO    
262                                                   
263     // The number of variables to be integrate    
264     //                                            
265     yOut[7] = yTemp[7] = yIn[7] = yInput[7];      
266                                                   
267     //  Saving yInput because yInput and yOut     
268     //                                            
269     for(i=0; i<numberOfVariables; ++i)            
270     {                                             
271         yIn[i]=yInput[i];                         
272     }                                             
273     // RightHandSide(yIn, dydx) ; // 1st Stage    
274                                                   
275     for(i=0; i<numberOfVariables; ++i)            
276     {                                             
277         yTemp[i] = yIn[i] + b21*Step*dydx[i] ;    
278     }                                             
279     RightHandSide(yTemp, ak2) ;              /    
280                                                   
281     for(i=0; i<numberOfVariables; ++i)            
282     {                                             
283         yTemp[i] = yIn[i] + Step*(b31*dydx[i]     
284     }                                             
285     RightHandSide(yTemp, ak3) ;              /    
286                                                   
287     for(i=0; i<numberOfVariables; ++i)            
288     {                                             
289         yTemp[i] = yIn[i] + Step*(b41*dydx[i]     
290     }                                             
291     RightHandSide(yTemp, ak4) ;              /    
292                                                   
293     for(i=0; i<numberOfVariables; ++i)            
294     {                                             
295         yTemp[i] = yIn[i] + Step*(b51*dydx[i]     
296                                   b54*ak4[i])     
297     }                                             
298     RightHandSide(yTemp, ak5) ;              /    
299                                                   
300     for(i=0; i<numberOfVariables; ++i)            
301     {                                             
302         yTemp[i] = yIn[i] + Step*(b61*dydx[i]     
303                                   b64*ak4[i] +    
304     }                                             
305     RightHandSide(yTemp, ak6) ;              /    
306                                                   
307     for(i=0; i<numberOfVariables; ++i)            
308     {                                             
309         yTemp[i] = yIn[i] + Step*(b71*dydx[i]     
310                                   b74*ak4[i] +    
311     }                                             
312     RightHandSide(yTemp, ak7);               /    
313                                                   
314     for(i=0; i<numberOfVariables; ++i)            
315     {                                             
316         yTemp[i] = yIn[i] + Step*(b81*dydx[i]     
317                                   b84*ak4[i] +    
318                                   b87*ak7[i]);    
319     }                                             
320     RightHandSide(yTemp, ak8);               /    
321                                                   
322     for(i=0; i<numberOfVariables; ++i)            
323     {                                             
324         yOut[i] = yIn[i] + Step*(b91*dydx[i] +    
325                                   b94*ak4[i] +    
326                                   b97*ak7[i] +    
327     }                                             
328     RightHandSide(yOut, ak9);                /    
329                                                   
330                                                   
331     for(i=0; i<numberOfVariables; ++i)            
332     {                                             
333         // Estimate error as difference betwee    
334         // 6th order methods                      
335         //                                        
336         yErr[i] = Step*(  dc1*dydx[i] + dc2*ak    
337                         + dc5*ak5[i] + dc6*ak6    
338                         + dc9*ak9[i] ) ;          
339                                                   
340         // Saving 'estimated' derivative at en    
341         // nextDydx[i] = ak9[i];                  
342                                                   
343         // Store Input and Final values, for p    
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() con    
358 {                                                 
359     G4double distLine, distChord;                 
360     G4ThreeVector initialPoint, finalPoint, mi    
361                                                   
362     // Store last initial and final points        
363     // (they will be overwritten in self-Stepp    
364     //                                            
365     initialPoint = G4ThreeVector( fLastInitial    
366                                  fLastInitialV    
367     finalPoint   = G4ThreeVector( fLastFinalVe    
368                                  fLastFinalVec    
369                                                   
370     // Do half a step using StepNoErr             
371                                                   
372     fAuxStepper->Stepper( fLastInitialVector,     
373                          fMidVector,   fMidErr    
374                                                   
375     midPoint = G4ThreeVector( fMidVector[0], f    
376                                                   
377     // Use stored values of Initial and Endpoi    
378     // distance of Chord                          
379     //                                            
380     if (initialPoint != finalPoint)               
381     {                                             
382         distLine  = G4LineSection::Distline( m    
383         distChord = distLine;                     
384     }                                             
385     else                                          
386     {                                             
387         distChord = (midPoint-initialPoint).ma    
388     }                                             
389     return distChord;                             
390 }                                                 
391                                                   
392 // The following interpolation scheme has been    
393 // Table 5. The RK6(5)9FM process and associat    
394 //                                                
395 // J. R. Dormand, M. A. Lockyer, N. E. McGorri    
396 // "Global error estimation with runge-kutta t    
397 // Computers & Mathematics with Applications,     
398 //                                                
399 // Fifth order interpolant with one extra func    
400 //                                                
401 void G4DormandPrinceRK56::SetupInterpolate_low    
402                                                   
403                                                   
404 {                                                 
405     const G4int numberOfVariables= this->GetNu    
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    
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]     
426                                 b_104*ak4[i] +    
427                                 b_107*ak7[i] +    
428     }                                             
429     RightHandSide(yTemp, ak10_low);          /    
430 }                                                 
431                                                   
432 void G4DormandPrinceRK56::Interpolate_low( con    
433                                            con    
434                                            con    
435                                                   
436                                                   
437 {                                                 
438   G4double bf1, bf4, bf5, bf6, bf7, bf8, bf9,     
439                                                   
440   G4double tau0 = tau;                            
441   const G4int numberOfVariables= this->GetNumb    
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    
454       / 28800.0 ;                                 
455   bf4 = -16.0*tau*(45312.0*tau_3 - 125933.0*ta    
456       / 70785.0 ;                                 
457   bf5 = -2187.0*tau*(19440.0*tau_3 - 45743.0*t    
458       / 1645600.0 ;                               
459   bf6 = tau*(12864.0*tau_3 - 30653.0*tau_2 + 2    
460       / 705.0 ;                                   
461   bf7 = -5764801.0*tau*(16464.0*tau_3 - 32797.    
462       / 7239323520.0 ;                            
463   bf8 =  37.0*tau*(336.0*tau_3 - 661.0*tau_2 +    
464       / 1440.0 ;                                  
465   bf9 = tau*(tau-1.0)*(16.0*tau_2 - 15.0*tau +    
466       / 4.0 ;                                     
467   bf10 = 8.0*tau*(tau - 1.0)*(tau - 1.0)*(2.0*    
468                                                   
469   for( G4int i=0; i<numberOfVariables; ++i)       
470   {                                               
471     yOut[i] = yIn[i] + Step*tau*( bf1*dydx[i]     
472                                   bf6*ak6[i] +    
473                                   bf9*ak9[i] +    
474   }                                               
475 }                                                 
476                                                   
477 // The following scheme and set of coefficient    
478 // Table 2. Sixth order dense formula based on    
479 // RK6(5)9FM with extra stages C1O= 1/2, C11 =    
480 //                                                
481 // T. S. Baker, J. R. Dormand, J. P. Gilmore,     
482 // "Continuous approximation with embedded Run    
483 // Applied Numerical Mathematics, vol. 22, no.    
484 //                                                
485 // --- Sixth order interpolant with 3 addition    
486 //                                                
487 // Function for calculating the additional sta    
488 //                                                
489 void G4DormandPrinceRK56::SetupInterpolate_hig    
490                                                   
491                                                   
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.    
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/83630285    
523              b128 =  2725085557.0/26167173120.    
524              b129 =  235429367.0/16354483200.0    
525              b1210 = -90924917.0/1040739840.0     
526              b1211 = -271149.0/21414400.0 ;       
527                                                   
528     const G4int numberOfVariables = GetNumberO    
529                                                   
530     // Saving yInput because yInput and yOut c    
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]    
543                                   b104*ak4[i]     
544                                   b107*ak7[i]     
545     }                                             
546     RightHandSide(yTemp, ak10);                   
547                                                   
548     for(G4int i=0; i<numberOfVariables; ++i)      
549     {                                             
550         yTemp[i] = yIn[i] + Step*(b111*dydx[i]    
551                                   b114*ak4[i]     
552                                   b117*ak7[i]     
553                                   b1110*ak10[i    
554     }                                             
555     RightHandSide(yTemp, ak11);                   
556                                                   
557     for(G4int i=0; i<numberOfVariables; ++i)      
558     {                                             
559         yTemp[i] = yIn[i] + Step*(b121*dydx[i]    
560                                   b124*ak4[i]     
561                                   b127*ak7[i]     
562                                   b1210*ak10[i    
563     }                                             
564     RightHandSide(yTemp, ak12);                   
565 }                                                 
566                                                   
567 // Function to interpolate to tau(passed in) f    
568 //                                                
569 void G4DormandPrinceRK56::Interpolate_high( co    
570                                             co    
571                                             co    
572                                                   
573                                                   
574 {                                                 
575     // Define the coefficients for the polynom    
576     //                                            
577     G4double bi[13][6], b[13];                    
578     G4int numberOfVariables = GetNumberOfVaria    
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 (coefficent    
697     //                                            
698     for(auto i=1; i<=12; ++i) // i is NOT the     
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 fr    
710     // the polynomial coefficients and the res    
711     //                                            
712     for(G4int i=0; i<numberOfVariables; ++i) /    
713     {                                             
714       yOut[i] = yIn[i] + Step*tau0*(b[1]*dydx[    
715                                     b[4]*ak4[i    
716                                     b[7]*ak7[i    
717                                  b[10]*ak10[i]    
718     }                                             
719 }                                                 
720