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>
8 namespace Genfun {
9 FUNCTION_OBJECT_IMP(RKIntegrator::RKFunction)
10 
11 RKIntegrator::RKFunction::RKFunction(RKData *data, unsigned int index)
12  :_data(data),
13  _index(index)
14 {
15  _data->ref();
16 }
17 
19 {
20  _data->unref();
21 }
22 
24  :AbsFunction(right),
25  _data(right._data),
26  _index(right._index)
27 {
28  _data->ref();
29 }
30 
31 
32 double RKIntegrator::RKFunction::operator() (double t) const {
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 
82 RKIntegrator::RKData::RKData():_locked(false) {
83 }
84 
85 RKIntegrator::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 
122 Parameter * 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)
Definition: RKIntegrator.cc:11
virtual double operator()(double argument) const
Definition: RKIntegrator.cc:32
virtual RKStepper * clone() const =0
virtual void step(const RKIntegrator::RKData *data, const RKIntegrator::RKData::Data &sdata, RKIntegrator::RKData::Data &ddata, double timeLimit=0) 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)
Definition: RKIntegrator.cc:92
const RKFunction * getFunction(unsigned int i) const
void f(void g())
Definition: excDblThrow.cc:38