00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #ifndef RKIntegrator_h
00030 #define RKIntegrator_h 1
00031 #include "CLHEP/GenericFunctions/AbsFunction.hh"
00032 #include "CLHEP/GenericFunctions/Parameter.hh"
00033 #include "CLHEP/GenericFunctions/RCBase.hh"
00034 #include <vector>
00035 #include <set>
00036 namespace Genfun {
00037
00043 class RKIntegrator {
00044
00045 public:
00046
00047
00048 class RKFunction;
00049 class RKData;
00050
00051
00052 RKIntegrator();
00053
00054
00055 virtual ~RKIntegrator();
00056
00057
00058
00059
00060
00061 Parameter * addDiffEquation (const AbsFunction * diffEquation,
00062 const std::string & variableName="anon",
00063 double defStartingValue=0.0,
00064 double startingValueMin=0.0,
00065 double startingValueMax=0.0);
00066
00067
00068
00069
00070 Parameter *createControlParameter (const std::string & variableName="anon",
00071 double defStartingValue=0.0,
00072 double startingValueMin=0.0,
00073 double startingValueMax=0.0);
00074
00075
00076
00077
00078 const RKFunction *getFunction(unsigned int i) const;
00079
00080
00081 private:
00082
00083
00084 const RKIntegrator & operator=(const RKIntegrator &right);
00085
00086
00087 RKIntegrator(const RKIntegrator &right);
00088
00089
00090
00091 RKData *_data;
00092
00093
00094
00095 std::vector<const RKFunction *> _fcn;
00096
00097
00098 };
00099
00100
00101 class RKIntegrator::RKData : public Genfun::RCBase {
00102
00103
00104 public:
00105
00106 struct Data{
00107
00108 std::vector<double> variable;
00109 mutable std::vector<double> firstDerivative;
00110 double time;
00111 mutable bool dcalc;
00112
00113 Data(int size): variable(size), firstDerivative(size), time(0), dcalc(false) {}
00114 bool operator < (const Data & right) const { return time < right.time; }
00115 bool operator == (const Data & right) const { return time==right.time; }
00116 };
00117
00118 RKData();
00119 void lock();
00120 void recache();
00121
00122 std::vector<Parameter *> _startingValParameter;
00123 std::vector<double> _startingValParameterCache;
00124
00125 std::vector <Parameter *> _controlParameter;
00126 std::vector <double> _controlParameterCache;
00127
00128 std::vector<const AbsFunction *> _diffEqn;
00129 std::set<Data > _fx;
00130 bool _locked;
00131
00132 private:
00133
00134 ~RKData();
00135 friend class ImaginaryFriend;
00136
00137 };
00138
00139 class RKIntegrator::RKFunction : public AbsFunction {
00140
00141 FUNCTION_OBJECT_DEF(RKFunction)
00142
00143 public:
00144
00145
00146 RKFunction(RKData *data, unsigned int index);
00147
00148
00149 virtual ~RKFunction();
00150
00151
00152 RKFunction(const RKFunction &right);
00153
00154
00155 virtual double operator ()(double argument) const;
00156 virtual double operator ()(const Argument & a) const {return operator() (a[0]);}
00157
00158 private:
00159
00160
00161 const RKFunction & operator=(const RKFunction &right);
00162
00163
00164 RKData *_data;
00165 const unsigned int _index;
00166
00167 void rk4 (const RKData::Data & sdata, RKData::Data & ddata) const;
00168
00169 void rkstep (const RKData::Data & sdata, RKData::Data & ddata) const;
00170 void rkck (const RKData::Data & sdata, RKData::Data & ddata, std::vector<double> & errors) const;
00171 };
00172
00173
00174 }
00175
00176 #endif