Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/magneticfield/src/G4BogackiShampine45.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 ]

  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 // G4BogackiShampine45 implementation
 27 //
 28 // Bogacki-Shampine's RK 5(4) non-FSAL interpolation method
 29 // Definition of the stepper() method that evaluates one step in
 30 // field propagation.
 31 //
 32 // The Butcher table of the Bogacki-Shampine-8-4-5 method is:
 33 //
 34 //   0  |
 35 //  1/6 |      1/6
 36 //  2/9 |      2/27          4/27
 37 //  3/7 |    183/1372     -162/343      1053/1372
 38 //  2/3 |     68/297        -4/11         42/143         1960/3861
 39 //  3/4 |    597/22528      81/352     63099/585728     58653/366080    4617/20480
 40 //   1  | 174197/959244 -30942/79937 8152137/19744439  666106/1039181 -29421/29068  482048/414219
 41 //   1  |    587/8064         0      4440339/15491840   24353/124800     387/44800    2152/5985   7267/94080
 42 //-------------------------------------------------------------------------------------------------------------------
 43 //           587/8064         0      4440339/15491840   24353/124800     387/44800    2152/5985   7267/94080       0
 44 //          2479/34992        0          123/416       612941/3411720    43/1440      2272/6561  79937/1113912  3293/556956
 45 //
 46 // Coefficients have been obtained from:
 47 // http://www.netlib.org/ode/rksuite/
 48 //
 49 // Note on meaning of label "non-FSAL version": 
 50 //   This method calculates the deriviative dy/dx at the endpoint of the
 51 //   integration interval at each step, as part of its evaluation of the
 52 //   endpoint and its error. So this value is available to be returned,
 53 //   for re-use in case of a successful step.
 54 //   (This is done in a 'later' version using a refined interface).
 55 //
 56 // Created: Somnath Banerjee, Google Summer of Code 2015, May-August 2015
 57 // Revision: John Apostolakis, CERN, May 2016
 58 // --------------------------------------------------------------------
 59 
 60 #include <cassert>
 61 
 62 #include "G4BogackiShampine45.hh"
 63 #include "G4LineSection.hh"
 64 
 65 G4bool G4BogackiShampine45::fPreparedConstants = false;
 66 G4double G4BogackiShampine45::bi[12][7];
 67 
 68 // Constructor
 69 //
 70 G4BogackiShampine45::G4BogackiShampine45(G4EquationOfMotion *EqRhs,
 71                                          G4int     noIntegrationVariables,
 72                                          G4bool    primary)   
 73    : G4MagIntegratorStepper(EqRhs, noIntegrationVariables)
 74 {
 75     const G4int numberOfVariables = noIntegrationVariables;
 76     
 77     // New Chunk of memory being created for use by the stepper
 78     
 79     // aki - for storing intermediate RHS
 80     ak2 = new G4double[numberOfVariables];
 81     ak3 = new G4double[numberOfVariables];
 82     ak4 = new G4double[numberOfVariables];
 83     ak5 = new G4double[numberOfVariables];
 84     ak6 = new G4double[numberOfVariables];
 85     ak7 = new G4double[numberOfVariables];
 86     ak8 = new G4double[numberOfVariables];
 87     ak9 = new G4double[numberOfVariables];
 88     ak10 = new G4double[numberOfVariables];
 89     ak11 = new G4double[numberOfVariables];    
 90 
 91     for (auto & i : p)
 92     {
 93        i= new G4double[numberOfVariables];
 94     }
 95 
 96     assert ( GetNumberOfStateVariables() >= 8 );    
 97     const G4int numStateVars = std::max(noIntegrationVariables,
 98                                         GetNumberOfStateVariables() );  
 99 
100     // Must ensure space extra 'state' variables exists - i.e. yIn[7]
101     yTemp = new G4double[numStateVars];
102     yIn = new G4double[numStateVars] ;
103     
104     fLastInitialVector = new G4double[numStateVars] ;
105     fLastFinalVector = new G4double[numStateVars] ;
106     fLastDyDx = new G4double[numberOfVariables];  // Only derivatives
107 
108     fMidVector = new G4double[numberOfVariables];
109     fMidError =  new G4double[numberOfVariables];
110 
111     if( ! fPreparedConstants )
112     {
113        PrepareConstants();
114     }
115     
116     if( primary )
117     {
118        fAuxStepper = new G4BogackiShampine45(EqRhs, numberOfVariables, false);
119     }
120 }
121 
122 // Destructor
123 //
124 G4BogackiShampine45::~G4BogackiShampine45()
125 {
126     // Clear all previously allocated memory for stepper and DistChord
127     //
128     delete [] ak2;
129     delete [] ak3;
130     delete [] ak4;
131     delete [] ak5;
132     delete [] ak6;
133     delete [] ak7;
134     delete [] ak8;
135     delete [] ak9;
136     delete [] ak10;
137     delete [] ak11;
138 
139     for (auto & i : p)
140     {
141        delete [] i;
142     }
143 
144     delete [] yTemp;
145     delete [] yIn;
146     
147     delete [] fLastInitialVector;
148     delete [] fLastFinalVector;
149     delete [] fLastDyDx;
150     delete [] fMidVector;
151     delete [] fMidError;
152     
153     delete fAuxStepper;
154 }
155 
156 void G4BogackiShampine45::GetLastDydx( G4double dyDxLast[] )
157 {
158    const G4int numberOfVariables = GetNumberOfVariables();
159    
160    for(G4int i=0; i < numberOfVariables; ++i )
161    {
162       dyDxLast[i] = ak9[i];
163    }
164 }       
165 
166 // Stepper
167 //
168 // Passing in the value of yInput[],the first time dydx[] and Step length
169 // Giving back yOut and yErr arrays for output and error respectively
170 //
171 void G4BogackiShampine45::Stepper( const G4double yInput[],
172                                    const G4double DyDx[],
173                                          G4double Step,
174                                          G4double yOut[],
175                                          G4double yErr[] )
176 {
177     G4int i;
178     
179     // Constants from the Butcher tableu
180     //
181     const G4double 
182        b21 = 1.0/6.0 ,
183        b31 = 2.0/27.0 ,     b32 = 4.0/27.0,
184        
185        b41 = 183.0/1372.0 , b42 = -162.0/343.0, b43 = 1053.0/1372.0,
186        
187        b51 = 68.0/297.0,    b52 = -4.0/11.0,
188        b53 = 42.0/143.0,    b54 = 1960.0/3861.0,
189        
190        b61 = 597.0/22528.0,    b62 = 81.0/352.0,
191        b63 = 63099.0/585728.0, b64 = 58653.0/366080.0,
192        b65 = 4617.0/20480.0,
193        
194        b71 = 174197.0/959244.0,    b72 = -30942.0/79937.0,
195        b73 = 8152137.0/19744439.0, b74 = 666106.0/1039181.0,
196        b75 = -29421.0/29068.0,     b76 = 482048.0/414219.0,
197        
198        b81 = 587.0/8064.0,         b82 = 0.0,
199        b83 = 4440339.0/15491840.0, b84 = 24353.0/124800.0,
200        b85 = 387.0/44800.0,        b86 = 2152.0/5985.0,
201        b87 = 7267.0/94080.0;
202     
203 //    c1 = 2479.0/34992.0,
204 //    c2 = 0.0,
205 //    c3 = 123.0/416.0,
206 //    c4 = 612941.0/3411720.0,
207 //    c5 = 43.0/1440.0,
208 //    c6 = 2272.0/6561.0,
209 //    c7 = 79937.0/1113912.0,
210 //    c8 = 3293.0/556956.0,
211     
212     // For the embedded higher order method only the difference of values
213     // taken and is used directly later (instead of defining the last row
214     // of Butcher table in separate constants and taking the
215     // difference)
216     //
217     const G4double 
218        dc1 = b81 -   2479.0 /   34992.0 ,
219        dc2 = 0.0,
220        dc3 = b83 -    123.0 /     416.0 ,
221        dc4 = b84 - 612941.0 / 3411720.0,
222        dc5 = b85 -     43.0 /    1440.0,
223        dc6 = b86 -   2272.0 /    6561.0,
224        dc7 = b87 -  79937.0 / 1113912.0,
225        dc8 =     -   3293.0 /  556956.0;
226 
227     const G4int numberOfVariables = GetNumberOfVariables();
228     
229     // The number of variables to be integrated over
230     //
231     yOut[7] = yTemp[7]  = yIn[7] = yInput[7];
232 
233     // Saving yInput because yInput and yOut can be aliases for same array
234     //
235     for(i=0; i<numberOfVariables; ++i)
236     {
237         yIn[i]=yInput[i];
238     }
239     
240     // RightHandSide(yIn, dydx) ;
241     // 1st Step - Not doing, getting passed
242     //
243     for(i=0; i<numberOfVariables; ++i)
244     {
245         yTemp[i] = yIn[i] + b21*Step*DyDx[i] ;
246     }
247     RightHandSide(yTemp, ak2) ;              // 2nd Step
248     
249     for(i=0; i<numberOfVariables; ++i)
250     {
251         yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + b32*ak2[i]) ;
252     }
253     RightHandSide(yTemp, ak3) ;              // 3rd Step
254     
255     for(i=0; i<numberOfVariables; ++i)
256     {
257         yTemp[i] = yIn[i] + Step*(b41*DyDx[i] + b42*ak2[i] + b43*ak3[i]) ;
258     }
259     RightHandSide(yTemp, ak4) ;              // 4th Step
260     
261     for(i=0; i<numberOfVariables; ++i)
262     {
263         yTemp[i] = yIn[i] + Step*(b51*DyDx[i] + b52*ak2[i] + b53*ak3[i] +
264                                   b54*ak4[i]) ;
265     }
266     RightHandSide(yTemp, ak5) ;              // 5th Step
267     
268     for(i=0; i<numberOfVariables; ++i)
269     {
270         yTemp[i] = yIn[i] + Step*(b61*DyDx[i] + b62*ak2[i] + b63*ak3[i] +
271                                   b64*ak4[i] + b65*ak5[i]) ;
272     }
273     RightHandSide(yTemp, ak6) ;              // 6th Step
274     
275     for(i=0; i<numberOfVariables; ++i)
276     {
277         yTemp[i] = yIn[i] + Step*(b71*DyDx[i] + b72*ak2[i] + b73*ak3[i] +
278                                   b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
279     }
280     RightHandSide(yTemp, ak7);               // 7th Step
281     
282     for(i=0; i<numberOfVariables; ++i)
283     {
284         yOut[i] = yIn[i] + Step*(b81*DyDx[i] + b82*ak2[i] + b83*ak3[i] +
285                                   b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
286                                   b87*ak7[i]);
287     }
288     RightHandSide(yOut, ak8);                // 8th Step - Final one Using FSAL
289 
290     for(i=0; i<numberOfVariables; ++i)
291     {
292         yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
293                         dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] + dc8*ak8[i]) ;
294         
295         // Store Input and Final values, for possible use in calculating chord
296         //
297         fLastInitialVector[i] = yIn[i] ;
298         fLastFinalVector[i]   = yOut[i];
299         fLastDyDx[i]          = DyDx[i];
300     }
301     
302     fLastStepLength = Step;
303     fPreparedInterpolation= false;
304 
305     return ;
306 }
307 
308 // DistChord
309 //
310 G4double G4BogackiShampine45::DistChord() const
311 {
312     G4double distLine, distChord;
313     G4ThreeVector initialPoint, finalPoint, midPoint;
314     
315     // Store last initial and final points
316     // (they will be overwritten in self-Stepper call!)
317     //
318     initialPoint = G4ThreeVector(fLastInitialVector[0],
319                                  fLastInitialVector[1], fLastInitialVector[2]);
320     finalPoint   = G4ThreeVector(fLastFinalVector[0],
321                                  fLastFinalVector[1],  fLastFinalVector[2]);
322 
323 #if 1
324     // Old method -- Do half a step using StepNoErr
325     //
326     fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5*fLastStepLength,
327                           fMidVector, fMidError);
328 #else
329     // New method -- Using interpolation,
330     // requires only 3 extra stages (ie 3 extra field evaluations )
331 
332     // Use Interpolation, instead of auxiliary stepper to evaluate midpoint
333     //
334     if( ! fPreparedInterpolation )
335     {
336        G4BogackiShampine45* cThis = const_cast<G4BogackiShampine45 *>(this);
337        cThis-> SetupInterpolationHigh();
338     }
339     // For calculating the output at the tau fraction of Step
340     //
341     G4double tau = 0.5;
342     InterpolateHigh( tau, fMidVector );    
343 #endif
344    
345     midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
346     
347     // Use stored values of Initial and Endpoint + new Midpoint to evaluate
348     // distance of Chord
349     
350     if (initialPoint != finalPoint)
351     {
352         distLine  = G4LineSection::Distline( midPoint,initialPoint,finalPoint );
353         distChord = distLine;
354     }
355     else
356     {
357         distChord = (midPoint-initialPoint).mag();
358     }
359     return distChord;
360 }
361 
362 void G4BogackiShampine45::SetupInterpolationHigh()
363 {
364     // Coefficients for the additional stages
365     //
366     const G4double
367     a91 = 455.0/6144.0 ,
368     a92 = 0.0 ,
369     a93 = 10256301.0/35409920.0 ,
370     a94 = 2307361.0/17971200.0 ,
371     a95 = -387.0/102400.0 ,
372     a96 = 73.0/5130.0 ,
373     a97 = -7267.0/215040.0 ,
374     a98 = 1.0/32.0 ,
375     
376     a101 = -837888343715.0/13176988637184.0 ,
377     a102 = 30409415.0/52955362.0 ,
378     a103 = -48321525963.0/759168069632.0 ,
379     a104 = 8530738453321.0/197654829557760.0 ,
380     a105 = 1361640523001.0/1626788720640.0 ,
381     a106 = -13143060689.0/38604458898.0 ,
382     a107 = 18700221969.0/379584034816.0 ,
383     a108 = -5831595.0/847285792.0 ,
384     a109 = -5183640.0/26477681.0 ,
385     
386     a111 = 98719073263.0/1551965184000.0 ,
387     a112 = 1307.0/123552.0 ,
388     a113 = 4632066559387.0/70181753241600.0 ,
389     a114 = 7828594302389.0/382182512025600.0 ,
390     a115 = 40763687.0/11070259200.0 ,
391     a116 = 34872732407.0/224610586200.0 ,
392     a117 = -2561897.0/30105600.0 ,
393     a118 = 1.0/10.0 ,
394     a119 = -1.0/10.0 ,
395     a1110 = -1403317093.0/11371610250.0 ;
396     
397     const G4int numberOfVariables= this->GetNumberOfVariables();
398     const G4double* dydx= fLastDyDx;
399     const G4double Step = fLastStepLength;
400     
401     yTemp[7]  = yIn[7];
402 
403     // Evaluate the extra stages
404     //
405     for(G4int i=0; i<numberOfVariables; ++i)
406     {
407         yTemp[i] = yIn[i] + Step*(a91*dydx[i] + a92*ak2[i] + a93*ak3[i] +
408                                   a94*ak4[i] + a95*ak5[i] + a96*ak6[i] +
409                                   a97*ak7[i] + a98*ak8[i] );
410     }
411     
412     RightHandSide(yTemp, ak9);   // 9th stage
413     
414     for(G4int i=0; i<numberOfVariables; ++i)
415     {
416         yTemp[i] = yIn[i] + Step*(a101*dydx[i] + a102*ak2[i] + a103*ak3[i] +
417                                   a104*ak4[i] + a105*ak5[i] + a106*ak6[i] +
418                                   a107*ak7[i] + a108*ak8[i] + a109*ak9[i] );
419     }
420     
421     RightHandSide(yTemp, ak10);  // 10th stage
422     
423     for(G4int i=0; i<numberOfVariables; ++i)
424     {
425         yTemp[i] = yIn[i] + Step*(a111*dydx[i] + a112*ak2[i] + a113*ak3[i] +
426                                   a114*ak4[i] + a115*ak5[i] + a116*ak6[i] +
427                                   a117*ak7[i] + a118*ak8[i] + a119*ak9[i] +
428                                   a1110*ak10[i] );
429     }
430     RightHandSide(yTemp, ak11);  // 11th stage
431 
432     //  In future we can restrict the number of variables interpolated
433     //
434     G4int nwant = numberOfVariables; 
435 
436     //  Form the coefficients of the interpolating polynomial in its shifted
437     //  and scaled form.  The terms are grouped to minimize the errors 
438     //  of the transformation, to cope with ill-conditioning. ( From RKSUITE )
439     //   
440     for (G4int l = 0; l < nwant; ++l)
441     {
442         //  Coefficient of tau^6        
443         p[5][l] =   bi[5][6]*ak5[l] +
444                   ((bi[10][6]*ak10[l] + bi[8][6]*ak8[l]) + 
445                    (bi[7][6]*ak7[l] + bi[6][6]*ak6[l]))  + 
446                   ((bi[4][6]*ak4[l] + bi[9][6]*ak9[l]) + 
447                    (bi[3][6]*ak3[l] + bi[11][6]*ak11[l]) + 
448                     bi[1][6]*dydx[l]);
449         //  Coefficient of tau^5
450         p[4][l] = (bi[10][5]*ak10[l] + bi[9][5]*ak9[l])  + 
451                  ((bi[7][5]*ak7[l] + bi[6][5]*ak6[l]) + 
452                    bi[5][5]*ak5[l])  +  ((bi[4][5]*ak4[l] + 
453                    bi[8][5]*ak8[l]) + (bi[3][5]*ak3[l] + 
454                    bi[11][5]*ak11[l]) + bi[1][5]*dydx[l]);
455         //  Coefficient of tau^4
456         p[3][l] = ((bi[4][4]*ak4[l] + bi[8][4]*ak8[l]) + 
457                    (bi[7][4]*ak7[l] + bi[6][4]*ak6[l]) + 
458                     bi[5][4]*ak5[l]) + ((bi[10][4]*ak10[l] + 
459                     bi[9][4]*ak9[l]) +  (bi[3][4]*ak3[l] + 
460                     bi[11][4]*ak11[l]) + bi[1][4]*dydx[l]);
461         //  Coefficient of tau^3
462         p[2][l] =  bi[5][3]*ak5[l] + bi[6][3]*ak6[l]  + 
463                  ((bi[3][3]*ak3[l] + bi[9][3]*ak9[l]) + 
464                  (bi[10][3]*ak10[l]+ bi[8][3]*ak8[l]) + bi[1][3]*dydx[l]) + 
465                  ((bi[4][3]*ak4[l] + bi[11][3]*ak11[l]) + bi[7][3]*ak7[l]);
466         //  Coefficient of tau^2
467         p[1][l] = bi[5][2]*ak5[l]  + ((bi[6][2]*ak6[l] + 
468                   bi[8][2]*ak8[l]) +   bi[1][2]*dydx[l])  + 
469                 ((bi[3][2]*ak3[l]  +   bi[9][2]*ak9[l]) + 
470                  bi[10][2]*ak10[l])+ ((bi[4][2]*ak4[l] + 
471                  bi[11][2]*ak2[l]) +   bi[7][2]*ak7[l]);
472       }
473 
474       //  Scale all the coefficients by the step size.
475       //
476       for (auto & i : p)
477       {
478          for (G4int l = 0; l < nwant; ++l)
479          {  
480             i[l] *= Step;
481          }
482       }
483 
484     fPreparedInterpolation = true;
485 }
486 
487 void G4BogackiShampine45::PrepareConstants()
488 {
489     for(auto i=1; i<= 11; ++i)
490     {
491         bi[i][1] = 0.0 ;
492     }
493     
494     for(auto i=1; i<=6; ++i)
495     {
496         bi[2][i] = 0.0 ;
497     }
498 
499     bi[1][6] = -12134338393.0 / 1050809760.0 ,
500     bi[1][5] =  -1620741229.0 / 50038560.0 ,
501     bi[1][4] =  -2048058893.0 / 59875200.0 ,
502     bi[1][3] = -87098480009.0 / 5254048800.0 ,
503     bi[1][2] = -11513270273.0 / 3502699200.0 ,
504     // 
505     bi[3][6] =  -33197340367.0 / 1218433216.0 ,
506     bi[3][5] = -539868024987.0 / 6092166080.0 ,
507     bi[3][4] =  -39991188681.0 / 374902528.0 ,
508     bi[3][3] =  -69509738227.0 / 1218433216.0 ,
509     bi[3][2] =  -29327744613.0 / 2436866432.0 ,
510     //
511     bi[4][6] =   -284800997201.0 /  19905339168.0 ,
512     bi[4][5] =  -7896875450471.0 / 165877826400.0 ,
513     bi[4][4] =   -333945812879.0 /   5671036800.0 ,
514     bi[4][3] = -16209923456237.0 / 497633479200.0 ,
515     bi[4][2] =  -2382590741699.0 / 331755652800.0 ,
516     //
517     bi[5][6] = -540919.0 / 741312.0 ,
518     bi[5][5] = -103626067.0 / 43243200.0 ,
519     bi[5][4] = -633779.0 / 211200.0 ,
520     bi[5][3] = -32406787.0 / 18532800.0 ,
521     bi[5][2] = -36591193.0 / 86486400.0 ,
522     //
523     bi[6][6] = 7157998304.0 / 374350977.0 ,
524     bi[6][5] = 30405842464.0 / 623918295.0 ,
525     bi[6][4] = 183022264.0 / 5332635.0 ,
526     bi[6][3] = -3357024032.0 / 1871754885.0 ,
527     bi[6][2] = -611586736.0 / 89131185.0 ,
528     //
529     bi[7][6] =  -138073.0 /  9408.0 ,
530     bi[7][5] =  -719433.0 / 15680.0 ,
531     bi[7][4] = -1620541.0 / 31360.0 ,
532     bi[7][3] =  -385151.0 / 15680.0 ,
533     bi[7][2] =  -65403.0  / 15680.0 ,
534     //
535     bi[8][6] = 1245.0 / 64.0 ,
536     bi[8][5] = 3991.0 / 64.0 ,
537     bi[8][4] = 4715.0 / 64.0 ,
538     bi[8][3] = 2501.0 / 64.0 ,
539     bi[8][2] =  149.0 / 16.0 ,
540     bi[8][1] = 1.0 ,
541     //
542     bi[9][6] = 55.0 / 3.0 ,
543     bi[9][5] = 71.0 ,
544     bi[9][4] = 103.0 ,
545     bi[9][3] = 199.0 / 3.0 ,
546     bi[9][2] = 16.0 ,
547     //
548     bi[10][6] =  -1774004627.0  /  75810735.0 ,
549     bi[10][5] =  -1774004627.0  /  25270245.0 ,
550     bi[10][4] =    -26477681.0  /    359975.0 ,
551     bi[10][3] = -11411880511.0  / 379053675.0 ,
552     bi[10][2] =   -423642896.0  / 126351225.0 ,
553     //
554     bi[11][6] =  35.0 ,
555     bi[11][5] = 105.0 ,
556     bi[11][4] = 117.0 ,
557     bi[11][3] =  59.0 ,
558     bi[11][2] =  12.0 ;
559 
560     fPreparedConstants = true;
561 }
562 
563 void G4BogackiShampine45::InterpolateHigh(G4double tau, G4double* yOut) const
564 {
565     G4int numberOfVariables = GetNumberOfVariables();
566     
567     G4Exception("G4BogackiShampine45::InterpolateHigh()", "GeomField0001",
568                 FatalException, "Method is not yet validated.");
569 
570     // const G4double *yIn=  fLastInitialVector;
571     // const G4double *dydx= fLastDyDx;
572     const G4double Step = fLastStepLength;
573     
574 #if 1
575     G4int nwant = numberOfVariables;
576     const G4int norder= 6;
577     G4int l, k;
578 
579     for (l = 0; l < nwant; ++l)
580     {
581       yOut[l] = p[norder-1][l] * tau;
582     }
583     for (k = norder - 2; k >= 1; --k)
584     {
585       for (l = 0; l < nwant; ++l)
586       {
587          yOut[l] = ( yOut[l] + p[k][l] ) * tau;
588       }
589     }
590     for (l = 0; l < nwant; ++l)
591     {
592       yOut[l] = ( yOut[l] + Step * ak8[l] ) * tau + yIn[l];
593     }
594     // The derivative at the end-point is nextDydx[i] = ak8[i];
595 #else
596     // The scheme tries to do the same as the DormandPrince745 routine,
597     // but fails 
598 
599     G4double b[12];
600     const G4double* dydx = fLastDyDx;
601     
602     G4double tau0 = tau;
603     
604     for(G4int iStage=1; iStage<=11; ++iStage)   // iStage = stage number
605     {
606         b[iStage] = 0.0;
607         tau = tau0;
608         for(G4int j=6; j>=1; --j)               // j reversed
609         {
610             b[iStage] += bi[iStage][j] * tau;
611             tau *= tau0;
612         }
613     }
614 
615     for(G4int i=0; i<numberOfVariables; ++i)
616     {
617         yOut[i] = yIn[i] + Step*(b[1]*dydx[i] + b[2]*ak2[i] + b[3]*ak3[i] +
618                                  b[4]*ak4[i] + b[5]*ak5[i] + b[6]*ak6[i] +
619                                  b[7]*ak7[i] + b[8]*ak8[i] + b[9]*ak9[i] +
620                                  b[10]*ak10[i] + b[11]*ak11[i] );
621     }
622 #endif    
623 }
624