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