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

CLHEP/GenericFunctions/RKIntegrator.hh
Go to the documentation of this file.
1 // -*- C++ -*-
2 // $Id:
3 //---------------------Runge-Kutte Integrator-------------------------------//
4 // //
5 // Class RKIntegrator //
6 // Joe Boudreau, November 2002 //
7 // //
8 // This is a Runge-Kutta Numerical integrator for a set of N autonomous //
9 // first order differential equations in N variables. The point is to //
10 // create one or more functions which are defined by A) the differential //
11 // equations governing their time evolution, and B) their values at time //
12 // t=0. //
13 // //
14 // You add differential eqns one at a time to this integrator. Each one //
15 // is a GENFUNCTION governing the time evolution of the i^th variable, and //
16 // should depend on all of the N variables, but not on the time //
17 // explicitly. You should add N differential equations in all. Each //
18 // time you add a differential equation the integrator creates a parameter //
19 // for you representing the starting value of the variable, and returns a //
20 // pointer. You may either set the values of that parameter to desired //
21 // values or else connect it to an external parameter if you wish to vary //
22 // the shape of the function by adjusting starting values. //
23 // //
24 // In addition, you may request the integrator to create a control //
25 // parameter. The control parameter may also be set, or connected. //
26 // It can be used in the equations that define the time evolution of the //
27 // variables. //
28 //--------------------------------------------------------------------------//
29 #ifndef RKIntegrator_h
30 #define RKIntegrator_h 1
34 #include <vector>
35 #include <set>
36 namespace Genfun {
37 
43  class RKIntegrator {
44 
45  public:
46 
47  // Some helper classes:
48  class RKFunction;
49  class RKData;
50 
51  // Constructor
52  RKIntegrator();
53 
54  // Destructor
55  virtual ~RKIntegrator();
56 
57  // Add a differential equation governing the time evolution of the next variable.
58  // Get back a parameter representing the starting value of that variable. You
59  // can either arrange for that parameter to have the right starting value, or you
60  // can connect it to another parameter so that you may change it.
61  Parameter * addDiffEquation (const AbsFunction * diffEquation,
62  const std::string & variableName="anon",
63  double defStartingValue=0.0,
64  double startingValueMin=0.0,
65  double startingValueMax=0.0);
66 
67 
68  // Create a control parameter. You can then connnect this to some other
69  // parameter.
70  Parameter *createControlParameter (const std::string & variableName="anon",
71  double defStartingValue=0.0,
72  double startingValueMin=0.0,
73  double startingValueMax=0.0);
74 
75  // Get back a function. This function will now actually change as parameters
76  // are changed; this includes both control parameters and starting value
77  // parameters.
78  const RKFunction *getFunction(unsigned int i) const;
79 
80 
81  private:
82 
83  // It is illegal to assign an adjustable constant
84  const RKIntegrator & operator=(const RKIntegrator &right);
85 
86  // It is illegal to copy an RKIntegrator
87  RKIntegrator(const RKIntegrator &right);
88 
89  // Here is the data, it belongs to the integrator and to the
90  // functions, and is reference counted:
91  RKData *_data;
92 
93 
94  // Here are the functions:
95  std::vector<const RKFunction *> _fcn;
96 
97 
98  };
99 
100 
102 
103 
104  public:
105 
106  struct Data{
107 
108  std::vector<double> variable;
109  mutable std::vector<double> firstDerivative;
110  double time;
111  mutable bool dcalc;
112 
113  Data(int size): variable(size), firstDerivative(size), time(0), dcalc(false) {}
114  bool operator < (const Data & right) const { return time < right.time; }
115  bool operator == (const Data & right) const { return time==right.time; }
116  };
117 
118  RKData();
119  void lock();
120  void recache();
121 
122  std::vector<Parameter *> _startingValParameter;
123  std::vector<double> _startingValParameterCache;
124 
125  std::vector <Parameter *> _controlParameter;
126  std::vector <double> _controlParameterCache;
127 
128  std::vector<const AbsFunction *> _diffEqn;
129  std::set<Data > _fx;
130  bool _locked;
131 
132  private:
133 
134  ~RKData();
135  friend class ImaginaryFriend; // Silence compiler warnings.
136 
137  };
138 
140 
142 
143  public:
144 
145  // Constructor
146  RKFunction(RKData *data, unsigned int index);
147 
148  // Destructor
149  virtual ~RKFunction();
150 
151  // Copy constructor
152  RKFunction(const RKFunction &right);
153 
154  // Retreive function value
155  virtual double operator ()(double argument) const;
156  virtual double operator ()(const Argument & a) const {return operator() (a[0]);}
157 
158  private:
159 
160  // It is illegal to assign a RKFunction
161  const RKFunction & operator=(const RKFunction &right);
162 
163  // The shared data:
164  RKData *_data;
165  const unsigned int _index;
166 
167  void rkstep (const RKData::Data & sdata, RKData::Data & ddata) const;
168  void rkck (const RKData::Data & sdata, RKData::Data & ddata, std::vector<double> & errors) const;
169 };
170 
171 
172 } // namespace Genfun
173 
174 #endif