Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-09-17 09:14:11 

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // G4TDormandPrince45
0027 //
0028 // Class desription:
0029 //
0030 //  An implementation of the 5th order embedded RK method from the paper:
0031 //  J. R. Dormand and P. J. Prince, "A family of embedded Runge-Kutta formulae"
0032 //  Journal of computational and applied Math., vol.6, no.1, pp.19-26, 1980.
0033 //
0034 //  DormandPrince7 - 5(4) embedded RK method
0035 
0036 // Author: Josh Xie (CERN, Google Summer of Code 2014), June 2014
0037 // Supervisors:  Sandro Wenzel, John Apostolakis (CERN)
0038 // --------------------------------------------------------------------
0039 #ifndef G4TDORMAND_PRINCE_45_HH
0040 #define G4TDORMAND_PRINCE_45_HH
0041 
0042 #include "G4MagIntegratorStepper.hh"
0043 #include "G4FieldUtils.hh"
0044 #include "G4LineSection.hh"
0045 
0046 #include <cstring>
0047 #include <cassert>
0048 
0049 /**
0050  * @brief G4TDormandPrince45 is a templated version of G4DormandPrince745
0051  * 5th order Runge-Kutta stepper.
0052  */
0053 
0054 template <class T_Equation, unsigned int N = 6 >
0055 class G4TDormandPrince45 : public G4MagIntegratorStepper
0056 {
0057   public:
0058 
0059     G4TDormandPrince45(T_Equation* equation );
0060     G4TDormandPrince45(T_Equation* equation, G4int numVar ); // must have numVar == N
0061    
0062     inline void StepWithError(const G4double yInput[],
0063                               const G4double dydx[],
0064                                     G4double hstep,
0065                                     G4double yOutput[],
0066                                     G4double yError[] ) ;
0067    
0068     void Stepper(const G4double yInput[],
0069                  const G4double dydx[],
0070                        G4double hstep,
0071                        G4double yOutput[],
0072                        G4double yError[]) final;
0073 
0074     inline void StepWithFinalDerivate(const G4double yInput[],   
0075                                       const G4double dydx[],
0076                                             G4double hstep,
0077                                             G4double yOutput[],
0078                                             G4double yError[],
0079                                             G4double dydxOutput[]);
0080    
0081     inline void SetupInterpolation() {}
0082 
0083     void Interpolate(G4double tau, G4double yOut[]) const
0084     {
0085       Interpolate4thOrder(yOut, tau);
0086     }
0087       // For calculating the output at the tau fraction of Step
0088 
0089     G4double DistChord() const final;
0090 
0091     inline G4int IntegratorOrder() const override { return 4; }
0092 
0093     G4StepperType StepperType() const override { return kTDormandPrince45; }
0094 
0095     inline const field_utils::ShortState<N>& GetYOut() const { return fyOut; }
0096 
0097     void Interpolate4thOrder(G4double yOut[], G4double tau) const;
0098 
0099     void SetupInterpolation5thOrder();
0100     void Interpolate5thOrder(G4double yOut[], G4double tau) const;
0101   
0102     // __attribute__((always_inline))
0103     inline void RightHandSideInl( const G4double y[],
0104                                         G4double dydx[] )
0105     {
0106       fEquation_Rhs->T_Equation::RightHandSide(y, dydx);
0107     }
0108 
0109     inline void Stepper(const G4double yInput[],
0110                         const G4double dydx[],
0111                               G4double hstep, G4double yOutput[],
0112                               G4double yError[], G4double dydxOutput[])
0113     {
0114       StepWithFinalDerivate(yInput, dydx, hstep, yOutput, yError, dydxOutput);
0115     }
0116 
0117     T_Equation* GetSpecificEquation() { return fEquation_Rhs; }
0118    
0119     static constexpr G4int N8 = N > 8 ? N : 8;  //  y[
0120    
0121   private:
0122 
0123     field_utils::ShortState<N>  ak2, ak3, ak4, ak5, ak6, ak7, ak8, ak9; 
0124     field_utils::ShortState<N8> fyIn;
0125     field_utils::ShortState<N>  fyOut, fdydxIn;
0126 
0127     // - Simpler :
0128     // field_utils::State ak2, ak3, ak4, ak5, ak6, ak7, ak8, ak9;
0129     // field_utils::State fyIn, fyOut, fdydxIn;
0130 
0131     G4double fLastStepLength = -1.0;
0132     T_Equation* fEquation_Rhs;
0133 };
0134 
0135 // --------------------------------------------------------------------
0136 // G4TDormandPrince745 implementation -- borrowed from G4DormandPrince745
0137 //
0138 // DormandPrince7 - 5(4) non-FSAL
0139 // definition of the stepper() method that evaluates one step in
0140 // field propagation.
0141 // The coefficients and the algorithm have been adapted from
0142 //
0143 // J. R. Dormand and P. J. Prince, "A family of embedded Runge-Kutta formulae"
0144 // Journal of computational and applied Math., vol.6, no.1, pp.19-26, 1980.
0145 //
0146 // The Butcher table of the Dormand-Prince-7-4-5 method is as follows :
0147 //
0148 //    0   |
0149 //    1/5 | 1/5
0150 //    3/10| 3/40       9/40
0151 //    4/5 | 44/45      56/15      32/9
0152 //    8/9 | 19372/6561 25360/2187 64448/6561  212/729
0153 //    1   | 9017/3168  355/33     46732/5247  49/176  5103/18656
0154 //    1   | 35/384     0          500/1113    125/192 2187/6784    11/84
0155 //    ------------------------------------------------------------------------
0156 //          35/384     0          500/1113    125/192 2187/6784    11/84    0
0157 //          5179/57600 0          7571/16695  393/640 92097/339200 187/2100 1/40
0158 //
0159 // --------------------------------------------------------------------
0160 
0161 // Constructor
0162 //
0163 template <class T_Equation, unsigned int N>
0164 G4TDormandPrince45<T_Equation,N>::G4TDormandPrince45(T_Equation* equation )
0165     : G4MagIntegratorStepper(dynamic_cast<G4EquationOfMotion*>(equation), N )
0166     , fEquation_Rhs(equation)
0167 {
0168   // assert( dynamic_cast<G4EquationOfMotion*>(equation) != nullptr );
0169   if( dynamic_cast<G4EquationOfMotion*>(equation) == nullptr )
0170   {
0171     G4Exception("G4TDormandPrince745: constructor", "GeomField0001",
0172                 FatalException, "T_Equation is not an G4EquationOfMotion.");
0173   }
0174 
0175   /***
0176   assert( equation->GetNumberOfVariables == N );
0177   if( equation->GetNumberOfVariables != N ){
0178     G4ExceptionDescription msg;
0179     msg << "Equation has an incompatible number of variables." ;
0180     msg << "   template N = " << N << " equation-Nvar= "
0181         << equation->GetNumberOfVariables;
0182     G4Exception("G4TCashKarpRKF45: constructor", "GeomField0001",
0183                 FatalException, msg );
0184    } ****/
0185 }
0186 
0187 template <class T_Equation, unsigned int N>
0188 inline G4TDormandPrince45<T_Equation,N>::
0189        G4TDormandPrince45(T_Equation* equation, G4int numVar )
0190   : G4TDormandPrince45<T_Equation,N>(equation )
0191 {
0192   if( numVar != G4int(N))
0193   {
0194     G4ExceptionDescription msg;
0195     msg << "Equation has an incompatible number of variables." ;
0196     msg << "   template N = " << N
0197         << "   argument numVar = " << numVar ;
0198     //    << " equation-Nvar= " << equation->GetNumberOfVariables(); // --> Expected later
0199     G4Exception("G4TCashKarpRKF45: constructor", "GeomField0001",
0200                    FatalErrorInArgument, msg );
0201   }
0202   assert( numVar == N );
0203 }
0204 
0205 template <class T_Equation, unsigned int N>
0206 inline void
0207 G4TDormandPrince45<T_Equation,N>::StepWithFinalDerivate(const G4double yInput[],
0208                                                         const G4double dydx[],
0209                                                               G4double hstep,
0210                                                               G4double yOutput[],
0211                                                               G4double yError[],
0212                                                               G4double dydxOutput[])
0213 {
0214   StepWithError(yInput, dydx, hstep, yOutput, yError);
0215   field_utils::copy(dydxOutput, ak7, N);
0216 }
0217 
0218 // Stepper
0219 //
0220 // Passing in the value of yInput[],the first time dydx[] and Step length
0221 // Giving back yOut and yErr arrays for output and error respectively
0222 //
0223 
0224 template <class T_Equation, unsigned int N>
0225 inline void
0226 G4TDormandPrince45<T_Equation,N>::StepWithError(const G4double yInput[],
0227                                                 const G4double dydx[],
0228                                                       G4double hstep,
0229                                                       G4double yOut[],
0230                                                       G4double yErr[] )
0231 {
0232   // The parameters of the Butcher tableu
0233   //
0234   constexpr G4double b21 = 0.2,
0235                  b31 = 3.0 / 40.0, b32 = 9.0 / 40.0,
0236                  b41 = 44.0 / 45.0, b42 = -56.0 / 15.0, b43 = 32.0/9.0,
0237 
0238   b51 = 19372.0 / 6561.0, b52 = -25360.0 / 2187.0, b53 = 64448.0 / 6561.0,
0239   b54 = -212.0 / 729.0,
0240         
0241   b61 = 9017.0 / 3168.0 , b62 =   -355.0 / 33.0,
0242   b63 = 46732.0 / 5247.0, b64 = 49.0 / 176.0,
0243   b65 = -5103.0 / 18656.0,
0244         
0245   b71 = 35.0 / 384.0, b72 = 0.,
0246   b73 = 500.0 / 1113.0, b74 = 125.0 / 192.0,
0247   b75 = -2187.0 / 6784.0, b76 = 11.0 / 84.0,
0248     
0249   // Sum of columns, sum(bij) = ei
0250   //    e1 = 0. ,
0251   //    e2 = 1.0/5.0 ,
0252   //    e3 = 3.0/10.0 ,
0253   //    e4 = 4.0/5.0 ,
0254   //    e5 = 8.0/9.0 ,
0255   //    e6 = 1.0 ,
0256   //    e7 = 1.0 ,
0257   
0258   // Difference between the higher and the lower order method coeff. :
0259   // b7j are the coefficients of higher order
0260     
0261   dc1 = -(b71 - 5179.0 / 57600.0),
0262   dc2 = -(b72 - .0),
0263   dc3 = -(b73 - 7571.0 / 16695.0),
0264   dc4 = -(b74 - 393.0 / 640.0),
0265   dc5 = -(b75 + 92097.0 / 339200.0),
0266   dc6 = -(b76 - 187.0 / 2100.0),
0267   dc7 = -(- 1.0 / 40.0);
0268     
0269   // const G4int numberOfVariables = GetNumberOfVariables();
0270   //   The number of variables to be integrated over    
0271   field_utils::ShortState<N8> yTemp;
0272     
0273   yOut[7] = yTemp[7]  = fyIn[7] = yInput[7];  // Pass along the time - used in RightHandSide
0274     
0275   //  Saving yInput because yInput and yOut can be aliases for same array
0276   //
0277   for(unsigned int i = 0; i < N; ++i)
0278   {
0279     fyIn[i]  = yInput[i];
0280     yTemp[i] = yInput[i] + b21 * hstep * dydx[i];
0281   }
0282   RightHandSideInl(yTemp, ak2);              // 2nd stage
0283 
0284   for(unsigned int i = 0; i < N; ++i)
0285   {
0286     yTemp[i] = fyIn[i] + hstep * (b31 * dydx[i] + b32 * ak2[i]);
0287   }
0288   RightHandSideInl(yTemp, ak3);              // 3rd stage
0289     
0290   for(unsigned int i = 0; i < N; ++i)
0291   {
0292     yTemp[i] = fyIn[i] + hstep * (
0293         b41 * dydx[i] + b42 * ak2[i] + b43 * ak3[i]);
0294   }
0295   RightHandSideInl(yTemp, ak4);              // 4th stage
0296     
0297   for(unsigned int i = 0; i < N; ++i)
0298   {
0299     yTemp[i] = fyIn[i] + hstep * (
0300         b51 * dydx[i] + b52 * ak2[i] + b53 * ak3[i] + b54 * ak4[i]);
0301   }
0302   RightHandSideInl(yTemp, ak5);              // 5th stage
0303     
0304   for(unsigned int i = 0; i < N; ++i)
0305   {
0306     yTemp[i] = fyIn[i] + hstep * (
0307         b61 * dydx[i] + b62 * ak2[i] + 
0308         b63 * ak3[i] + b64 * ak4[i] + b65 * ak5[i]);
0309   }
0310   RightHandSideInl(yTemp, ak6);              // 6th stage
0311     
0312   for(unsigned int i = 0; i < N; ++i)
0313   {
0314     yOut[i] = fyIn[i] + hstep * (
0315         b71 * dydx[i] + b72 * ak2[i] + b73 * ak3[i] +
0316         b74 * ak4[i] + b75 * ak5[i] + b76 * ak6[i]);
0317   }
0318   RightHandSideInl(yOut, ak7);               // 7th and Final stage
0319     
0320   for(unsigned int i = 0; i < N; ++i)
0321   {
0322     yErr[i] = hstep * (
0323         dc1 * dydx[i] + dc2 * ak2[i] + 
0324         dc3 * ak3[i] + dc4 * ak4[i] +
0325         dc5 * ak5[i] + dc6 * ak6[i] + dc7 * ak7[i]
0326     ) + 1.5e-18;
0327 
0328     // Store Input and Final values, for possible use in calculating chord
0329     //
0330     fyOut[i] = yOut[i];
0331     fdydxIn[i] = dydx[i];        
0332   }
0333     
0334   fLastStepLength = hstep;
0335 }
0336 
0337 template <class T_Equation, unsigned int N >
0338 inline void
0339 G4TDormandPrince45<T_Equation,N>::Stepper(const G4double yInput[],
0340                                           const G4double dydx[],
0341                                                 G4double Step,
0342                                                 G4double yOutput[],
0343                                                 G4double yError[])
0344 {
0345   assert( yOutput != yInput );
0346   assert( yError  != yInput );
0347   
0348   StepWithError( yInput, dydx, Step, yOutput, yError);
0349 }
0350 
0351 template <class T_Equation, unsigned int N>
0352 inline G4double G4TDormandPrince45<T_Equation,N>::DistChord() const
0353 {
0354   // Coefficients were taken from Some Practical Runge-Kutta Formulas
0355   // by Lawrence F. Shampine, page 149, c*
0356   //
0357   const G4double hf1 = 6025192743.0 / 30085553152.0,
0358                  hf3 = 51252292925.0 / 65400821598.0,
0359                  hf4 = - 2691868925.0 / 45128329728.0,
0360                  hf5 = 187940372067.0 / 1594534317056.0,
0361                  hf6 = - 1776094331.0 / 19743644256.0,
0362                  hf7 = 11237099.0 / 235043384.0;
0363 
0364   G4ThreeVector mid;
0365 
0366   for(unsigned int i = 0; i < 3; ++i) 
0367   {
0368     mid[i] = fyIn[i] + 0.5 * fLastStepLength * (
0369         hf1 * fdydxIn[i] + hf3 * ak3[i] + 
0370         hf4 * ak4[i] + hf5 * ak5[i] + hf6 * ak6[i] + hf7 * ak7[i]);
0371   }
0372 
0373   const G4ThreeVector begin = makeVector(fyIn, field_utils::Value3D::Position);
0374   const G4ThreeVector end = makeVector(fyOut,  field_utils::Value3D::Position);
0375 
0376   return G4LineSection::Distline(mid, begin, end);
0377 }
0378 
0379 // The lower (4th) order interpolant given by Dormand and Prince:
0380 //        J. R. Dormand and P. J. Prince, "Runge-Kutta triples"
0381 //        Computers & Mathematics with Applications, vol. 12, no. 9,
0382 //        pp. 1007-1017, 1986.
0383 //
0384 template <class T_Equation, unsigned int N>
0385 inline void
0386 G4TDormandPrince45<T_Equation,N>::Interpolate4thOrder(G4double yOut[],
0387                                                       G4double tau) const
0388 {    
0389   const G4double tau2 = tau * tau,
0390                  tau3 = tau * tau2,
0391                  tau4 = tau2 * tau2;
0392 
0393   const G4double bf1 = 1.0 / 11282082432.0 * (
0394       157015080.0 * tau4 - 13107642775.0 * tau3 + 34969693132.0 * tau2 - 
0395       32272833064.0 * tau + 11282082432.0);
0396 
0397   const G4double bf3 = - 100.0 / 32700410799.0 * tau * (
0398       15701508.0 * tau3 - 914128567.0 * tau2 + 2074956840.0 * tau - 
0399       1323431896.0);
0400     
0401   const G4double bf4 = 25.0 / 5641041216.0 * tau * (
0402       94209048.0 * tau3 - 1518414297.0 * tau2 + 2460397220.0 * tau - 
0403       889289856.0);
0404     
0405   const G4double bf5 = - 2187.0 / 199316789632.0 * tau * (
0406       52338360.0 * tau3 - 451824525.0 * tau2 + 687873124.0 * tau - 
0407       259006536.0);
0408 
0409   const G4double bf6 = 11.0 / 2467955532.0 * tau * (
0410       106151040.0 * tau3 - 661884105.0 * tau2 + 
0411       946554244.0 * tau - 361440756.0);
0412     
0413   const G4double bf7 = 1.0 / 29380423.0 * tau * (1.0 - tau) * (
0414       8293050.0 * tau2 - 82437520.0 * tau + 44764047.0);
0415 
0416   for(unsigned int i = 0; i < N; ++i)
0417   {
0418       yOut[i] = fyIn[i] + fLastStepLength * tau * (
0419           bf1 * fdydxIn[i] + bf3 * ak3[i] + bf4 * ak4[i] + 
0420           bf5 * ak5[i] + bf6 * ak6[i] + bf7 * ak7[i]);
0421   }
0422 }
0423 
0424 // Following interpolant of order 5 was given by Baker,Dormand,Gilmore, Prince :
0425 //        T. S. Baker, J. R. Dormand, J. P. Gilmore, and P. J. Prince,
0426 //        "Continuous approximation with embedded Runge-Kutta methods"
0427 //        Applied Numerical Mathematics, vol. 22, no. 1, pp. 51-62, 1996.
0428 //
0429 // Calculating the extra stages for the interpolant
0430 //
0431 template <class T_Equation, unsigned int N>
0432 inline void G4TDormandPrince45<T_Equation,N>::SetupInterpolation5thOrder()
0433 {
0434   // Coefficients for the additional stages
0435   //
0436   const G4double b81 = 6245.0 / 62208.0,
0437                  b82 = 0.0,
0438                  b83 = 8875.0 / 103032.0,
0439                  b84 = -125.0 / 1728.0,
0440                  b85 = 801.0 / 13568.0,
0441                  b86 = -13519.0 / 368064.0,
0442                  b87 = 11105.0 / 368064.0,
0443 
0444                  b91 = 632855.0 / 4478976.0,
0445                  b92 = 0.0,
0446                  b93 = 4146875.0 / 6491016.0,
0447                  b94 = 5490625.0 /14183424.0,
0448                  b95 = -15975.0 / 108544.0,
0449                  b96 = 8295925.0 / 220286304.0,
0450                  b97 = -1779595.0 / 62938944.0,
0451                  b98 = -805.0 / 4104.0;
0452 
0453   field_utils::ShortState<N> yTemp;
0454 
0455   // Evaluate the extra stages
0456   //
0457   for(unsigned int i = 0; i < N; ++i)
0458   {
0459     yTemp[i] = fyIn[i] + fLastStepLength * (
0460         b81 * fdydxIn[i] + b82 * ak2[i] + b83 * ak3[i] +
0461         b84 * ak4[i] + b85 * ak5[i] + b86 * ak6[i] +
0462         b87 * ak7[i]
0463     );
0464   }
0465   RightHandSideInl(yTemp, ak8);          // 8th Stage
0466     
0467   for(unsigned int i = 0; i < N; ++i)
0468   {
0469     yTemp[i] = fyIn[i] + fLastStepLength * (
0470         b91 * fdydxIn[i] + b92 * ak2[i] + b93 * ak3[i] +
0471         b94 * ak4[i] + b95 * ak5[i] + b96 * ak6[i] +
0472         b97 * ak7[i] + b98 * ak8[i]
0473     );
0474   }
0475   RightHandSideInl(yTemp, ak9);          // 9th Stage
0476 }
0477 
0478 // Calculating the interpolated result yOut with the coefficients
0479 //
0480 template <class T_Equation, unsigned int N>
0481 inline void G4TDormandPrince45<T_Equation,N>::
0482 Interpolate5thOrder(G4double yOut[], G4double tau) const
0483 {
0484   // Define the coefficients for the polynomials
0485   //
0486   G4double bi[10][5];
0487     
0488   //  COEFFICIENTS OF   bi[1]
0489   bi[1][0] = 1.0,
0490   bi[1][1] = -38039.0 / 7040.0,
0491   bi[1][2] = 125923.0 / 10560.0,
0492   bi[1][3] = -19683.0 / 1760.0,
0493   bi[1][4] = 3303.0 / 880.0,
0494   //  --------------------------------------------------------
0495   //
0496   //  COEFFICIENTS OF  bi[2]
0497   bi[2][0] = 0.0,
0498   bi[2][1] = 0.0,
0499   bi[2][2] = 0.0,
0500   bi[2][3] = 0.0,
0501   bi[2][4] = 0.0,
0502   //  --------------------------------------------------------
0503   //
0504   //  COEFFICIENTS OF  bi[3]
0505   bi[3][0] = 0.0,
0506   bi[3][1] = -12500.0 / 4081.0,
0507   bi[3][2] = 205000.0 / 12243.0,
0508   bi[3][3] = -90000.0 / 4081.0,
0509   bi[3][4] = 36000.0 / 4081.0,
0510   //  --------------------------------------------------------
0511   //
0512   //  COEFFICIENTS OF  bi[4]
0513   bi[4][0] = 0.0,
0514   bi[4][1] = -3125.0 / 704.0,
0515   bi[4][2] = 25625.0 / 1056.0,
0516   bi[4][3] = -5625.0 / 176.0,
0517   bi[4][4] = 1125.0 / 88.0,
0518   //  --------------------------------------------------------
0519   //
0520   //  COEFFICIENTS OF  bi[5]
0521   bi[5][0] = 0.0,
0522   bi[5][1] = 164025.0 / 74624.0,
0523   bi[5][2] = -448335.0 / 37312.0,
0524   bi[5][3] = 295245.0 / 18656.0,
0525   bi[5][4] = -59049.0 / 9328.0,
0526   //  --------------------------------------------------------
0527   //
0528   //  COEFFICIENTS OF  bi[6]
0529   bi[6][0] = 0.0,
0530   bi[6][1] = -25.0 / 28.0,
0531   bi[6][2] = 205.0 / 42.0,
0532   bi[6][3] = -45.0 / 7.0,
0533   bi[6][4] = 18.0 / 7.0,
0534   //  --------------------------------------------------------
0535   //
0536   //  COEFFICIENTS OF  bi[7]
0537   bi[7][0] = 0.0,
0538   bi[7][1] = -2.0 / 11.0,
0539   bi[7][2] = 73.0 / 55.0,
0540   bi[7][3] = -171.0 / 55.0,
0541   bi[7][4] = 108.0 / 55.0,
0542   //  --------------------------------------------------------
0543   //
0544   //  COEFFICIENTS OF  bi[8]
0545   bi[8][0] = 0.0,
0546   bi[8][1] = 189.0 / 22.0,
0547   bi[8][2] = -1593.0 / 55.0,
0548   bi[8][3] = 3537.0 / 110.0,
0549   bi[8][4] = -648.0 / 55.0,
0550   //  --------------------------------------------------------
0551   //
0552   //  COEFFICIENTS OF  bi[9]
0553   bi[9][0] = 0.0,
0554   bi[9][1] = 351.0 / 110.0,
0555   bi[9][2] = -999.0 / 55.0,
0556   bi[9][3] = 2943.0 / 110.0,
0557   bi[9][4] = -648.0 / 55.0;
0558   //  --------------------------------------------------------
0559         
0560   // Calculating the polynomials
0561  
0562   G4double b[10];
0563   std::memset(b, 0.0, sizeof(b));
0564 
0565   G4double tauPower = 1.0;
0566   for(G4int j = 0; j <= 4; ++j)
0567   {       
0568     for(G4int iStage = 1; iStage <= 9; ++iStage)
0569     {
0570       b[iStage] += bi[iStage][j] * tauPower;
0571     }
0572     tauPower *= tau;       
0573   }
0574 
0575   const G4double stepLen = fLastStepLength * tau;
0576   for(G4int i = 0; i < N; ++i)
0577   {
0578     yOut[i] = fyIn[i] + stepLen * (
0579        b[1] * fdydxIn[i] + b[2] * ak2[i] + b[3] * ak3[i] +
0580        b[4] * ak4[i] + b[5] * ak5[i] + b[6] * ak6[i] +
0581        b[7] * ak7[i] + b[8] * ak8[i] + b[9] * ak9[i]
0582     );
0583   }
0584 }
0585 
0586 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!