include/Geant4/G4TSimpleRunge.hh

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 // G4TSimpleRunge
0027 //
0028 // Class description:
0029 //
0030 // Templated version of G4SimpleRunge.
0031 // Adapted from G4G4TSimpleRunge class.
0032
0033 // Author: Josh Xie (CERN, Google Summer of Code 2014), June 2014
0034 // Supervisors:  Sandro Wenzel, John Apostolakis (CERN)
0035 // --------------------------------------------------------------------
0036 #ifndef G4TSimpleRunge_HH
0037 #define G4TSimpleRunge_HH
0038
0039 #include <cassert>
0040 #include "G4TMagErrorStepper.hh"
0041 #include "G4ThreeVector.hh"
0042
0043 /**
0044  * @brief G4TSimpleRunge is a templated version of G4SimpleRunge.
0045  */
0046
0047 template <class T_Equation, int N>
0048 class G4TSimpleRunge
0049   : public G4TMagErrorStepper<G4TSimpleRunge<T_Equation, N>, T_Equation, N>
0050 {
0051   public:
0052
0053     static constexpr double IntegratorCorrection = 1. / ((1 << 2) - 1);
0054
0055     G4TSimpleRunge(T_Equation* EqRhs, G4int numberOfVariables = 6)
0056       : G4TMagErrorStepper<G4TSimpleRunge<T_Equation, N>, T_Equation, N>(
0057           EqRhs, numberOfVariables)
0058       , fNumberOfVariables(numberOfVariables)
0059       , fEquation_Rhs(EqRhs)
0060
0061     {
0062       // default GetNumberOfStateVariables() == 12
0063       assert(this->GetNumberOfStateVariables() <= 12);
0064     }
0065
0066     ~G4TSimpleRunge() = default;
0067
0068     inline void RightHandSide(G4double y[],
0069                               G4double dydx[])
0070     {
0071       fEquation_Rhs->T_Equation::RightHandSide(y, dydx);
0072     }
0073
0074     inline void DumbStepper(const G4double yIn[],
0075                             const G4double dydx[],
0076                             G4double h, G4double yOut[]) // override final
0077     {
0078       // Initialise time to t0, needed when it is not updated by the integration.
0079       yTemp[7] = yOut[7] = yIn[7];  //  Better to set it to NaN;  // TODO
0080
0081       for(G4int i = 0; i < N; ++i)
0082       {
0083         yTemp[i] = yIn[i] + 0.5 * h * dydx[i];
0084       }
0085
0086       this->RightHandSide(yTemp, dydxTemp);
0087
0088       for(G4int i = 0; i < N; ++i)
0089       {
0090         yOut[i] = yIn[i] + h * (dydxTemp[i]);
0091       }
0092     }
0093
0094     inline G4int IntegratorOrder() const { return 2; }
0095
0096   private:
0097
0098     G4int fNumberOfVariables;
0099     G4double dydxTemp[N > 12 ? N : 12];
0100     G4double yTemp[N > 12 ? N : 12];
0101
0102     T_Equation* fEquation_Rhs;
0103     // scratch space
0104 };
0105
0106 #endif