Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/magneticfield/src/G4FSALBogackiShampine45.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/G4FSALBogackiShampine45.cc (Version 11.3.0) and /geometry/magneticfield/src/G4FSALBogackiShampine45.cc (Version 11.0.p4)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // G4FSALBogackiShampine45 implementation          26 // G4FSALBogackiShampine45 implementation
 27 //                                                 27 //
 28 // The Butcher table of the Bogacki-Shampine-8     28 // The Butcher table of the Bogacki-Shampine-8-4-5 method is as follows:
 29 //                                                 29 //
 30 // 0   |                                           30 // 0   |
 31 // 1/6 | 1/6                                       31 // 1/6 | 1/6
 32 // 2/9 | 2/27          4/27                        32 // 2/9 | 2/27          4/27
 33 // 3/7 | 183/1372      -162/343     1053/1372      33 // 3/7 | 183/1372      -162/343     1053/1372
 34 // 2/3 | 68/297        -4/11        42/143 196     34 // 2/3 | 68/297        -4/11        42/143 1960/3861
 35 // 3/4 | 597/22528     81/352       63099/5857     35 // 3/4 | 597/22528     81/352       63099/585728     58653/366080   4617/20480
 36 // 1   | 174197/959244 -30942/79937 8152137/19     36 // 1   | 174197/959244 -30942/79937 8152137/19744439 666106/1039181 -29421/29068 482048/414219
 37 // 1   | 587/8064      0            4440339/15     37 // 1   | 587/8064      0            4440339/15491840 24353/124800   387/44800    2152/5985     7267/94080
 38 // -------------------------------------------     38 // -------------------------------------------------------------------------------------------------------------------
 39 //       587/8064      0            4440339/15     39 //       587/8064      0            4440339/15491840 24353/124800    387/44800   2152/5985     7267/94080     0
 40 //       2479/34992    0            123/416        40 //       2479/34992    0            123/416          612941/3411720  43/1440     2272/6561     79937/1113912  3293/556956
 41 //                                                 41 //
 42 // Created: Somnath Banerjee, Google Summer of     42 // Created: Somnath Banerjee, Google Summer of Code 2015, 26 May 2015
 43 // Supervision: John Apostolakis, CERN             43 // Supervision: John Apostolakis, CERN
 44 // -------------------------------------------     44 // --------------------------------------------------------------------
 45                                                    45 
 46 // Plan is that this source file / class will      46 // Plan is that this source file / class will be merged with the updated
 47 // BogackiShampine45 class, which contains imp     47 // BogackiShampine45 class, which contains improvements (May 2016) 
 48                                                    48 
 49 #include <cassert>                                 49 #include <cassert>
 50                                                    50 
 51 #include "G4FSALBogackiShampine45.hh"              51 #include "G4FSALBogackiShampine45.hh"
 52 #include "G4LineSection.hh"                        52 #include "G4LineSection.hh"
 53                                                    53 
 54 G4bool   G4FSALBogackiShampine45::fPreparedCon     54 G4bool   G4FSALBogackiShampine45::fPreparedConstants = false;
 55 G4double G4FSALBogackiShampine45::bi[12][7];       55 G4double G4FSALBogackiShampine45::bi[12][7];
 56                                                    56 
 57 // Constructor                                     57 // Constructor
 58 //                                                 58 //
 59 G4FSALBogackiShampine45::G4FSALBogackiShampine     59 G4FSALBogackiShampine45::G4FSALBogackiShampine45(G4EquationOfMotion* EqRhs,
 60                                                    60                                                  G4int noIntegrationVariables,
 61                                                    61                                                  G4bool primary)
 62    : G4VFSALIntegrationStepper(EqRhs, noIntegr     62    : G4VFSALIntegrationStepper(EqRhs, noIntegrationVariables)
 63 {                                                  63 {
 64     const G4int numberOfVariables = noIntegrat     64     const G4int numberOfVariables = noIntegrationVariables;
 65                                                    65     
 66     // New Chunk of memory being created for u     66     // New Chunk of memory being created for use by the stepper
 67                                                    67     
 68     // aki - for storing intermediate RHS          68     // aki - for storing intermediate RHS
 69     //                                             69     //
 70     ak2 = new G4double[numberOfVariables];         70     ak2 = new G4double[numberOfVariables];
 71     ak3 = new G4double[numberOfVariables];         71     ak3 = new G4double[numberOfVariables];
 72     ak4 = new G4double[numberOfVariables];         72     ak4 = new G4double[numberOfVariables];
 73     ak5 = new G4double[numberOfVariables];         73     ak5 = new G4double[numberOfVariables];
 74     ak6 = new G4double[numberOfVariables];         74     ak6 = new G4double[numberOfVariables];
 75     ak7 = new G4double[numberOfVariables];         75     ak7 = new G4double[numberOfVariables];
 76     ak8 = new G4double[numberOfVariables];         76     ak8 = new G4double[numberOfVariables];
 77                                                    77 
 78     ak9 = new G4double[numberOfVariables];         78     ak9 = new G4double[numberOfVariables];
 79     ak10 = new G4double[numberOfVariables];        79     ak10 = new G4double[numberOfVariables];
 80     ak11 = new G4double[numberOfVariables];        80     ak11 = new G4double[numberOfVariables];
 81     DyDx = new G4double[numberOfVariables];        81     DyDx = new G4double[numberOfVariables];
 82                                                    82     
 83     assert ( GetNumberOfStateVariables() >= 8      83     assert ( GetNumberOfStateVariables() >= 8 );
 84     const G4int numStateVars = std::max(noInte     84     const G4int numStateVars = std::max(noIntegrationVariables,
 85                                         GetNum     85                                         GetNumberOfStateVariables() );  
 86                                                    86 
 87     // Must ensure space extra 'state' variabl     87     // Must ensure space extra 'state' variables exists - i.e. yIn[7]
 88     //                                             88     //
 89     yTemp = new G4double[numStateVars];            89     yTemp = new G4double[numStateVars];
 90     yIn = new G4double[numStateVars] ;             90     yIn = new G4double[numStateVars] ;
 91                                                    91     
 92     fLastInitialVector = new G4double[numState     92     fLastInitialVector = new G4double[numStateVars] ;
 93     fLastFinalVector = new G4double[numStateVa     93     fLastFinalVector = new G4double[numStateVars] ;
 94     fLastDyDx = new G4double[numberOfVariables     94     fLastDyDx = new G4double[numberOfVariables];  // Only derivatives
 95                                                    95     
 96     fMidVector = new G4double[numStateVars];       96     fMidVector = new G4double[numStateVars];
 97     fMidError =  new G4double[numStateVars];       97     fMidError =  new G4double[numStateVars];
 98                                                    98     
 99     pseudoDydx_for_DistChord = new G4double[nu     99     pseudoDydx_for_DistChord = new G4double[numberOfVariables];
100                                                   100     
101     fMidVector = new G4double[numberOfVariable    101     fMidVector = new G4double[numberOfVariables];
102     fMidError =  new G4double[numberOfVariable    102     fMidError =  new G4double[numberOfVariables];
103     if( primary )                                 103     if( primary )
104     {                                             104     {
105       fAuxStepper = new G4FSALBogackiShampine4    105       fAuxStepper = new G4FSALBogackiShampine45(EqRhs, numberOfVariables,
106                                                   106                                                 !primary);
107     }                                             107     }
108     if( !fPreparedConstants )                     108     if( !fPreparedConstants )
109     {                                             109     {
110        PrepareConstants();                        110        PrepareConstants();
111     }                                             111     }
112 }                                                 112 }
113                                                   113 
114 // Destructor                                     114 // Destructor
115 //                                                115 //
116 G4FSALBogackiShampine45::~G4FSALBogackiShampin    116 G4FSALBogackiShampine45::~G4FSALBogackiShampine45()
117 {                                                 117 {
118     // Clear all previously allocated memory f    118     // Clear all previously allocated memory for stepper and DistChord
119                                                   119 
120     delete [] ak2;                                120     delete [] ak2;
121     delete [] ak3;                                121     delete [] ak3;
122     delete [] ak4;                                122     delete [] ak4;
123     delete [] ak5;                                123     delete [] ak5;
124     delete [] ak6;                                124     delete [] ak6;
125     delete [] ak7;                                125     delete [] ak7;
126     delete [] ak8;                                126     delete [] ak8;
127     delete [] ak9;                                127     delete [] ak9;
128     delete [] ak10;                               128     delete [] ak10;
129     delete [] ak11;                               129     delete [] ak11;
130     delete [] DyDx;                               130     delete [] DyDx;    
131     delete [] yTemp;                              131     delete [] yTemp;
132     delete [] yIn;                                132     delete [] yIn;
133                                                   133     
134     delete [] fLastInitialVector;                 134     delete [] fLastInitialVector;
135     delete [] fLastFinalVector;                   135     delete [] fLastFinalVector;
136     delete [] fLastDyDx;                          136     delete [] fLastDyDx;
137     delete [] fMidVector;                         137     delete [] fMidVector;
138     delete [] fMidError;                          138     delete [] fMidError;
139                                                   139     
140     delete fAuxStepper;                           140     delete fAuxStepper;
141                                                   141     
142     delete [] pseudoDydx_for_DistChord;           142     delete [] pseudoDydx_for_DistChord;
143 }                                                 143 }
144                                                   144 
145 // Stepper                                        145 // Stepper
146 //                                                146 //
147 // Passing in the value of yInput[],the first     147 // Passing in the value of yInput[],the first time dydx[] and Step length
148 // Giving back yOut and yErr arrays for output    148 // Giving back yOut and yErr arrays for output and error respectively
149 //                                                149 //
150 void G4FSALBogackiShampine45::Stepper(const G4    150 void G4FSALBogackiShampine45::Stepper(const G4double yInput[],
151                                       const G4    151                                       const G4double dydx[],
152                                             G4    152                                             G4double Step,
153                                             G4    153                                             G4double yOut[],
154                                             G4    154                                             G4double yErr[],
155                                             G4    155                                             G4double nextDydx[])
156 {                                                 156 {
157     G4int i;                                      157     G4int i;
158                                                   158     
159     // The various constants defined on the ba    159     // The various constants defined on the basis of butcher tableu
160                                                   160 
161     const G4double b21 = 1.0/6.0 ,                161     const G4double b21 = 1.0/6.0 ,
162                    b31 = 2.0/27.0 , b32 = 4.0/    162                    b31 = 2.0/27.0 , b32 = 4.0/27.0,
163                                                   163     
164                    b41 = 183.0/1372.0 , b42 =     164                    b41 = 183.0/1372.0 , b42 = -162.0/343.0, b43 = 1053.0/1372.0,
165                                                   165     
166                    b51 = 68.0/297.0, b52 = -4.    166                    b51 = 68.0/297.0, b52 = -4.0/11.0,
167                    b53 = 42.0/143.0, b54 = 196    167                    b53 = 42.0/143.0, b54 = 1960.0/3861.0,
168                                                   168     
169                    b61 = 597.0/22528.0, b62 =     169                    b61 = 597.0/22528.0, b62 = 81.0/352.0,
170                    b63 = 63099.0/585728.0, b64    170                    b63 = 63099.0/585728.0, b64 = 58653.0/366080.0,
171                    b65 = 4617.0/20480.0,          171                    b65 = 4617.0/20480.0,
172                                                   172     
173                    b71 = 174197.0/959244.0, b7    173                    b71 = 174197.0/959244.0, b72 = -30942.0/79937.0,
174                    b73 = 8152137.0/19744439.0,    174                    b73 = 8152137.0/19744439.0, b74 = 666106.0/1039181.0,
175                    b75 = -29421.0/29068.0,  b7    175                    b75 = -29421.0/29068.0,  b76 = 482048.0/414219.0,
176                                                   176     
177                    b81 = 587.0/8064.0,  b82 =     177                    b81 = 587.0/8064.0,  b82 = 0.0,
178                    b83 = 4440339.0/15491840.0,    178                    b83 = 4440339.0/15491840.0,  b84 = 24353.0/124800.0,
179                    b85 = 387.0/44800.0, b86 =     179                    b85 = 387.0/44800.0, b86 = 2152.0/5985.0,
180                    b87 = 7267.0/94080.0,          180                    b87 = 7267.0/94080.0,
181                                                   181     
182                                                   182     
183              //    c1 = 2479.0/34992.0,           183              //    c1 = 2479.0/34992.0,
184              //    c2 = 0.0,                      184              //    c2 = 0.0,
185              //    c3 = 123.0/416.0,              185              //    c3 = 123.0/416.0,
186              //    c4 = 612941.0/3411720.0,       186              //    c4 = 612941.0/3411720.0,
187              //    c5 = 43.0/1440.0,              187              //    c5 = 43.0/1440.0,
188              //    c6 = 2272.0/6561.0,            188              //    c6 = 2272.0/6561.0,
189              //    c7 = 79937.0/1113912.0,        189              //    c7 = 79937.0/1113912.0,
190              //    c8 = 3293.0/556956.0,          190              //    c8 = 3293.0/556956.0,
191                                                   191     
192     // For the embedded higher order method on    192     // For the embedded higher order method only the difference of values
193     // taken and is used directly later instea    193     // taken and is used directly later instead of defining the last row
194     // of butcher table in a separate set of v    194     // of butcher table in a separate set of variables and taking the
195     // difference there                           195     // difference there
196                                                   196     
197                    dc1 = b81 - 2479.0/34992.0     197                    dc1 = b81 - 2479.0/34992.0 ,
198                    dc2 = 0.0,                     198                    dc2 = 0.0,
199                    dc3 = b83 - 123.0/416.0 ,      199                    dc3 = b83 - 123.0/416.0 ,
200                    dc4 = b84 - 612941.0/341172    200                    dc4 = b84 - 612941.0/3411720.0,
201                    dc5 = b85 - 43.0/1440.0,       201                    dc5 = b85 - 43.0/1440.0,
202                    dc6 = b86 - 2272.0/6561.0,     202                    dc6 = b86 - 2272.0/6561.0,
203                    dc7 = b87 - 79937.0/1113912    203                    dc7 = b87 - 79937.0/1113912.0,
204                    dc8 = -3293.0/556956.0;   /    204                    dc8 = -3293.0/556956.0;   // end of declaration
205                                                   205 
206     const G4int numberOfVariables = GetNumberO    206     const G4int numberOfVariables = GetNumberOfVariables();
207                                                   207     
208     // The number of variables to be integrate    208     // The number of variables to be integrated over
209     //                                            209     //
210     yOut[7] = yTemp[7]  = yIn[7];                 210     yOut[7] = yTemp[7]  = yIn[7];
211                                                   211 
212     //  Saving yInput because yInput and yOut     212     //  Saving yInput because yInput and yOut can be aliases for same array
213     //                                            213     //
214     for(i=0; i<numberOfVariables; ++i)            214     for(i=0; i<numberOfVariables; ++i)
215     {                                             215     {
216         yIn[i]=yInput[i];                         216         yIn[i]=yInput[i];
217         DyDx[i] = dydx[i];                        217         DyDx[i] = dydx[i];
218     }                                             218     }
219     // RightHandSide(yIn, dydx) ;   // 1st Ste    219     // RightHandSide(yIn, dydx) ;   // 1st Step - Not doing, getting passed
220                                                   220     
221     for(i=0; i<numberOfVariables; ++i)            221     for(i=0; i<numberOfVariables; ++i)
222     {                                             222     {
223         yTemp[i] = yIn[i] + b21*Step*DyDx[i] ;    223         yTemp[i] = yIn[i] + b21*Step*DyDx[i] ;
224     }                                             224     }
225     RightHandSide(yTemp, ak2) ;              /    225     RightHandSide(yTemp, ak2) ;              // 2nd Step
226                                                   226     
227     for(i=0; i<numberOfVariables; ++i)            227     for(i=0; i<numberOfVariables; ++i)
228     {                                             228     {
229         yTemp[i] = yIn[i] + Step*(b31*DyDx[i]     229         yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + b32*ak2[i]) ;
230     }                                             230     }
231     RightHandSide(yTemp, ak3) ;              /    231     RightHandSide(yTemp, ak3) ;              // 3rd Step
232                                                   232     
233     for(i=0; i<numberOfVariables; ++i)            233     for(i=0; i<numberOfVariables; ++i)
234     {                                             234     {
235         yTemp[i] = yIn[i] + Step*(b41*DyDx[i]     235         yTemp[i] = yIn[i] + Step*(b41*DyDx[i] + b42*ak2[i] + b43*ak3[i]) ;
236     }                                             236     }
237     RightHandSide(yTemp, ak4) ;              /    237     RightHandSide(yTemp, ak4) ;              // 4th Step
238                                                   238     
239     for(i=0; i<numberOfVariables; ++i)            239     for(i=0; i<numberOfVariables; ++i)
240     {                                             240     {
241         yTemp[i] = yIn[i] + Step*(b51*DyDx[i]     241         yTemp[i] = yIn[i] + Step*(b51*DyDx[i] + b52*ak2[i] + b53*ak3[i] +
242                                   b54*ak4[i])     242                                   b54*ak4[i]) ;
243     }                                             243     }
244     RightHandSide(yTemp, ak5) ;              /    244     RightHandSide(yTemp, ak5) ;              // 5th Step
245                                                   245     
246     for(i=0; i<numberOfVariables; ++i)            246     for(i=0; i<numberOfVariables; ++i)
247     {                                             247     {
248         yTemp[i] = yIn[i] + Step*(b61*DyDx[i]     248         yTemp[i] = yIn[i] + Step*(b61*DyDx[i] + b62*ak2[i] + b63*ak3[i] +
249                                   b64*ak4[i] +    249                                   b64*ak4[i] + b65*ak5[i]) ;
250     }                                             250     }
251     RightHandSide(yTemp, ak6) ;              /    251     RightHandSide(yTemp, ak6) ;              // 6th Step
252                                                   252     
253     for(i=0; i<numberOfVariables; ++i)            253     for(i=0; i<numberOfVariables; ++i)
254     {                                             254     {
255         yTemp[i] = yIn[i] + Step*(b71*DyDx[i]     255         yTemp[i] = yIn[i] + Step*(b71*DyDx[i] + b72*ak2[i] + b73*ak3[i] +
256                                   b74*ak4[i] +    256                                   b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
257     }                                             257     }
258     RightHandSide(yTemp, ak7);               /    258     RightHandSide(yTemp, ak7);               // 7th Step
259                                                   259     
260     for(i=0; i<numberOfVariables; ++i)            260     for(i=0; i<numberOfVariables; ++i)
261     {                                             261     {
262         yOut[i] = yIn[i] + Step*(b81*DyDx[i] +    262         yOut[i] = yIn[i] + Step*(b81*DyDx[i] + b82*ak2[i] + b83*ak3[i] +
263                                   b84*ak4[i] +    263                                   b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
264                                   b87*ak7[i]);    264                                   b87*ak7[i]);
265     }                                             265     }
266     RightHandSide(yOut, ak8);                /    266     RightHandSide(yOut, ak8);                // 8th Step - Final one Using FSAL
267                                                   267 
268                                                   268     
269     for(i=0; i<numberOfVariables; ++i)            269     for(i=0; i<numberOfVariables; ++i)
270     {                                             270     {
271                                                   271         
272         yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[    272         yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
273                         dc5*ak5[i] + dc6*ak6[i    273                         dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] + dc8*ak8[i]) ;
274                                                   274         
275                                                   275         
276         // FSAL stepper : Must pass the last D    276         // FSAL stepper : Must pass the last DyDx for the next step, here ak8
277         //                                        277         //
278         nextDydx[i] = ak8[i];                     278         nextDydx[i] = ak8[i];
279                                                   279         
280         // Store Input and Final values, for p    280         // Store Input and Final values, for possible use in calculating chord
281         //                                        281         //
282         fLastInitialVector[i] = yIn[i] ;          282         fLastInitialVector[i] = yIn[i] ;
283         fLastFinalVector[i]   = yOut[i];          283         fLastFinalVector[i]   = yOut[i];
284         fLastDyDx[i]          = DyDx[i];          284         fLastDyDx[i]          = DyDx[i];
285     }                                             285     }
286     fLastStepLength = Step;                       286     fLastStepLength = Step;
287                                                   287     
288     return;                                       288     return;
289 }                                                 289 }
290                                                   290 
291 // DistChord                                      291 // DistChord
292 //                                                292 //
293 G4double G4FSALBogackiShampine45::DistChord()     293 G4double G4FSALBogackiShampine45::DistChord() const
294 {                                                 294 {
295     G4double distLine, distChord;                 295     G4double distLine, distChord;
296     G4ThreeVector initialPoint, finalPoint, mi    296     G4ThreeVector initialPoint, finalPoint, midPoint;
297                                                   297 
298     // Store last initial and final points        298     // Store last initial and final points
299     // (they will be overwritten in self-Stepp    299     // (they will be overwritten in self-Stepper call!)
300     //                                            300     //
301     initialPoint = G4ThreeVector( fLastInitial    301     initialPoint = G4ThreeVector( fLastInitialVector[0],
302                                  fLastInitialV    302                                  fLastInitialVector[1], fLastInitialVector[2]);
303     finalPoint   = G4ThreeVector( fLastFinalVe    303     finalPoint   = G4ThreeVector( fLastFinalVector[0],
304                                  fLastFinalVec    304                                  fLastFinalVector[1],  fLastFinalVector[2]);
305                                                   305     
306     // Do half a step using StepNoErr             306     // Do half a step using StepNoErr
307                                                   307     
308     fAuxStepper->Stepper( fLastInitialVector,     308     fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
309                           fMidVector, fMidErro    309                           fMidVector, fMidError, pseudoDydx_for_DistChord );
310                                                   310     
311     midPoint = G4ThreeVector( fMidVector[0], f    311     midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2] );
312                                                   312     
313     // Use stored values of Initial and Endpoi    313     // Use stored values of Initial and Endpoint + new Midpoint to evaluate
314     // distance of Chord                          314     // distance of Chord
315     //                                            315     //
316     if (initialPoint != finalPoint)               316     if (initialPoint != finalPoint)
317     {                                             317     {
318         distLine = G4LineSection::Distline(mid    318         distLine = G4LineSection::Distline(midPoint, initialPoint, finalPoint);
319         distChord = distLine;                     319         distChord = distLine;
320     }                                             320     }
321     else                                          321     else
322     {                                             322     {
323         distChord = (midPoint-initialPoint).ma    323         distChord = (midPoint-initialPoint).mag();
324     }                                             324     }
325     return distChord;                             325     return distChord;
326 }                                                 326 }
327                                                   327 
328 // PrepareConstants                               328 // PrepareConstants
329 //                                                329 //
330 void G4FSALBogackiShampine45::PrepareConstants    330 void G4FSALBogackiShampine45::PrepareConstants()
331 {                                                 331 {
332     //  --------------------------------------    332     //  --------------------------------------------------------
333     //  COEFFICIENTS FOR INTERPOLANT  bi  WITH    333     //  COEFFICIENTS FOR INTERPOLANT  bi  WITH  11  STAGES
334     //  --------------------------------------    334     //  --------------------------------------------------------
335                                                   335     
336     // Initialise all values of G4double bi[12    336     // Initialise all values of G4double bi[12][7]
337     //                                            337     //
338     for(auto i=1; i<12; ++i)                      338     for(auto i=1; i<12; ++i)
339     {                                             339     {
340        for(auto j=1; j<7; ++j)                    340        for(auto j=1; j<7; ++j)
341        {                                          341        {
342           bi[i][j] = 0.0 ;                        342           bi[i][j] = 0.0 ;
343        }                                          343        }
344     }                                             344     }
345                                                   345   
346     bi[1][6] = -12134338393.0/1050809760.0 ,      346     bi[1][6] = -12134338393.0/1050809760.0 ,
347     bi[1][5] = -1620741229.0/50038560.0 ,         347     bi[1][5] = -1620741229.0/50038560.0 ,
348     bi[1][4] = -2048058893.0/59875200.0 ,         348     bi[1][4] = -2048058893.0/59875200.0 ,
349     bi[1][3] = -87098480009.0/5254048800.0 ,      349     bi[1][3] = -87098480009.0/5254048800.0 ,
350     bi[1][2] = -11513270273.0/3502699200.0 ,      350     bi[1][2] = -11513270273.0/3502699200.0 ,
351     //                                            351     //
352     bi[3][6] = -33197340367.0/1218433216.0 ,      352     bi[3][6] = -33197340367.0/1218433216.0 ,
353     bi[3][5] = -539868024987.0/6092166080.0 ,     353     bi[3][5] = -539868024987.0/6092166080.0 ,
354     bi[3][4] = -39991188681.0/374902528.0 ,       354     bi[3][4] = -39991188681.0/374902528.0 ,
355     bi[3][3] = -69509738227.0/1218433216.0 ,      355     bi[3][3] = -69509738227.0/1218433216.0 ,
356     bi[3][2] = -29327744613.0/2436866432.0 ,      356     bi[3][2] = -29327744613.0/2436866432.0 ,
357     //                                            357     //
358     bi[4][6] = -284800997201.0/19905339168.0 ,    358     bi[4][6] = -284800997201.0/19905339168.0 ,
359     bi[4][5] = -7896875450471.0/165877826400.0    359     bi[4][5] = -7896875450471.0/165877826400.0 ,
360     bi[4][4] = -333945812879.0/5671036800.0 ,     360     bi[4][4] = -333945812879.0/5671036800.0 ,
361     bi[4][3] = -16209923456237.0/497633479200.    361     bi[4][3] = -16209923456237.0/497633479200.0 ,
362     bi[4][2] = -2382590741699.0/331755652800.0    362     bi[4][2] = -2382590741699.0/331755652800.0 ,
363     //                                            363     //
364     bi[5][6] = -540919.0/741312.0 ,               364     bi[5][6] = -540919.0/741312.0 ,
365     bi[5][5] = -103626067.0/43243200.0 ,          365     bi[5][5] = -103626067.0/43243200.0 ,
366     bi[5][4] = -633779.0/211200.0 ,               366     bi[5][4] = -633779.0/211200.0 ,
367     bi[5][3] = -32406787.0/18532800.0 ,           367     bi[5][3] = -32406787.0/18532800.0 ,
368     bi[5][2] = -36591193.0/86486400.0 ,           368     bi[5][2] = -36591193.0/86486400.0 ,
369     //                                            369     //
370     bi[6][6] = 7157998304.0/374350977.0 ,         370     bi[6][6] = 7157998304.0/374350977.0 ,
371     bi[6][5] = 30405842464.0/623918295.0 ,        371     bi[6][5] = 30405842464.0/623918295.0 ,
372     bi[6][4] = 183022264.0/5332635.0 ,            372     bi[6][4] = 183022264.0/5332635.0 ,
373     bi[6][3] = -3357024032.0/1871754885.0 ,       373     bi[6][3] = -3357024032.0/1871754885.0 ,
374     bi[6][2] = -611586736.0/89131185.0 ,          374     bi[6][2] = -611586736.0/89131185.0 ,
375     //                                            375     //
376     bi[7][6] = -138073.0/9408.0 ,                 376     bi[7][6] = -138073.0/9408.0 ,
377     bi[7][5] = -719433.0/15680.0 ,                377     bi[7][5] = -719433.0/15680.0 ,
378     bi[7][4] = -1620541.0/31360.0 ,               378     bi[7][4] = -1620541.0/31360.0 ,
379     bi[7][3] = -385151.0/15680.0 ,                379     bi[7][3] = -385151.0/15680.0 ,
380     bi[7][2] = -65403.0/15680.0 ,                 380     bi[7][2] = -65403.0/15680.0 ,
381     //                                            381     //
382     bi[8][6] = 1245.0/64.0 ,                      382     bi[8][6] = 1245.0/64.0 ,
383     bi[8][5] = 3991.0/64.0 ,                      383     bi[8][5] = 3991.0/64.0 ,
384     bi[8][4] = 4715.0/64.0 ,                      384     bi[8][4] = 4715.0/64.0 ,
385     bi[8][3] = 2501.0/64.0 ,                      385     bi[8][3] = 2501.0/64.0 ,
386     bi[8][2] = 149.0/16.0 ,                       386     bi[8][2] = 149.0/16.0 ,
387     bi[8][1] = 1.0 ,                              387     bi[8][1] = 1.0 ,
388     //                                            388     //
389     bi[9][6] = 55.0/3.0 ,                         389     bi[9][6] = 55.0/3.0 ,
390     bi[9][5] = 71.0 ,                             390     bi[9][5] = 71.0 ,
391     bi[9][4] = 103.0 ,                            391     bi[9][4] = 103.0 ,
392     bi[9][3] = 199.0/3.0 ,                        392     bi[9][3] = 199.0/3.0 ,
393     bi[9][2] = 16.0 ,                             393     bi[9][2] = 16.0 ,
394     //                                            394     //
395     bi[10][6] = -1774004627.0/75810735.0 ,        395     bi[10][6] = -1774004627.0/75810735.0 ,
396     bi[10][5] = -1774004627.0/25270245.0 ,        396     bi[10][5] = -1774004627.0/25270245.0 ,
397     bi[10][4] = -26477681.0/359975.0 ,            397     bi[10][4] = -26477681.0/359975.0 ,
398     bi[10][3] = -11411880511.0/379053675.0 ,      398     bi[10][3] = -11411880511.0/379053675.0 ,
399     bi[10][2] = -423642896.0/126351225.0 ,        399     bi[10][2] = -423642896.0/126351225.0 ,
400     //                                            400     //
401     bi[11][6] = 35.0 ,                            401     bi[11][6] = 35.0 ,
402     bi[11][5] = 105.0 ,                           402     bi[11][5] = 105.0 ,
403     bi[11][4] = 117.0 ,                           403     bi[11][4] = 117.0 ,
404     bi[11][3] = 59.0 ,                            404     bi[11][3] = 59.0 ,
405     bi[11][2] = 12.0 ;                            405     bi[11][2] = 12.0 ;
406 }                                                 406 }
407                                                   407 
408 // -------------------------------------------    408 // ---------------------------------------------------------------------------------------
409                                                   409 
410 void G4FSALBogackiShampine45::interpolate( con    410 void G4FSALBogackiShampine45::interpolate( const G4double yInput[],
411                                            con    411                                            const G4double dydx[],
412                                                   412                                                  G4double yOut[],
413                                                   413                                                  G4double Step,
414                                                   414                                                  G4double tau )
415 {                                                 415 {
416     const G4double a91 = 455.0/6144.0 ,           416     const G4double a91 = 455.0/6144.0 ,
417                    a92 = 0.0 ,                    417                    a92 = 0.0 ,
418                    a93 = 10256301.0/35409920.0    418                    a93 = 10256301.0/35409920.0 ,
419                    a94 = 2307361.0/17971200.0     419                    a94 = 2307361.0/17971200.0 ,
420                    a95 = -387.0/102400.0 ,        420                    a95 = -387.0/102400.0 ,
421                    a96 = 73.0/5130.0 ,            421                    a96 = 73.0/5130.0 ,
422                    a97 = -7267.0/215040.0 ,       422                    a97 = -7267.0/215040.0 ,
423                    a98 = 1.0/32.0 ,               423                    a98 = 1.0/32.0 ,
424                                                   424 
425                    a101 = -837888343715.0/1317    425                    a101 = -837888343715.0/13176988637184.0 ,
426                    a102 = 30409415.0/52955362.    426                    a102 = 30409415.0/52955362.0 ,
427                    a103 = -48321525963.0/75916    427                    a103 = -48321525963.0/759168069632.0 ,
428                    a104 = 8530738453321.0/1976    428                    a104 = 8530738453321.0/197654829557760.0 ,
429                    a105 = 1361640523001.0/1626    429                    a105 = 1361640523001.0/1626788720640.0 ,
430                    a106 = -13143060689.0/38604    430                    a106 = -13143060689.0/38604458898.0 ,
431                    a107 = 18700221969.0/379584    431                    a107 = 18700221969.0/379584034816.0 ,
432                    a108 = -5831595.0/847285792    432                    a108 = -5831595.0/847285792.0 ,
433                    a109 = -5183640.0/26477681.    433                    a109 = -5183640.0/26477681.0 ,
434                                                   434 
435                    a111 = 98719073263.0/155196    435                    a111 = 98719073263.0/1551965184000.0 ,
436                    a112 = 1307.0/123552.0 ,       436                    a112 = 1307.0/123552.0 ,
437                    a113 = 4632066559387.0/7018    437                    a113 = 4632066559387.0/70181753241600.0 ,
438                    a114 = 7828594302389.0/3821    438                    a114 = 7828594302389.0/382182512025600.0 ,
439                    a115 = 40763687.0/110702592    439                    a115 = 40763687.0/11070259200.0 ,
440                    a116 = 34872732407.0/224610    440                    a116 = 34872732407.0/224610586200.0 ,
441                    a117 = -2561897.0/30105600.    441                    a117 = -2561897.0/30105600.0 ,
442                    a118 =  1.0/10.0 ,             442                    a118 =  1.0/10.0 ,
443                    a119 = -1.0/10.0 ,             443                    a119 = -1.0/10.0 ,
444                    a1110 = -1403317093.0/11371    444                    a1110 = -1403317093.0/11371610250.0 ;
445                                                   445 
446     const G4int numberOfVariables = GetNumberO    446     const G4int numberOfVariables = GetNumberOfVariables();
447                                                   447 
448     // Saving yInput because yInput and yOut c    448     // Saving yInput because yInput and yOut can be aliases for same array
449     //                                            449     //
450     for(G4int i=0; i<numberOfVariables; ++i)      450     for(G4int i=0; i<numberOfVariables; ++i)
451     {                                             451     {
452         yIn[i]=yInput[i];                         452         yIn[i]=yInput[i];
453     }                                             453     }
454                                                   454     
455     // The number of variables to be integrate    455     // The number of variables to be integrated over
456     //                                            456     //
457     yOut[7] = yTemp[7]  = yIn[7];                 457     yOut[7] = yTemp[7]  = yIn[7];
458                                                   458 
459     // Calculating extra stages                   459     // Calculating extra stages
460     //                                            460     //
461     for(G4int i=0; i<numberOfVariables; ++i)      461     for(G4int i=0; i<numberOfVariables; ++i)
462     {                                             462     {
463         yTemp[i] = yIn[i] + Step*(a91*dydx[i]     463         yTemp[i] = yIn[i] + Step*(a91*dydx[i] + a92*ak2[i] + a93*ak3[i] +
464                                   a94*ak4[i] +    464                                   a94*ak4[i] + a95*ak5[i] + a96*ak6[i] +
465                                   a97*ak7[i] +    465                                   a97*ak7[i] + a98*ak8[i] );
466     }                                             466     }
467                                                   467     
468     RightHandSide(yTemp, ak9);                    468     RightHandSide(yTemp, ak9);
469                                                   469     
470     for(G4int i=0; i<numberOfVariables; ++i)      470     for(G4int i=0; i<numberOfVariables; ++i)
471     {                                             471     {
472         yTemp[i] = yIn[i] + Step*(a101*dydx[i]    472         yTemp[i] = yIn[i] + Step*(a101*dydx[i] + a102*ak2[i] + a103*ak3[i] +
473                                   a104*ak4[i]     473                                   a104*ak4[i] + a105*ak5[i] + a106*ak6[i] +
474                                   a107*ak7[i]     474                                   a107*ak7[i] + a108*ak8[i] + a109*ak9[i] );
475     }                                             475     }
476                                                   476     
477     RightHandSide(yTemp, ak10);                   477     RightHandSide(yTemp, ak10);
478                                                   478 
479     for(G4int i=0; i<numberOfVariables; ++i)      479     for(G4int i=0; i<numberOfVariables; ++i)
480     {                                             480     {
481         yTemp[i] = yIn[i] + Step*(a111*dydx[i]    481         yTemp[i] = yIn[i] + Step*(a111*dydx[i] + a112*ak2[i] + a113*ak3[i] +
482                                   a114*ak4[i]     482                                   a114*ak4[i] + a115*ak5[i] + a116*ak6[i] +
483                                   a117*ak7[i]     483                                   a117*ak7[i] + a118*ak8[i] + a119*ak9[i] +
484                                   a1110*ak10[i    484                                   a1110*ak10[i] );
485     }                                             485     }
486                                                   486     
487     RightHandSide(yTemp, ak11);                   487     RightHandSide(yTemp, ak11);
488                                                   488     
489     G4double tau0 = tau;                          489     G4double tau0 = tau;
490                                                   490 
491     // Calculating the polynomials                491     // Calculating the polynomials
492     //                                            492     //
493     for(auto i=1; i<=11; ++i) // i is NOT the     493     for(auto i=1; i<=11; ++i) // i is NOT the coordinate no., it's stage no.
494     {                                             494     {
495         b[i] = 0.0;                               495         b[i] = 0.0;
496         tau = tau0;                               496         tau = tau0;
497         for(auto j=1; j<=6; ++j)                  497         for(auto j=1; j<=6; ++j)
498         {                                         498         {
499             b[i] += bi[i][j]*tau;                 499             b[i] += bi[i][j]*tau;
500             tau*=tau0;                            500             tau*=tau0;
501         }                                         501         }
502     }                                             502     }
503                                                   503     
504     for(G4int i=0; i<numberOfVariables; ++i)      504     for(G4int i=0; i<numberOfVariables; ++i)
505     {                                             505     {
506         yOut[i] = yIn[i] + Step*(b[1]*dydx[i]     506         yOut[i] = yIn[i] + Step*(b[1]*dydx[i] + b[2]*ak2[i] + b[3]*ak3[i] +
507                                  b[4]*ak4[i] +    507                                  b[4]*ak4[i] + b[5]*ak5[i] + b[6]*ak6[i] +
508                                  b[7]*ak7[i] +    508                                  b[7]*ak7[i] + b[8]*ak8[i] + b[9]*ak9[i] +
509                                  b[10]*ak10[i]    509                                  b[10]*ak10[i] + b[11]*ak11[i] );
510     }                                             510     }  
511 }                                                 511 }
512                                                   512