Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // 27 // 28 /// \file ODESolver.hh 29 /// \brief Definition of the ODESolver class 30 31 #ifndef ODESolver_h 32 #define ODESolver_h 1 33 34 #include <vector> 35 #include <map> 36 #include <functional> 37 38 39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 40 41 using myODEs = std::vector<double>(*)(double ,std::vector<double>) ; 42 std::vector<double> operator*(const std::vector<double> v, double alfa); 43 std::vector<double> operator+(const std::vector<double> v, double alfa); 44 std::vector<double> operator+(const std::vector<double> v1, const std::vector<double> v2); 45 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 47 48 class ODESolver 49 { 50 public: 51 ODESolver(); 52 ~ODESolver() = default; 53 void Embedded_RungeKutta_Fehlberg( 54 std::function<std::vector<double>(double,std::vector<double>)>, std::vector<double> &y, 55 double start,double end,double stepsize=-1, double epsilon = 1e-3, 56 std::vector<double> *time_observer=nullptr,std::vector<std::vector<double>> *state_observer=nullptr); 57 void SetNstepsForObserver(unsigned int ndt) {fNstepsForObserver = ndt;} 58 void RungeKutta4( 59 std::function<std::vector<double>(double,std::vector<double>)>, std::vector<double> &y, 60 double start,double end,double stepsize=-1, 61 std::vector<double> *time_observer=nullptr,std::vector<std::vector<double>> *state_observer=nullptr); 62 private: 63 double RungeKutta_Fehlberg(std::function<std::vector<double>(double,std::vector<double>)>, 64 std::vector<double> &y,double t, double stepsize=-1); 65 void absValuesVector(std::vector<double> &vIn); 66 unsigned int fNstepsForObserver{1}; 67 }; 68 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 70 71 inline void ODESolver::absValuesVector(std::vector<double> &vIn) 72 { 73 for (double &val : vIn) { 74 if (val < 0) val *= -1.; 75 } 76 } 77 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 79 80 #endif