CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

RKIntegrator.cc
Go to the documentation of this file.
1// -*- C++ -*-
2// $Id:
6#include <climits>
7#include <stdexcept>
8namespace Genfun {
9FUNCTION_OBJECT_IMP(RKIntegrator::RKFunction)
10
11RKIntegrator::RKFunction::RKFunction(RKData *data, unsigned int index)
12 :_data(data),
13 _index(index)
14{
15 _data->ref();
16}
17
22
24 :AbsFunction(right),
25 _data(right._data),
26 _index(right._index)
27{
28 _data->ref();
29}
30
31
33 if (t<0) return 0;
34 if (!_data->_locked) _data->lock();
35
36 // Do this first, thereafter, just read the cache
37 _data->recache();
38
39
40 // If the cache is empty, make an entry for t=0;
41 size_t nvar = _data->_startingValParameter.size();
42 if (_data->_fx.empty()) {
43 RKData::Data d(nvar);
44 d.time=0;
45 Argument arg(nvar);
46 for (size_t f=0;f<nvar;f++) {
48 arg[f]=d.variable[f];
49 }
50 _data->_fx.insert(d);
51 }
52
53 if (t==0) return (*_data->_fx.begin()).variable[_index];
54
55 RKData::Data dt(nvar);
56 dt.time=t;
57 std::set<RKData::Data>::iterator l =_data->_fx.lower_bound(dt);
58
59 // We did find an exact value (measure 0), just return it.
60 if (l!=_data->_fx.end() && (*l).time==t) {
61 return (*l).variable[_index];
62 }
63
64 else {
65 std::set<RKData::Data>::iterator u =_data->_fx.upper_bound(dt);
66
67 while (u==_data->_fx.end()) {
68 u--;
69 RKData::Data newData(nvar);;
70 _data->_stepper->step(_data,*u, newData, 0);
71 _data->_fx.insert(l,newData);
72 if (newData.time==t) return newData.variable[_index];
73 u = _data->_fx.upper_bound(dt);
74 }
75 u--;
76 _data->_stepper->step(_data,*u, dt, t);
77 return dt.variable[_index];
78 }
79}
80
81
83}
84
85RKIntegrator::RKData::~RKData() {
86 for (size_t i=0;i<_startingValParameter.size();i++) delete _startingValParameter[i];
87 for (size_t i=0;i<_controlParameter.size();i++) delete _controlParameter[i];
88 for (size_t i=0;i<_diffEqn.size(); i++) delete _diffEqn[i];
89 delete _stepper;
90}
91
93 _data(new RKData())
94{
95 if (stepper) _data->_stepper=stepper->clone();
96 else _data->_stepper= new AdaptiveRKStepper();
97 _data->ref();
98}
99
101 _data->unref();
102 for (size_t i=0;i<_fcn.size();i++) delete _fcn[i];
103}
104
106 const std::string &variableName,
107 double defStartingValue,
108 double defValueMin,
109 double defValueMax) {
110 Parameter *par = new Parameter(variableName, defStartingValue, defValueMin, defValueMax);
111 _data->_startingValParameter.push_back(par);
112 _data->_diffEqn.push_back(diffEquation->clone());
113 _data->_startingValParameterCache.push_back(defStartingValue);
114 _fcn.push_back(new RKFunction(_data,_fcn.size()));
115 return par;
116}
117
118
119
120
121
122Parameter * RKIntegrator::createControlParameter (const std::string & variableName,
123 double defStartingValue,
124 double startingValueMin,
125 double startingValueMax) {
126
127 Parameter *par = new Parameter(variableName, defStartingValue, startingValueMin, startingValueMax);
128 _data->_controlParameter.push_back(par);
129 _data->_controlParameterCache.push_back(defStartingValue);
130 return par;
131
132}
133
134
135
137 return _fcn[i];
138}
139
140
141
143 if (!_locked) {
144 unsigned int size = _diffEqn.size();
145 for (size_t i=0;i<size;i++) {
146 if (!(_diffEqn[i]->dimensionality()==size)) throw std::runtime_error("Runtime error in RKIntegrator");
147 }
148 _locked=true;
149 }
150}
151
152
153
155
156 bool stale=false;
157 if (!stale) {
158 for (size_t p=0;p<_startingValParameter.size();p++) {
159 if (_startingValParameter[p]->getValue()!=_startingValParameterCache[p]) {
160 _startingValParameterCache[p]=_startingValParameter[p]->getValue();
161 stale=true;
162 break;
163 }
164 }
165 }
166
167 if (!stale) {
168 for (size_t p=0;p<_controlParameter.size();p++) {
169 if (_controlParameter[p]->getValue()!=_controlParameterCache[p]) {
170 _controlParameterCache[p]=_controlParameter[p]->getValue();
171 stale=true;
172 break;
173 }
174 }
175 }
176
177 if (stale) {
178 _fx.erase(_fx.begin(),_fx.end());
179 }
180
181}
182
183
184
186
187
188} // namespace Genfun
#define FUNCTION_OBJECT_IMP(classname)
virtual AbsFunction * clone() const =0
void unref() const
Definition RCBase.cc:20
void ref() const
Definition RCBase.cc:15
std::vector< const AbsFunction * > _diffEqn
RKFunction(RKData *data, unsigned int index)
virtual double operator()(double argument) const
virtual void step(const RKIntegrator::RKData *data, const RKIntegrator::RKData::Data &sdata, RKIntegrator::RKData::Data &ddata, double timeLimit=0) const =0
virtual RKStepper * clone() const =0
Parameter * addDiffEquation(const AbsFunction *diffEquation, const std::string &variableName="anon", double defStartingValue=0.0, double startingValueMin=0.0, double startingValueMax=0.0)
Parameter * createControlParameter(const std::string &variableName="anon", double defStartingValue=0.0, double startingValueMin=0.0, double startingValueMax=0.0)
RKIntegrator(const RKStepper *stepper=NULL)
const RKFunction * getFunction(unsigned int i) const
void f(void g())