Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:59:09 

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