Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/magneticfield/src/G4FSALDormandPrince745.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/G4FSALDormandPrince745.cc (Version 11.3.0) and /geometry/magneticfield/src/G4FSALDormandPrince745.cc (Version 10.6.p1)


  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 // G4FSALDormandPrince745 implementation           26 // G4FSALDormandPrince745 implementation
 27 //                                                 27 //
 28 // The Butcher table of the FDormand-Prince-7-     28 // The Butcher table of the FDormand-Prince-7-4-5 method is as follows:
 29 //                                                 29 //
 30 //  0   |                                          30 //  0   |
 31 //  1/5 | 1/5                                      31 //  1/5 | 1/5
 32 //  3/10| 3/40        9/40                         32 //  3/10| 3/40        9/40
 33 //  4/5 | 44/45       56/15      32/9              33 //  4/5 | 44/45       56/15      32/9
 34 //  8/9 | 19372/6561  25360/2187 64448/6561  2     34 //  8/9 | 19372/6561  25360/2187 64448/6561  212/729
 35 //  1   | 9017/3168   355/33     46732/5247  4     35 //  1   | 9017/3168   355/33     46732/5247  49/176  5103/18656
 36 //  1   | 35/384      0          500/1113    1     36 //  1   | 35/384      0          500/1113    125/192 2187/6784    11/84
 37 //  ------------------------------------------     37 //  ---------------------------------------------------------------------------
 38 //        35/384      0          500/1113    1     38 //        35/384      0          500/1113    125/192 2187/6784    11/84    0
 39 //        5179/57600  0          7571/16695  3     39 //        5179/57600  0          7571/16695  393/640 92097/339200 187/2100 1/40
 40 //                                                 40 //
 41 // Created: Somnath Banerjee, Google Summer of     41 // Created: Somnath Banerjee, Google Summer of Code 2015, 25 May 2015
 42 // Supervision: John Apostolakis, CERN             42 // Supervision: John Apostolakis, CERN
 43 // -------------------------------------------     43 // --------------------------------------------------------------------
 44                                                    44 
 45 #include "G4FSALDormandPrince745.hh"               45 #include "G4FSALDormandPrince745.hh"
 46 #include "G4LineSection.hh"                        46 #include "G4LineSection.hh"
 47 #include <cmath>                                   47 #include <cmath>
 48                                                    48 
 49 // Constructor                                     49 // Constructor
 50 //                                                 50 //
 51 G4FSALDormandPrince745::G4FSALDormandPrince745     51 G4FSALDormandPrince745::G4FSALDormandPrince745(G4EquationOfMotion* EqRhs,
 52                                                    52                                                G4int noIntegrationVariables,
 53                                                    53                                                G4bool primary)
 54   : G4VFSALIntegrationStepper(EqRhs, noIntegra     54   : G4VFSALIntegrationStepper(EqRhs, noIntegrationVariables)
 55 {                                                  55 {
 56   const G4int numberOfVariables = noIntegratio     56   const G4int numberOfVariables = noIntegrationVariables;
 57                                                    57     
 58   // New Chunk of memory being created for use     58   // New Chunk of memory being created for use by the stepper
 59                                                    59     
 60   // aki - for storing intermediate RHS            60   // aki - for storing intermediate RHS
 61   //                                               61   //
 62   ak2 = new G4double[numberOfVariables];           62   ak2 = new G4double[numberOfVariables];
 63   ak3 = new G4double[numberOfVariables];           63   ak3 = new G4double[numberOfVariables];
 64   ak4 = new G4double[numberOfVariables];           64   ak4 = new G4double[numberOfVariables];
 65   ak5 = new G4double[numberOfVariables];           65   ak5 = new G4double[numberOfVariables];
 66   ak6 = new G4double[numberOfVariables];           66   ak6 = new G4double[numberOfVariables];
 67   ak7 = new G4double[numberOfVariables];           67   ak7 = new G4double[numberOfVariables];
 68                                                    68 
 69   // Also always allocate arrays for interpola     69   // Also always allocate arrays for interpolation stages    
 70   //                                               70   //
 71   ak8 = new G4double[numberOfVariables];           71   ak8 = new G4double[numberOfVariables];
 72   ak9 = new G4double[numberOfVariables];           72   ak9 = new G4double[numberOfVariables];
 73                                                    73     
 74   yTemp = new G4double[numberOfVariables] ;        74   yTemp = new G4double[numberOfVariables] ;
 75   yIn = new G4double[numberOfVariables] ;          75   yIn = new G4double[numberOfVariables] ;
 76                                                    76     
 77   pseudoDydx_for_DistChord = new G4double[numb     77   pseudoDydx_for_DistChord = new G4double[numberOfVariables];
 78                                                    78 
 79   fInitialDyDx = new G4double[numberOfVariable     79   fInitialDyDx = new G4double[numberOfVariables];    
 80   fLastInitialVector = new G4double[numberOfVa     80   fLastInitialVector = new G4double[numberOfVariables] ;
 81   fLastFinalVector = new G4double[numberOfVari     81   fLastFinalVector = new G4double[numberOfVariables] ;
 82   fLastDyDx = new G4double[numberOfVariables];     82   fLastDyDx = new G4double[numberOfVariables];
 83                                                    83     
 84   fMidVector = new G4double[numberOfVariables]     84   fMidVector = new G4double[numberOfVariables];
 85   fMidError =  new G4double[numberOfVariables]     85   fMidError =  new G4double[numberOfVariables];
 86                                                    86     
 87   if( primary )                                    87   if( primary )
 88   {                                                88   {
 89     fAuxStepper = new G4FSALDormandPrince745(E     89     fAuxStepper = new G4FSALDormandPrince745(EqRhs,numberOfVariables,!primary);
 90   }                                                90   }
 91 }                                                  91 }
 92                                                    92 
 93 // Destructor                                      93 // Destructor
 94 //                                                 94 //
 95 G4FSALDormandPrince745::~G4FSALDormandPrince74     95 G4FSALDormandPrince745::~G4FSALDormandPrince745()
 96 {                                                  96 {
 97     // Clear all previously allocated memory f     97     // Clear all previously allocated memory for stepper and DistChord
 98                                                    98 
 99     delete [] ak2;  ak2 = nullptr;                 99     delete [] ak2;  ak2 = nullptr;
100     delete [] ak3;  ak3 = nullptr;                100     delete [] ak3;  ak3 = nullptr;
101     delete [] ak4;  ak4 = nullptr;                101     delete [] ak4;  ak4 = nullptr;
102     delete [] ak5;  ak5 = nullptr;                102     delete [] ak5;  ak5 = nullptr;
103     delete [] ak6;  ak6 = nullptr;                103     delete [] ak6;  ak6 = nullptr;
104     delete [] ak7;  ak7 = nullptr;                104     delete [] ak7;  ak7 = nullptr;
105     delete [] ak8;  ak8 = nullptr;                105     delete [] ak8;  ak8 = nullptr;
106     delete [] ak9;  ak9 = nullptr;                106     delete [] ak9;  ak9 = nullptr;
107                                                   107     
108     delete [] yTemp; yTemp = nullptr;             108     delete [] yTemp; yTemp = nullptr;
109     delete [] yIn;   yIn = nullptr;               109     delete [] yIn;   yIn = nullptr;
110                                                   110 
111     delete [] pseudoDydx_for_DistChord;  pseud    111     delete [] pseudoDydx_for_DistChord;  pseudoDydx_for_DistChord = nullptr;
112     delete [] fInitialDyDx;              fInit    112     delete [] fInitialDyDx;              fInitialDyDx = nullptr;
113                                                   113     
114     delete [] fLastInitialVector;    fLastInit    114     delete [] fLastInitialVector;    fLastInitialVector = nullptr;
115     delete [] fLastFinalVector;      fLastFina    115     delete [] fLastFinalVector;      fLastFinalVector  = nullptr;
116     delete [] fLastDyDx;             fLastDyDx    116     delete [] fLastDyDx;             fLastDyDx  = nullptr;
117     delete [] fMidVector;            fMidVecto    117     delete [] fMidVector;            fMidVector = nullptr;
118     delete [] fMidError;             fMidError    118     delete [] fMidError;             fMidError  = nullptr;
119                                                   119     
120     delete fAuxStepper; fAuxStepper = nullptr;    120     delete fAuxStepper; fAuxStepper = nullptr;
121 }                                                 121 }
122                                                   122 
123 // Stepper                                        123 // Stepper
124 //                                                124 //
125 // Passing in the value of yInput[],the first     125 // Passing in the value of yInput[],the first time dydx[] and Step length
126 // Giving back yOut and yErr arrays for output    126 // Giving back yOut and yErr arrays for output and error respectively
127 //                                                127 //
128 void G4FSALDormandPrince745::Stepper(const G4d    128 void G4FSALDormandPrince745::Stepper(const G4double yInput[],
129                                      const G4d    129                                      const G4double dydx[],
130                                            G4d    130                                            G4double Step,
131                                            G4d    131                                            G4double yOut[],
132                                            G4d    132                                            G4double yErr[],
133                                            G4d    133                                            G4double nextDydx[] )
134 {                                                 134 {
135     G4int i;                                      135     G4int i;
136                                                   136     
137     // The various constants defined on the ba    137     // The various constants defined on the basis of butcher tableu
138                                                   138 
139     const G4double b21 = 0.2 ,                    139     const G4double b21 = 0.2 ,
140                    b31 = 3.0/40.0, b32 = 9.0/4    140                    b31 = 3.0/40.0, b32 = 9.0/40.0 ,
141                                                   141     
142                    b41 = 44.0/45.0, b42 = -56.    142                    b41 = 44.0/45.0, b42 = -56.0/15.0, b43 = 32.0/9.0,
143                                                   143     
144                    b51 = 19372.0/6561.0, b52 =    144                    b51 = 19372.0/6561.0, b52 = -25360.0/2187.0,
145                    b53 = 64448.0/6561.0, b54 =    145                    b53 = 64448.0/6561.0, b54 = -212.0/729.0 ,
146                                                   146     
147                    b61 = 9017.0/3168.0 , b62 =    147                    b61 = 9017.0/3168.0 , b62 =   -355.0/33.0,
148                    b63 =  46732.0/5247.0    ,     148                    b63 =  46732.0/5247.0    , b64 = 49.0/176.0 ,
149                    b65 = -5103.0/18656.0 ,        149                    b65 = -5103.0/18656.0 ,
150                                                   150     
151                    b71 = 35.0/384.0, b72 = 0.,    151                    b71 = 35.0/384.0, b72 = 0.,
152                    b73 = 500.0/1113.0, b74 = 1    152                    b73 = 500.0/1113.0, b74 = 125.0/192.0,
153                    b75 = -2187.0/6784.0, b76 =    153                    b75 = -2187.0/6784.0, b76 = 11.0/84.0,
154                                                   154     
155              //    c1 = 35.0/384.0, c2 = .0,      155              //    c1 = 35.0/384.0, c2 = .0,
156              //    c3 = 500.0/1113.0, c4 = 125    156              //    c3 = 500.0/1113.0, c4 = 125.0/192.0,
157              //    c5 = -2187.0/6784.0, c6 = 1    157              //    c5 = -2187.0/6784.0, c6 = 11.0/84.0,
158              //    c7 = 0,                        158              //    c7 = 0,
159                                                   159     
160                    dc1 = b71 - 5179.0/57600.0,    160                    dc1 = b71 - 5179.0/57600.0,
161                    dc2 = b72 - .0,                161                    dc2 = b72 - .0,
162                    dc3 = b73 - 7571.0/16695.0,    162                    dc3 = b73 - 7571.0/16695.0,
163                    dc4 = b74 - 393.0/640.0,       163                    dc4 = b74 - 393.0/640.0,
164                    dc5 = b75 + 92097.0/339200.    164                    dc5 = b75 + 92097.0/339200.0,
165                    dc6 = b76 - 187.0/2100.0,      165                    dc6 = b76 - 187.0/2100.0,
166                    dc7 = - 1.0/40.0 ;  //end o    166                    dc7 = - 1.0/40.0 ;  //end of declaration
167                                                   167 
168     const G4int numberOfVariables = GetNumberO    168     const G4int numberOfVariables = GetNumberOfVariables();
169       // The number of variables to be integra    169       // The number of variables to be integrated over
170                                                   170 
171     // Saving yInput because yInput and yOut c    171     // Saving yInput because yInput and yOut can be aliases for same array
172     //                                            172     //
173     for(i=0; i<numberOfVariables; ++i)            173     for(i=0; i<numberOfVariables; ++i)
174     {                                             174     {
175         yIn[i]          = yInput[i];              175         yIn[i]          = yInput[i];
176         fInitialDyDx[i] = dydx[i];                176         fInitialDyDx[i] = dydx[i];
177     }                                             177     }
178     // Ensure that time is initialised - in ca    178     // Ensure that time is initialised - in case it is not integrated
179     //                                            179     //
180     yOut[7] = yTemp[7]  = yInput[7];              180     yOut[7] = yTemp[7]  = yInput[7];
181     // RightHandSide(yIn, DyDx) ;   // 1st Ste    181     // RightHandSide(yIn, DyDx) ;   // 1st Step - Not doing, getting passed
182                                                   182     
183     for(i=0; i<numberOfVariables; ++i)            183     for(i=0; i<numberOfVariables; ++i)
184     {                                             184     {
185         yTemp[i] = yIn[i] + b21*Step*fInitialD    185         yTemp[i] = yIn[i] + b21*Step*fInitialDyDx[i] ;
186     }                                             186     }
187     RightHandSide(yTemp, ak2) ;              /    187     RightHandSide(yTemp, ak2) ;              // 2nd Step
188                                                   188     
189     for(i=0; i<numberOfVariables; ++i)            189     for(i=0; i<numberOfVariables; ++i)
190     {                                             190     {
191         yTemp[i] = yIn[i] + Step*(b31*fInitial    191         yTemp[i] = yIn[i] + Step*(b31*fInitialDyDx[i] + b32*ak2[i]) ;
192     }                                             192     }
193     RightHandSide(yTemp, ak3) ;              /    193     RightHandSide(yTemp, ak3) ;              // 3rd Step
194                                                   194     
195     for(i=0; i<numberOfVariables; ++i)            195     for(i=0; i<numberOfVariables; ++i)
196     {                                             196     {
197         yTemp[i] = yIn[i] + Step*(b41*fInitial    197         yTemp[i] = yIn[i] + Step*(b41*fInitialDyDx[i]
198                                 + b42*ak2[i] +    198                                 + b42*ak2[i] + b43*ak3[i]) ;
199     }                                             199     }
200     RightHandSide(yTemp, ak4) ;              /    200     RightHandSide(yTemp, ak4) ;              // 4th Step
201                                                   201     
202     for(i=0; i<numberOfVariables; ++i)            202     for(i=0; i<numberOfVariables; ++i)
203     {                                             203     {
204         yTemp[i] = yIn[i] + Step*(b51*fInitial    204         yTemp[i] = yIn[i] + Step*(b51*fInitialDyDx[i]
205                                 + b52*ak2[i] +    205                                 + b52*ak2[i] + b53*ak3[i] + b54*ak4[i]) ;
206     }                                             206     }
207     RightHandSide(yTemp, ak5) ;              /    207     RightHandSide(yTemp, ak5) ;              // 5th Step
208                                                   208     
209     for(i=0; i<numberOfVariables; ++i)            209     for(i=0; i<numberOfVariables; ++i)
210     {                                             210     {
211         yTemp[i] = yIn[i] + Step*(b61*fInitial    211         yTemp[i] = yIn[i] + Step*(b61*fInitialDyDx[i] + b62*ak2[i]
212                                 + b63*ak3[i] +    212                                 + b63*ak3[i] + b64*ak4[i] + b65*ak5[i]) ;
213     }                                             213     }
214     RightHandSide(yTemp, ak6) ;              /    214     RightHandSide(yTemp, ak6) ;              // 6th Step
215                                                   215     
216     for(i=0; i<numberOfVariables; ++i)            216     for(i=0; i<numberOfVariables; ++i)
217     {                                             217     {
218         yOut[i] = yIn[i] + Step*(b71*fInitialD    218         yOut[i] = yIn[i] + Step*(b71*fInitialDyDx[i] + b72*ak2[i] + b73*ak3[i]
219                                + b74*ak4[i] +     219                                + b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
220     }                                             220     }
221     RightHandSide(yOut, ak7);               //    221     RightHandSide(yOut, ak7);               //7th and Final step
222                                                   222     
223     for(i=0; i<numberOfVariables; ++i)            223     for(i=0; i<numberOfVariables; ++i)
224     {                                             224     {
225                                                   225         
226         yErr[i] = Step*(dc1*fInitialDyDx[i] +     226         yErr[i] = Step*(dc1*fInitialDyDx[i] + dc2*ak2[i] + dc3*ak3[i]
227                       + dc4*ak4[i] + dc5*ak5[i    227                       + dc4*ak4[i] + dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] ) ;
228                                                   228 
229         // Store Input and Final values, for p    229         // Store Input and Final values, for possible use in calculating chord
230         //                                        230         //
231         fLastInitialVector[i] = yIn[i] ;          231         fLastInitialVector[i] = yIn[i] ;
232         fLastFinalVector[i]   = yOut[i];          232         fLastFinalVector[i]   = yOut[i];
233         fLastDyDx[i]          = fInitialDyDx[i    233         fLastDyDx[i]          = fInitialDyDx[i];
234         nextDydx[i] = ak7[i];                     234         nextDydx[i] = ak7[i];
235     }                                             235     }
236     fLastStepLength = Step;                       236     fLastStepLength = Step;
237                                                   237     
238     return ;                                      238     return ;
239 }                                                 239 }
240                                                   240 
241 // DistChord                                      241 // DistChord
242 //                                                242 //
243 G4double G4FSALDormandPrince745::DistChord() c    243 G4double G4FSALDormandPrince745::DistChord() const
244 {                                                 244 {
245     G4double distLine, distChord;                 245     G4double distLine, distChord;
246     G4ThreeVector initialPoint, finalPoint, mi    246     G4ThreeVector initialPoint, finalPoint, midPoint;
247                                                   247     
248     // Store last initial and final points        248     // Store last initial and final points
249     // (they will be overwritten in self-Stepp    249     // (they will be overwritten in self-Stepper call!)
250     //                                            250     //
251     initialPoint = G4ThreeVector(fLastInitialV    251     initialPoint = G4ThreeVector(fLastInitialVector[0],
252                                  fLastInitialV    252                                  fLastInitialVector[1], fLastInitialVector[2]);
253     finalPoint   = G4ThreeVector(fLastFinalVec    253     finalPoint   = G4ThreeVector(fLastFinalVector[0],
254                                  fLastFinalVec    254                                  fLastFinalVector[1],  fLastFinalVector[2]);
255                                                   255     
256     // Do half a step using StepNoErr             256     // Do half a step using StepNoErr
257                                                   257     
258     fAuxStepper->Stepper( fLastInitialVector,     258     fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
259                           fMidVector, fMidErro    259                           fMidVector, fMidError, pseudoDydx_for_DistChord );
260                                                   260     
261     midPoint = G4ThreeVector( fMidVector[0], f    261     midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2] );
262                                                   262     
263     // Use stored values of Initial and Endpoi    263     // Use stored values of Initial and Endpoint + new Midpoint to evaluate
264     // distance of Chord                          264     // distance of Chord
265     //                                            265     //
266     if (initialPoint != finalPoint)               266     if (initialPoint != finalPoint)
267     {                                             267     {
268         distLine  = G4LineSection::Distline( m    268         distLine  = G4LineSection::Distline( midPoint,initialPoint,finalPoint );
269         distChord = distLine;                     269         distChord = distLine;
270     }                                             270     }
271     else                                          271     else
272     {                                             272     {
273         distChord = (midPoint-initialPoint).ma    273         distChord = (midPoint-initialPoint).mag();
274     }                                             274     }
275     return distChord;                             275     return distChord;
276 }                                                 276 }
277                                                   277 
278 // interpolate                                    278 // interpolate
279 //                                                279 //
280 void G4FSALDormandPrince745::interpolate( cons    280 void G4FSALDormandPrince745::interpolate( const G4double yInput[],
281                                           cons    281                                           const G4double dydx[],
282                                                   282                                                 G4double yOut[],
283                                                   283                                                 G4double Step,
284                                                   284                                                 G4double tau)
285 {                                                 285 {    
286     G4double bf1, bf2, bf3, bf4, bf5, bf6, bf7    286     G4double bf1, bf2, bf3, bf4, bf5, bf6, bf7;
287                                                   287 
288     const G4int numberOfVariables = GetNumberO    288     const G4int numberOfVariables = GetNumberOfVariables();
289                                                   289     
290     G4double tau0 = tau;                          290     G4double tau0 = tau;
291                                                   291     
292     for(G4int i=0;i<numberOfVariables; ++i)       292     for(G4int i=0;i<numberOfVariables; ++i)
293     {                                             293     {
294         yIn[i]=yInput[i];                         294         yIn[i]=yInput[i];
295     }                                             295     }
296                                                   296     
297     G4double tau_2 = tau0*tau0 ,                  297     G4double tau_2 = tau0*tau0 ,
298              tau_3 = tau0*tau_2,                  298              tau_3 = tau0*tau_2,
299              tau_4 = tau_2*tau_2;                 299              tau_4 = tau_2*tau_2;
300                                                   300     
301     bf1 = (157015080.0*tau_4 - 13107642775.0*t    301     bf1 = (157015080.0*tau_4 - 13107642775.0*tau_3
302          + 34969693132.0*tau_2- 32272833064.0*    302          + 34969693132.0*tau_2- 32272833064.0*tau + 11282082432.0)
303          / 11282082432.0;                         303          / 11282082432.0;
304     bf2 = 0.0;                                    304     bf2 = 0.0;
305     bf3 = - 100.0*tau*(15701508.0*tau_3 - 9141    305     bf3 = - 100.0*tau*(15701508.0*tau_3 - 914128567.0*tau_2
306                      + 2074956840.0*tau - 1323    306                      + 2074956840.0*tau - 1323431896.0) / 32700410799.0;
307     bf4 = 25.0*tau*(94209048.0*tau_3- 15184142    307     bf4 = 25.0*tau*(94209048.0*tau_3- 1518414297.0*tau_2
308                   + 2460397220.0*tau - 8892898    308                   + 2460397220.0*tau - 889289856.0)
309                   / 5641041216.0;                 309                   / 5641041216.0;
310     bf5 = -2187.0*tau*(52338360.0*tau_3 - 4518    310     bf5 = -2187.0*tau*(52338360.0*tau_3 - 451824525.0*tau_2
311                      + 687873124.0*tau - 25900    311                      + 687873124.0*tau - 259006536.0)
312                      / 199316789632.0;            312                      / 199316789632.0;
313     bf6 =  11.0*tau*(106151040.0*tau_3- 661884    313     bf6 =  11.0*tau*(106151040.0*tau_3- 661884105.0*tau_2
314                    + 946554244.0*tau - 3614407    314                    + 946554244.0*tau - 361440756.0)
315                    / 2467955532.0;                315                    / 2467955532.0;
316     bf7 = tau*(1.0 - tau)*(8293050.0*tau_2 - 8    316     bf7 = tau*(1.0 - tau)*(8293050.0*tau_2 - 82437520.0*tau + 44764047.0)
317                          / 29380423.0;            317                          / 29380423.0;
318                                                   318 
319     for(G4int i=0; i<numberOfVariables; ++i)      319     for(G4int i=0; i<numberOfVariables; ++i)
320     {                                             320     {
321       yOut[i] = yIn[i] + Step*tau*(bf1*dydx[i]    321       yOut[i] = yIn[i] + Step*tau*(bf1*dydx[i] + bf2*ak2[i] + bf3*ak3[i]
322                                  + bf4*ak4[i]     322                                  + bf4*ak4[i] + bf5*ak5[i] + bf6*ak6[i]
323                                  + bf7*ak7[i]     323                                  + bf7*ak7[i] );
324     }                                             324     }
325 }                                                 325 }
326                                                   326 
327 // SetupInterpolate                               327 // SetupInterpolate
328 //                                                328 //
329 void G4FSALDormandPrince745::SetupInterpolate(    329 void G4FSALDormandPrince745::SetupInterpolate(const G4double yInput[],
330                                                   330                                               const G4double dydx[],
331                                                   331                                               const G4double Step )
332 {                                                 332 {    
333     // Coefficients for the additional stages     333     // Coefficients for the additional stages
334     //                                            334     //
335     G4double b81 =  6245.0/62208.0 ,              335     G4double b81 =  6245.0/62208.0 ,
336              b82 =  0.0 ,                         336              b82 =  0.0 ,
337              b83 =  8875.0/103032.0 ,             337              b83 =  8875.0/103032.0 ,
338              b84 = -125.0/1728.0 ,                338              b84 = -125.0/1728.0 ,
339              b85 =  801.0/13568.0 ,               339              b85 =  801.0/13568.0 ,
340              b86 = -13519.0/368064.0 ,            340              b86 = -13519.0/368064.0 ,
341              b87 =  11105.0/368064.0 ,            341              b87 =  11105.0/368064.0 ,
342                                                   342     
343              b91 =  632855.0/4478976.0 ,          343              b91 =  632855.0/4478976.0 ,
344              b92 =  0.0 ,                         344              b92 =  0.0 ,
345              b93 =  4146875.0/6491016.0 ,         345              b93 =  4146875.0/6491016.0 ,
346              b94 =  5490625.0/14183424.0 ,        346              b94 =  5490625.0/14183424.0 ,
347              b95 = -15975.0/108544.0 ,            347              b95 = -15975.0/108544.0 ,
348              b96 =  8295925.0/220286304.0 ,       348              b96 =  8295925.0/220286304.0 ,
349              b97 = -1779595.0/62938944.0 ,        349              b97 = -1779595.0/62938944.0 ,
350              b98 = -805.0/4104.0 ;                350              b98 = -805.0/4104.0 ;
351                                                   351     
352     const G4int numberOfVariables = GetNumberO    352     const G4int numberOfVariables = GetNumberOfVariables();
353                                                   353     
354     // Saving yInput because yInput and yOut c    354     // Saving yInput because yInput and yOut can be aliases for same array
355     //                                            355     //
356     for(G4int i=0; i<numberOfVariables; ++i)      356     for(G4int i=0; i<numberOfVariables; ++i)
357     {                                             357     {
358       yIn[i] = yInput[i];                         358       yIn[i] = yInput[i];
359     }                                             359     }
360                                                   360     
361     yTemp[7] = yIn[7];                            361     yTemp[7] = yIn[7];
362                                                   362     
363     // Evaluate the extra stages                  363     // Evaluate the extra stages
364     //                                            364     //
365     for(G4int i=0; i<numberOfVariables; ++i)      365     for(G4int i=0; i<numberOfVariables; ++i)
366     {                                             366     {
367         yTemp[i] = yIn[i] + Step*( b81*dydx[i]    367         yTemp[i] = yIn[i] + Step*( b81*dydx[i] + b82*ak2[i] + b83*ak3[i] +
368                                    b84*ak4[i]     368                                    b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
369                                    b87*ak7[i]     369                                    b87*ak7[i] );
370     }                                             370     }
371     RightHandSide( yTemp, ak8 );                  371     RightHandSide( yTemp, ak8 );              // 8th Stage
372                                                   372     
373     for(G4int i=0; i<numberOfVariables; ++i)      373     for(G4int i=0; i<numberOfVariables; ++i)
374     {                                             374     {
375         yTemp[i] = yIn[i] + Step * ( b91*dydx[    375         yTemp[i] = yIn[i] + Step * ( b91*dydx[i] + b92*ak2[i] + b93*ak3[i] +
376                                    b94*ak4[i]     376                                    b94*ak4[i] + b95*ak5[i] + b96*ak6[i] +
377                                    b97*ak7[i]     377                                    b97*ak7[i] + b98*ak8[i] );
378     }                                             378     }
379     RightHandSide( yTemp, ak9 );                  379     RightHandSide( yTemp, ak9 );              // 9th Stage
380 }                                                 380 }
381                                                   381 
382 // Interpolate                                    382 // Interpolate
383 //                                                383 //
384 void G4FSALDormandPrince745::Interpolate( cons    384 void G4FSALDormandPrince745::Interpolate( const G4double yInput[],
385                                           cons    385                                           const G4double dydx[],
386                                           cons    386                                           const G4double Step,
387                                                   387                                                 G4double yOut[],
388                                                   388                                                 G4double tau )
389 {                                                 389 {
390     // Define the coefficients for the polynom    390     // Define the coefficients for the polynomials
391                                                   391 
392     G4double bi[10][5], b[10];                    392     G4double bi[10][5], b[10];
393     G4int numberOfVariables = GetNumberOfVaria    393     G4int numberOfVariables = GetNumberOfVariables();
394                                                   394     
395     //  COEFFICIENTS OF   bi[1]                   395     //  COEFFICIENTS OF   bi[1]
396     bi[1][0] =  1.0 ,                             396     bi[1][0] =  1.0 ,
397     bi[1][1] = -38039.0/7040.0 ,                  397     bi[1][1] = -38039.0/7040.0 ,
398     bi[1][2] =  125923.0/10560.0 ,                398     bi[1][2] =  125923.0/10560.0 ,
399     bi[1][3] = -19683.0/1760.0 ,                  399     bi[1][3] = -19683.0/1760.0 ,
400     bi[1][4] =  3303.0/880.0 ,                    400     bi[1][4] =  3303.0/880.0 ,
401     //  --------------------------------------    401     //  --------------------------------------------------------
402     //                                            402     //
403     //  COEFFICIENTS OF  bi[2]                    403     //  COEFFICIENTS OF  bi[2]
404     bi[2][0] =  0.0 ,                             404     bi[2][0] =  0.0 ,
405     bi[2][1] =  0.0 ,                             405     bi[2][1] =  0.0 ,
406     bi[2][2] =  0.0 ,                             406     bi[2][2] =  0.0 ,
407     bi[2][3] =  0.0 ,                             407     bi[2][3] =  0.0 ,
408     bi[2][4] =  0.0 ,                             408     bi[2][4] =  0.0 ,
409     //  --------------------------------------    409     //  --------------------------------------------------------
410     //                                            410     //
411     //  COEFFICIENTS OF  bi[3]                    411     //  COEFFICIENTS OF  bi[3]
412     bi[3][0] =  0.0 ,                             412     bi[3][0] =  0.0 ,
413     bi[3][1] = -12500.0/4081.0 ,                  413     bi[3][1] = -12500.0/4081.0 ,
414     bi[3][2] =  205000.0/12243.0 ,                414     bi[3][2] =  205000.0/12243.0 ,
415     bi[3][3] = -90000.0/4081.0 ,                  415     bi[3][3] = -90000.0/4081.0 ,
416     bi[3][4] =  36000.0/4081.0 ,                  416     bi[3][4] =  36000.0/4081.0 ,
417     //  --------------------------------------    417     //  --------------------------------------------------------
418     //                                            418     //
419     //  COEFFICIENTS OF  bi[4]                    419     //  COEFFICIENTS OF  bi[4]
420     bi[4][0] =  0.0 ,                             420     bi[4][0] =  0.0 ,
421     bi[4][1] = -3125.0/704.0 ,                    421     bi[4][1] = -3125.0/704.0 ,
422     bi[4][2] =  25625.0/1056.0 ,                  422     bi[4][2] =  25625.0/1056.0 ,
423     bi[4][3] = -5625.0/176.0 ,                    423     bi[4][3] = -5625.0/176.0 ,
424     bi[4][4] =  1125.0/88.0 ,                     424     bi[4][4] =  1125.0/88.0 ,
425     //  --------------------------------------    425     //  --------------------------------------------------------
426     //                                            426     //
427     //  COEFFICIENTS OF  bi[5]                    427     //  COEFFICIENTS OF  bi[5]
428     bi[5][0] =  0.0 ,                             428     bi[5][0] =  0.0 ,
429     bi[5][1] =  164025.0/74624.0 ,                429     bi[5][1] =  164025.0/74624.0 ,
430     bi[5][2] = -448335.0/37312.0 ,                430     bi[5][2] = -448335.0/37312.0 ,
431     bi[5][3] =  295245.0/18656.0 ,                431     bi[5][3] =  295245.0/18656.0 ,
432     bi[5][4] = -59049.0/9328.0 ,                  432     bi[5][4] = -59049.0/9328.0 ,
433     //  --------------------------------------    433     //  --------------------------------------------------------
434     //                                            434     //
435     //  COEFFICIENTS OF  bi[6]                    435     //  COEFFICIENTS OF  bi[6]
436     bi[6][0] =  0.0 ,                             436     bi[6][0] =  0.0 ,
437     bi[6][1] = -25.0/28.0 ,                       437     bi[6][1] = -25.0/28.0 ,
438     bi[6][2] =  205.0/42.0 ,                      438     bi[6][2] =  205.0/42.0 ,
439     bi[6][3] = -45.0/7.0 ,                        439     bi[6][3] = -45.0/7.0 ,
440     bi[6][4] =  18.0/7.0 ,                        440     bi[6][4] =  18.0/7.0 ,
441     //  --------------------------------------    441     //  --------------------------------------------------------
442     //                                            442     //
443     //  COEFFICIENTS OF  bi[7]                    443     //  COEFFICIENTS OF  bi[7]
444     bi[7][0] =  0.0 ,                             444     bi[7][0] =  0.0 ,
445     bi[7][1] = -2.0/11.0 ,                        445     bi[7][1] = -2.0/11.0 ,
446     bi[7][2] =  73.0/55.0 ,                       446     bi[7][2] =  73.0/55.0 ,
447     bi[7][3] = -171.0/55.0 ,                      447     bi[7][3] = -171.0/55.0 ,
448     bi[7][4] =  108.0/55.0 ,                      448     bi[7][4] =  108.0/55.0 ,
449     //  --------------------------------------    449     //  --------------------------------------------------------
450     //                                            450     //
451     //  COEFFICIENTS OF  bi[8]                    451     //  COEFFICIENTS OF  bi[8]
452     bi[8][0] =  0.0 ,                             452     bi[8][0] =  0.0 ,
453     bi[8][1] =  189.0/22.0 ,                      453     bi[8][1] =  189.0/22.0 ,
454     bi[8][2] = -1593.0/55.0 ,                     454     bi[8][2] = -1593.0/55.0 ,
455     bi[8][3] =  3537.0/110.0 ,                    455     bi[8][3] =  3537.0/110.0 ,
456     bi[8][4] = -648.0/55.0 ,                      456     bi[8][4] = -648.0/55.0 ,
457     //  --------------------------------------    457     //  --------------------------------------------------------
458     //                                            458     //
459     //  COEFFICIENTS OF  bi[9]                    459     //  COEFFICIENTS OF  bi[9]
460     bi[9][0] =  0.0 ,                             460     bi[9][0] =  0.0 ,
461     bi[9][1] =  351.0/110.0 ,                     461     bi[9][1] =  351.0/110.0 ,
462     bi[9][2] = -999.0/55.0 ,                      462     bi[9][2] = -999.0/55.0 ,
463     bi[9][3] =  2943.0/110.0 ,                    463     bi[9][3] =  2943.0/110.0 ,
464     bi[9][4] = -648.0/55.0 ;                      464     bi[9][4] = -648.0/55.0 ;
465     //  --------------------------------------    465     //  --------------------------------------------------------
466                                                   466 
467     for(G4int i = 0; i< numberOfVariables; ++i    467     for(G4int i = 0; i< numberOfVariables; ++i)
468     {                                             468     {
469         yIn[i] = yInput[i];                       469         yIn[i] = yInput[i];
470     }                                             470     }
471                                                   471     
472     G4double tau0 = tau;                          472     G4double tau0 = tau;
473                                                   473 
474     // Calculating the polynomials                474     // Calculating the polynomials
475     //                                            475     //    
476     for(auto i=1; i<=9; ++i)  // i is NOT the     476     for(auto i=1; i<=9; ++i)  // i is NOT the coordinate no., it's stage no.
477     {                                             477     {
478         b[i] = 0;                                 478         b[i] = 0;
479         tau = 1.0;                                479         tau = 1.0;
480         for(auto j=0; j<=4; ++j)                  480         for(auto j=0; j<=4; ++j)
481         {                                         481         {
482             b[i] += bi[i][j]*tau;                 482             b[i] += bi[i][j]*tau;
483             tau*=tau0;                            483             tau*=tau0;
484         }                                         484         }
485     }                                             485     }
486                                                   486     
487     for(G4int i=0; i<numberOfVariables; ++i) /    487     for(G4int i=0; i<numberOfVariables; ++i) // Here i IS the coordinate no.
488     {                                             488     {
489         yOut[i] = yIn[i] + Step*tau0*(b[1]*dyd    489         yOut[i] = yIn[i] + Step*tau0*(b[1]*dydx[i] + b[2]*ak2[i] + b[3]*ak3[i] +
490                                       b[4]*ak4    490                                       b[4]*ak4[i] + b[5]*ak5[i] + b[6]*ak6[i] +
491                                       b[7]*ak7    491                                       b[7]*ak7[i] + b[8]*ak8[i] + b[9]*ak9[i] );
492     }                                             492     }
493 }                                                 493 }
494                                                   494