Geant4 Cross Reference

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


  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 // G4TsitourasRK45 implementation              << 
 27 //                                                 26 //
 28 // Author: Somnath Banerjee, Google Summer of  <<  27 //  Tsitouras - 5(4) RK steppers   ( non-FSAL version )
 29 // Supervision: John Apostolakis, CERN         <<  28 //
                                                   >>  29 //  Implements RK tableau from 'Table 1' of 
                                                   >>  30 //    C. Tsitouras, “Runge–Kutta pairs of order 5(4) satisfying only
                                                   >>  31 //    the first column simplifying assumption,”
                                                   >>  32 //   Computers & Mathematics with Applications,
                                                   >>  33 //   vol. 62, no. 2, pp. 770–775, 2011.
                                                   >>  34 //
                                                   >>  35 //  Adaptation / Geant4 implementation by Somnath Banerjee
                                                   >>  36 //  Supervision / code review: John Apostolakis
                                                   >>  37 //
                                                   >>  38 // Sponsored by Google in Google Summer of Code 2015.
                                                   >>  39 //
                                                   >>  40 // First version: 12 June 2015
                                                   >>  41 //
 30 // -------------------------------------------     42 // -------------------------------------------------------------------
 31                                                    43 
 32 #include "G4TsitourasRK45.hh"                      44 #include "G4TsitourasRK45.hh"
 33 #include "G4LineSection.hh"                        45 #include "G4LineSection.hh"
 34                                                    46 
 35 //////////////////////////////////////////////     47 /////////////////////////////////////////////////////////////////////
 36 //                                                 48 //
 37 // Constructor                                     49 // Constructor
 38 //                                             <<  50 
 39 G4TsitourasRK45::G4TsitourasRK45(G4EquationOfM     51 G4TsitourasRK45::G4TsitourasRK45(G4EquationOfMotion *EqRhs, 
 40                                  G4int noInteg <<  52          G4int noIntegrationVariables, 
 41                                  G4bool primar <<  53          G4bool primary)
 42   : G4MagIntegratorStepper(EqRhs, noIntegratio <<  54   : G4MagIntegratorStepper(EqRhs, noIntegrationVariables),
                                                   >>  55     fLastStepLength(0.), fAuxStepper(0)
 43 {                                                  56 {
 44   const G4int numberOfVariables = noIntegratio     57   const G4int numberOfVariables = noIntegrationVariables;
                                                   >>  58   // G4cout << "G4TsitourasRK45 constructor called." << G4endl;
 45                                                    59   
 46   ak2 = new G4double[numberOfVariables] ;          60   ak2 = new G4double[numberOfVariables] ;  
 47   ak3 = new G4double[numberOfVariables] ;          61   ak3 = new G4double[numberOfVariables] ; 
 48   ak4 = new G4double[numberOfVariables] ;          62   ak4 = new G4double[numberOfVariables] ; 
 49   ak5 = new G4double[numberOfVariables] ;          63   ak5 = new G4double[numberOfVariables] ; 
 50   ak6 = new G4double[numberOfVariables] ;          64   ak6 = new G4double[numberOfVariables] ; 
 51   ak7 = new G4double[numberOfVariables] ;          65   ak7 = new G4double[numberOfVariables] ;
 52   ak8 = new G4double[numberOfVariables] ;          66   ak8 = new G4double[numberOfVariables] ;
 53                                                    67 
                                                   >>  68 
 54   // Must ensure space extra 'state' variables     69   // Must ensure space extra 'state' variables exists - i.e. yIn[7]
 55   //                                           << 
 56   const G4int numStateMax  = std::max(GetNumbe     70   const G4int numStateMax  = std::max(GetNumberOfStateVariables(), 8);  
 57   const G4int numStateVars = std::max(noIntegr     71   const G4int numStateVars = std::max(noIntegrationVariables,
 58                                       numState     72                                       numStateMax );
                                                   >>  73                                    // GetNumberOfStateVariables() ); 
 59                                                    74                                       
 60   yTemp = new G4double[numStateVars] ;             75   yTemp = new G4double[numStateVars] ;
 61   yIn = new G4double[numStateVars] ;               76   yIn = new G4double[numStateVars] ;
 62                                                    77 
 63   fLastInitialVector = new G4double[numberOfVa     78   fLastInitialVector = new G4double[numberOfVariables] ;
 64   fLastFinalVector = new G4double[numberOfVari     79   fLastFinalVector = new G4double[numberOfVariables] ;
 65                                                    80 
 66   fLastDyDx = new G4double[numberOfVariables];     81   fLastDyDx = new G4double[numberOfVariables];
 67                                                    82 
 68   fMidVector = new G4double[numberOfVariables]     83   fMidVector = new G4double[numberOfVariables];
 69   fMidError =  new G4double[numberOfVariables]     84   fMidError =  new G4double[numberOfVariables];
 70                                                << 
 71   if( primary )                                    85   if( primary )
 72   {                                                86   { 
 73     fAuxStepper = new G4TsitourasRK45(EqRhs, n     87     fAuxStepper = new G4TsitourasRK45(EqRhs, numberOfVariables, !primary);
 74   }                                                88   }
 75 }                                                  89 }
 76                                                    90 
 77 //////////////////////////////////////////////     91 /////////////////////////////////////////////////////////////////////
 78 //                                                 92 //
 79 // Destructor                                      93 // Destructor
 80 //                                             <<  94 
 81 G4TsitourasRK45::~G4TsitourasRK45()                95 G4TsitourasRK45::~G4TsitourasRK45()
 82 {                                                  96 {
 83   delete [] ak2;                               <<  97   delete[] ak2;
 84   delete [] ak3;                               <<  98   delete[] ak3;
 85   delete [] ak4;                               <<  99   delete[] ak4;
 86   delete [] ak5;                               << 100   delete[] ak5;
 87   delete [] ak6;                               << 101   delete[] ak6;
 88   delete [] ak7;                               << 102   delete[] ak7;
 89   delete [] ak8;                               << 103   delete[] ak8;
 90                                                << 104 
 91   delete [] yTemp;                             << 105   delete[] yTemp;
 92   delete [] yIn;                               << 106   delete[] yIn;
 93                                                << 107 
 94   delete [] fLastInitialVector;                << 108   delete[] fLastInitialVector;
 95   delete [] fLastFinalVector;                  << 109   delete[] fLastFinalVector;
 96   delete [] fLastDyDx;                         << 110   delete[] fLastDyDx;
 97   delete [] fMidVector;                        << 111   delete[] fMidVector;
 98   delete [] fMidError;                         << 112   delete[] fMidError; 
 99                                                   113 
100   delete fAuxStepper;                             114   delete fAuxStepper;
101 }                                                 115 }
102                                                   116 
103 // The following coefficients have been obtain << 117 //The following coefficients have been obtained from
104 // Table 1: The Coefficients of the new pair      118 // Table 1: The Coefficients of the new pair
105 //                                             << 119 //---Ref---
106 // C. Tsitouras, "Runge–Kutta pairs of order << 120 // C. Tsitouras, “Runge–Kutta pairs of order 5(4) satisfying only
107 //                the first column simplifying << 121 // the first column simplifying assumption,”
108 // Computers & Mathematics with Applications,  << 122 // Computers & Mathematics with Applications,
109 //                                             << 123 // vol. 62, no. 2, pp. 770–775, 2011.
110 // A corresponding matlab code was also found  << 124 //-----------------------------------
111 //   http://users.ntua.gr/tsitoura/new54.m     << 125 // A corresponding matlab code was also found @ http://users.ntua.gr/tsitoura/new54.m
112 //                                             << 126 
113 // Doing a step                                   127 // Doing a step
114 //                                             << 
115 void                                              128 void
116 G4TsitourasRK45::Stepper( const G4double yInpu << 129 G4TsitourasRK45::Stepper(  const G4double yInput[],
117                           const G4double dydx[ << 130                          const G4double dydx[],
118                                 G4double Step, << 131                                G4double Step,
119                                 G4double yOut[ << 132                                G4double yOut[],
120                                 G4double yErr[ << 133                                G4double yErr[])
121 {                                                 134 {
122     const G4double b21 = 0.161 ,               << 135     G4int i;
123                    b31 = -0.008480655492356988 << 136 
124                    b32 = 0.335480655492356989  << 137     const G4double
125                                                << 138     b21 = 0.161 ,
126                    b41 =  2.89715305710549343  << 139     
127                    b42 = -6.35944848997507484  << 140     b31 = -0.00848065549235698854 ,
128                    b43 = 4.36229543286958141 , << 141     b32 = 0.335480655492356989 ,
129                                                << 
130                    b51 = 5.325864828439257,    << 
131                    b52 = -11.748883564062828,  << 
132                    b53 = 7.49553934288983621 , << 
133                    b54 = -0.09249506636175525, << 
134                                                << 
135                    b61 = 5.8614554429464200,   << 
136                    b62 = -12.9209693178471093  << 
137                    b63 = 8.1593678985761586 ,  << 
138                    b64 = -0.071584973281400997 << 
139                    b65 = -0.028269050394068382 << 
140                                                << 
141                    b71 = 0.0964607668180652295 << 
142                    b72 = 0.01,                 << 
143                    b73 = 0.479889650414499575, << 
144                    b74 = 1.37900857410374189,  << 
145                    b75 = -3.2900695154360807,  << 
146                    b76 = 2.32471052409977398,  << 
147                                                << 
148              //    c1 = 0.001780011052226 ,    << 
149              //    c2 = 0.000816434459657 ,    << 
150              //    c3 = -0.007880878010262 ,   << 
151              //    c4 = 0.144711007173263 ,    << 
152              //    c5 = -0.582357165452555 ,   << 
153              //    c6 = 0.458082105929187 ,    << 
154              //    c7 = 1.0/66.0 ;             << 
155                                                << 
156                    dc1 = 0.0935237485818927066 << 
157                    dc2 = 0.0086528831415663676 << 
158                    dc3 = 0.492893099131431868  << 
159                    dc4 = 1.14023541226785810 - << 
160                    dc5 = - 2.3291801924393646  << 
161                    dc6 = 1.56887504931661552 - << 
162                    dc7 = 0.025; //- 1.0/66.0 ; << 
163                                                << 
164              //    dc1 = -3.0/1280.0,          << 
165              //    dc2 = 0.0,                  << 
166              //    dc3 = 6561.0/632320.0,      << 
167              //    dc4 = -343.0/20800.0,       << 
168              //    dc5 = 243.0/12800.0,        << 
169              //    dc6 = -1.0/95.0,            << 
170              //    dc7 = 0.0   ;               << 
171                                                   142     
172     const G4int numberOfVariables = GetNumberO << 143     b41 =  2.89715305710549343 ,
                                                   >> 144     b42 = -6.35944848997507484 ,
                                                   >> 145     b43 = 4.36229543286958141 ,
                                                   >> 146 
                                                   >> 147     b51 = 5.325864828439257,
                                                   >> 148     b52 = -11.748883564062828,
                                                   >> 149     b53 = 7.49553934288983621 ,
                                                   >> 150     b54 = -0.09249506636175525,
                                                   >> 151 
                                                   >> 152     b61 = 5.8614554429464200,
                                                   >> 153     b62 = -12.9209693178471093 ,
                                                   >> 154     b63 = 8.1593678985761586 ,
                                                   >> 155     b64 = -0.071584973281400997,
                                                   >> 156     b65 = -0.0282690503940683829,
                                                   >> 157 
                                                   >> 158     
                                                   >> 159     b71 = 0.0964607668180652295 ,
                                                   >> 160     b72 = 0.01, 
                                                   >> 161     b73 = 0.479889650414499575,
                                                   >> 162     b74 = 1.37900857410374189,
                                                   >> 163     b75 = -3.2900695154360807,
                                                   >> 164     b76 = 2.32471052409977398,
                                                   >> 165     
                                                   >> 166 //    c1 = 0.001780011052226 ,
                                                   >> 167 //    c2 = 0.000816434459657 ,
                                                   >> 168 //    c3 = -0.007880878010262 ,
                                                   >> 169 //    c4 = 0.144711007173263 ,
                                                   >> 170 //    c5 = -0.582357165452555 ,
                                                   >> 171 //    c6 = 0.458082105929187 ,
                                                   >> 172 //    c7 = 1.0/66.0 ;
                                                   >> 173 
                                                   >> 174     dc1 = 0.0935237485818927066 - b71 , // - 0.001780011052226,
                                                   >> 175     dc2 = 0.00865288314156636761 - b72, // - 0.000816434459657,
                                                   >> 176     dc3 = 0.492893099131431868 - b73 ,  // + 0.007880878010262,
                                                   >> 177     dc4 = 1.14023541226785810 - b74 ,   //   0.144711007173263,
                                                   >> 178     dc5 = - 2.3291801924393646 - b75,   // + 0.582357165452555,
                                                   >> 179     dc6 = 1.56887504931661552 - b76 ,   // - 0.458082105929187,
                                                   >> 180     dc7 = 0.025; //- 1.0/66.0 ;
                                                   >> 181     
                                                   >> 182 //    dc1 = -3.0/1280.0,
                                                   >> 183 //    dc2 = 0.0,
                                                   >> 184 //    dc3 = 6561.0/632320.0,
                                                   >> 185 //    dc4 = -343.0/20800.0,
                                                   >> 186 //    dc5 = 243.0/12800.0,
                                                   >> 187 //    dc6 = -1.0/95.0,
                                                   >> 188 //    dc7 = 0.0   ;
                                                   >> 189     
                                                   >> 190     const G4int numberOfVariables= this->GetNumberOfVariables();
173                                                   191     
174     // The number of variables to be integrate    192     // The number of variables to be integrated over
175     //                                         << 
176     yOut[7] = yTemp[7]  = yIn[7] = yInput[7];     193     yOut[7] = yTemp[7]  = yIn[7] = yInput[7];
177                                                << 
178     //  Saving yInput because yInput and yOut     194     //  Saving yInput because yInput and yOut can be aliases for same array
179     //                                         << 195     
180     for(G4int i=0; i<numberOfVariables; ++i)   << 196     for(i=0;i<numberOfVariables;i++)
181     {                                             197     {
182         yIn[i]=yInput[i];                         198         yIn[i]=yInput[i];
183     }                                             199     }
184     // RightHandSide(yIn, dydx) ;   // 1st Ste << 200 
                                                   >> 201     // RightHandSide(yIn, dydx) ;
                                                   >> 202     // 1st Step - Not doing, getting passed
185                                                   203     
186     for(G4int i=0; i<numberOfVariables; ++i)   << 204     for(i=0;i<numberOfVariables;i++)
187     {                                             205     {
188         yTemp[i] = yIn[i] + b21*Step*dydx[i] ;    206         yTemp[i] = yIn[i] + b21*Step*dydx[i] ;
189     }                                             207     }
190     RightHandSide(yTemp, ak2) ;              /    208     RightHandSide(yTemp, ak2) ;              // 2nd Stage
191                                                   209     
192     for(G4int i=0; i<numberOfVariables; ++i)   << 210     for(i=0;i<numberOfVariables;i++)
193     {                                             211     {
194         yTemp[i] = yIn[i] + Step*(b31*dydx[i]     212         yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b32*ak2[i]) ;
195     }                                             213     }
196     RightHandSide(yTemp, ak3) ;              /    214     RightHandSide(yTemp, ak3) ;              // 3rd Stage
197                                                   215     
198     for(G4int i=0; i<numberOfVariables; ++i)   << 216     for(i=0;i<numberOfVariables;i++)
199     {                                             217     {
200         yTemp[i] = yIn[i] + Step*(b41*dydx[i]     218         yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]) ;
201     }                                             219     }
202     RightHandSide(yTemp, ak4) ;              /    220     RightHandSide(yTemp, ak4) ;              // 4th Stage
203                                                   221     
204     for(G4int i=0; i<numberOfVariables; ++i)   << 222     for(i=0;i<numberOfVariables;i++)
205     {                                             223     {
206         yTemp[i] = yIn[i] + Step*(b51*dydx[i]     224         yTemp[i] = yIn[i] + Step*(b51*dydx[i] + b52*ak2[i] + b53*ak3[i] +
207                                   b54*ak4[i])     225                                   b54*ak4[i]) ;
208     }                                             226     }
209     RightHandSide(yTemp, ak5) ;              /    227     RightHandSide(yTemp, ak5) ;              // 5th Stage
210                                                   228     
211     for(G4int i=0; i<numberOfVariables; ++i)   << 229     for(i=0;i<numberOfVariables;i++)
212     {                                             230     {
213         yTemp[i] = yIn[i] + Step*(b61*dydx[i]     231         yTemp[i] = yIn[i] + Step*(b61*dydx[i] + b62*ak2[i] + b63*ak3[i] +
214                                   b64*ak4[i] +    232                                   b64*ak4[i] + b65*ak5[i]) ;
215     }                                             233     }
216     RightHandSide(yTemp, ak6) ;              /    234     RightHandSide(yTemp, ak6) ;              // 6th Stage
217                                                   235     
218     for(G4int i=0; i<numberOfVariables; ++i)   << 236     for(i=0;i<numberOfVariables;i++)
219     {                                             237     {
220         yOut[i] = yIn[i] + Step*(b71*dydx[i] +    238         yOut[i] = yIn[i] + Step*(b71*dydx[i] + b72*ak2[i] + b73*ak3[i] +
221                                  b74*ak4[i] +     239                                  b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
222     }                                             240     }
223     RightHandSide(yOut, ak7);                / << 241     RightHandSide(yOut, ak7);       //7th Stage
224                                                   242 
225     //Calculate the error in the step:            243     //Calculate the error in the step:
226     for(G4int i=0; i<numberOfVariables; ++i)   << 244     for(i=0;i<numberOfVariables;i++)
227     {                                             245     {
228         yErr[i] = Step*(dc1*dydx[i] + dc2*ak2[    246         yErr[i] = Step*(dc1*dydx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
229                         dc5*ak5[i]  + dc6*ak6[    247                         dc5*ak5[i]  + dc6*ak6[i] + dc7*ak7[i] ) ;
230                                                   248         
231         // Store Input and Final values, for p    249         // Store Input and Final values, for possible use in calculating chord
232         //                                     << 
233         fLastInitialVector[i] = yIn[i] ;          250         fLastInitialVector[i] = yIn[i] ;
234         fLastFinalVector[i]   = yOut[i];          251         fLastFinalVector[i]   = yOut[i];
235         fLastDyDx[i]          = dydx[i];          252         fLastDyDx[i]          = dydx[i];
236     }                                             253     }
237                                                   254     
238     fLastStepLength = Step;                       255     fLastStepLength = Step;
239                                                   256     
240     return ;                                      257     return ;
241 }                                                 258 }
242                                                   259 
243 void G4TsitourasRK45::SetupInterpolation()     << 260 void G4TsitourasRK45::SetupInterpolation() // (const G4double *yInput, const G4double *dydx, const G4double Step)
244   // (const G4double *yInput, const G4double * << 
245 {                                                 261 {
246   // Nothing to be done                        << 262     //Nothing to be done
247 }                                                 263 }
248                                                   264 
249 void G4TsitourasRK45::Interpolate(const G4doub << 265 
250                                   const G4doub << 266 void G4TsitourasRK45::Interpolate(const G4double *yInput, const G4double *dydx, const G4double Step, G4double *yOut, G4double tau){
251                                   const G4doub << 267 
252                                         G4doub << 268     
253                                         G4doub << 
254 {                                              << 
255     G4double bf1, bf2, bf3, bf4, bf5, bf6, bf7    269     G4double bf1, bf2, bf3, bf4, bf5, bf6, bf7;
256       // Coefficients for all the seven stages << 270     // Coefficients for all the seven stages.
257                                                   271     
258     const G4int numberOfVariables = GetNumberO << 272     const G4int numberOfVariables= this->GetNumberOfVariables();
259                                                   273     
260     G4double tau0 = tau;                          274     G4double tau0 = tau;
261                                                   275     
262     for(G4int i=0; i<numberOfVariables; ++i)   << 276     for(int i=0;i<numberOfVariables;i++)
263     {                                             277     {
264        yIn[i] = yInput[i];                     << 278         yIn[i]=yInput[i];
265     }                                             279     }
266                                                   280     
267     G4double tau_2 = tau0*tau0 ;               << 281     G4double
268        //    tau_3 = tau0*tau_2,               << 282     tau_2 = tau0*tau0 ;
269        //    tau_4 = tau_2*tau_2;              << 283 //    tau_3 = tau0*tau_2,
                                                   >> 284 //    tau_4 = tau_2*tau_2;
270                                                   285 
271     bf1 = -1.0530884977290216*tau*(tau - 1.329    286     bf1 = -1.0530884977290216*tau*(tau - 1.3299890189751412)*(tau_2 -
272                 1.4364028541716351*tau + 0.713 << 287                 1.4364028541716351*tau + 0.7139816917074209),
273     bf2 = 0.1017*tau_2*(tau_2 - 2.196656833824    288     bf2 = 0.1017*tau_2*(tau_2 - 2.1966568338249754*tau +
274                         1.2949852507374631);   << 289                         1.2949852507374631),
275     bf3 = 2.490627285651252793*tau_2*(tau_2 -     290     bf3 = 2.490627285651252793*tau_2*(tau_2 - 2.38535645472061657*tau
276                                       + 1.5780 << 291                                       + 1.57803468208092486) ,
277     bf4 = -16.54810288924490272*(tau - 1.21712    292     bf4 = -16.54810288924490272*(tau - 1.21712927295533244)*
278                                     (tau - 0.6 << 293             (tau - 0.61620406037800089)*tau_2,
279     bf5 = 47.37952196281928122*(tau - 1.203071    294     bf5 = 47.37952196281928122*(tau - 1.203071208372362603)*
280                                     (tau - 0.6 << 295             (tau - 0.658047292653547382)*tau_2,
281     bf6 = -34.87065786149660974*(tau - 1.2)*(t    296     bf6 = -34.87065786149660974*(tau - 1.2)*(tau -
282                                 0.666666666666 << 297                                 0.666666666666666667)*tau_2,
283     bf7 = 2.5*(tau - 1.0)*(tau - 0.6)*tau_2;      298     bf7 = 2.5*(tau - 1.0)*(tau - 0.6)*tau_2;
284                                                   299     
285     // Putting together the coefficients calcu << 300     //Putting together the coefficients calculated as the respective stage coefficients
286     // stage coefficients                      << 301     for( int i=0; i<numberOfVariables; i++){
287     //                                         << 302         yOut[i] = yIn[i] + Step*(  bf1*dydx[i] + bf2*ak2[i] + bf3*ak3[i] + bf4*ak4[i]
288     for(G4int i=0; i<numberOfVariables; ++i)   << 303                                  + bf5*ak5[i]  + bf6*ak6[i] + bf7*ak7[i]  ) ;
289     {                                          << 
290         yOut[i] = yIn[i] + Step*(  bf1*dydx[i] << 
291                                  + bf4*ak4[i]  << 
292                                  + bf7*ak7[i]  << 
293     }                                             304     }
294 }                                                 305 }
295                                                   306 
                                                   >> 307 ///////////////////////////////////////////////////////////////////////////////
                                                   >> 308 
296 G4double  G4TsitourasRK45::DistChord() const      309 G4double  G4TsitourasRK45::DistChord() const
297 {                                                 310 {
298   G4double distLine, distChord;                   311   G4double distLine, distChord; 
299   G4ThreeVector initialPoint, finalPoint, midP    312   G4ThreeVector initialPoint, finalPoint, midPoint;
300                                                   313 
301   // Store last initial and final points (they << 314   // Store last initial and final points (they will be overwritten in self-Stepper call!)
302   // overwritten in self-Stepper call!)        << 
303   //                                           << 
304   initialPoint = G4ThreeVector( fLastInitialVe    315   initialPoint = G4ThreeVector( fLastInitialVector[0], 
305                                 fLastInitialVe    316                                 fLastInitialVector[1], fLastInitialVector[2]); 
306   finalPoint   = G4ThreeVector( fLastFinalVect    317   finalPoint   = G4ThreeVector( fLastFinalVector[0],  
307                                 fLastFinalVect    318                                 fLastFinalVector[1],  fLastFinalVector[2]); 
308                                                   319 
309   // Do half a step using StepNoErr               320   // Do half a step using StepNoErr
310   //                                           << 321 
311   fAuxStepper->Stepper( fLastInitialVector, fL    322   fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength, 
312            fMidVector,   fMidError );             323            fMidVector,   fMidError );
313                                                   324 
314   midPoint = G4ThreeVector( fMidVector[0], fMi    325   midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);       
315                                                   326 
316   // Use stored values of Initial and Endpoint    327   // Use stored values of Initial and Endpoint + new Midpoint to evaluate
317   //  distance of Chord                           328   //  distance of Chord
318   //                                           << 329 
                                                   >> 330 
319   if (initialPoint != finalPoint)                 331   if (initialPoint != finalPoint) 
320   {                                               332   {
321      distLine  = G4LineSection::Distline( midP    333      distLine  = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
322      distChord = distLine;                        334      distChord = distLine;
323   }                                               335   }
324   else                                            336   else
325   {                                               337   {
326      distChord = (midPoint-initialPoint).mag()    338      distChord = (midPoint-initialPoint).mag();
327   }                                               339   }
328   return distChord;                               340   return distChord;
329 }                                                 341 }
330                                                   342