OSHUN  beta
Arbitrary Order Spherical-Harmonic 1D-3P Vlasov-Fokker-Planck-Maxwell code
lib-algorithms.h
Go to the documentation of this file.
1 
13 //--------------------------------------------------------------
14 #ifndef ALGORITHM_LIBRARY_H
15 #define ALGORITHM_LIBRARY_H
16 
17 
18 //**************************************************************
19 namespace Algorithms {
20 
21 
22 
23 //--------------------------------------------------------------
24  template<typename T>
25  valarray<T> MakeAxis(const T min, const T max, const size_t N){
26  valarray<T> v(N);
27  for (size_t i(0); i < N; ++i) {
28  v[i] = static_cast<T>(i);
29  }
30  v *= (max-min)/(static_cast<T>(N-1));
31  v += min;
32  return v;
33  }
34 
35  template<typename T>
36  valarray<T> MakeCAxis(const T min, const T max, const size_t N){
37  valarray<T> v(N);
38  for (size_t i(0); i < N; ++i) {
39  v[i] = static_cast<T>(i);
40  }
41  v *= (max-min)/(static_cast<T>(N));
42  v += min+0.5*(max-min)/(static_cast<T>(N));
43  return v;
44  }
45 
46  template<typename T> class Axis {
47  public:
48  valarray<T>* v;
49 
50 // Constructor
51  Axis(const T min=0, const T max=0, const size_t N=1){
52  if (N > 1) {
53  v = new valarray<T>(MakeAxis(min,max,N));
54  }
55  else v = NULL;
56  }
57 // Copy constructor
58 
59  Axis(const Axis& other){
60  if (other.v != NULL) {
61  v = new valarray<T>(*(other.v));
62  }
63  else v = NULL;
64  }
65 
66 // Destructor
67  ~Axis() { delete v; }
68  };
69 
70  template<typename T> class CAxis {
71  public:
72  valarray<T>* v;
73 
74 // Constructor
75  CAxis(const T min=0, const T max=0, const size_t N=1){
76  if (N > 1) {
77  v = new valarray<T>(MakeCAxis(min,max,N));
78  }
79  else v = NULL;
80  }
81 // Copy constructor
82 
83  CAxis(const CAxis& other){
84  if (other.v != NULL) {
85  v = new valarray<T>(*(other.v));
86  }
87  else v = NULL;
88  }
89 
90 // Destructor
91  ~CAxis() { delete v; }
92  };
93 
94  template<typename T> class AxisBundle {
95  public:
96 // Local axes x(0) --> x, x(1) --> y, x(2) --> z
97  valarray<T> x(size_t i) const { return *(_x[i].v); }
98  size_t Nx(size_t i) const { return (*(_x[i].v)).size(); }
99  T xmin(size_t i) const { return (*(_x[i].v))[0]-0.5*_dx[i]; }
100  T xmax(size_t i) const { return (*(_x[i].v))[Nx(i)-1]+0.5*_dx[i]; }
101  size_t xdim() const { return _x.size(); }
102  T dx(size_t i) const { return _dx[i]; }
103 
104 // Global axes xg(0) --> x, xg(1) --> y, xg(2) --> z
105  valarray<T> xg(size_t i) const { return *(_xg[i].v); }
106  size_t Nxg(size_t i) const { return (*(_xg[i].v)).size(); }
107  T xgmin(size_t i) const { return (*(_xg[i].v))[0]-0.5*_dx[i]; }
108  T xgmax(size_t i) const { return (*(_xg[i].v))[Nxg(i)-1]+0.5*_dx[i]; }
109  size_t xgdim() const { return _xg.size(); }
110 
111 // Global axes p(0) --> species1, p(1) --> species2 ...
112  valarray<T> p(size_t i) const { return *(_p[i].v); }
113  size_t Np(size_t i) const { return (*(_p[i].v)).size(); }
114  T pmin(size_t i) const { return (*(_p[i].v))[0]-0.5*_dp[i]; }
115  T pmax(size_t i) const { return (*(_p[i].v))[Np(i)-1]+0.5*_dp[i]; }
116  size_t pdim() const { return _p.size(); }
117  T dp(size_t i) const { return _dp[i]; }
118 
119 // Global axes px(0) --> species1, px(1) --> species2 ...
120  valarray<T> px(size_t i) const { return *(_px[i].v); }
121  size_t Npx(size_t i) const { return (*(_px[i].v)).size(); }
122  T pxmin(size_t i) const { return (*(_px[i].v))[0]; }
123  T pxmax(size_t i) const { return (*(_px[i].v))[Npx(i)-1]; }
124  size_t pxdim() const { return _px.size(); }
125 
126  // Global axes py(0) --> species1, py(1) --> species2 ...
127  valarray<T> py(size_t i) const { return *(_py[i].v); }
128  size_t Npy(size_t i) const { return (*(_py[i].v)).size(); }
129  T pymin(size_t i) const { return (*(_py[i].v))[0]; }
130  T pymax(size_t i) const { return (*(_py[i].v))[Npy(i)-1]; }
131  size_t pydim() const { return _py.size(); }
132 
133  // Global axes py(0) --> species1, py(1) --> species2 ...
134  valarray<T> pz(size_t i) const { return *(_pz[i].v); }
135  size_t Npz(size_t i) const { return (*(_pz[i].v)).size(); }
136  T pzmin(size_t i) const { return (*(_pz[i].v))[0]; }
137  T pzmax(size_t i) const { return (*(_pz[i].v))[Npz(i)-1]; }
138  size_t pzdim() const { return _pz.size(); }
139 
140 // Constructor
141  AxisBundle(const vector<T> _xmin, const vector<T> _xmax, const vector<size_t> _Nx, // local spatial axes
142  const vector<T> _xgmin, const vector<T> _xgmax, const vector<size_t> _Nxg, // global spatial axes
143  const vector<T> _pmax, const vector<size_t> _Np, // momentum axis for each species
144  const vector<size_t> _Npx, const vector<size_t> _Npy, const vector<size_t> _Npz ){ // dimensions in configuration space
145  for (size_t i(0); i < _Nx.size(); ++i) {
146  _x.push_back( CAxis<T>( _xmin[i], _xmax[i], _Nx[i]));
147  }
148  for (size_t i(0); i < _Nxg.size(); ++i) {
149  _xg.push_back( CAxis<T>( _xgmin[i], _xgmax[i], _Nxg[i]));
150  }
151  for (size_t i(0); i < _Nxg.size(); ++i) {
152  _dx.push_back((_xmax[i]-_xmin[i])/(static_cast<T>(_Nx[i])));
153  }
154  // for (size_t i(0); i < _Np.size(); ++i) {
155  // _p.push_back( Axis<T>( _pmax[i]/(static_cast<T>(2 * _Np[i])), _pmax[i], _Np[i]) );
156  // }
157  for (size_t i(0); i < _Np.size(); ++i) {
158  _p.push_back( CAxis<T>(static_cast<T>(0.0), _pmax[i], _Np[i]) );
159  }
160 
161  for (size_t i(0); i < _Np.size(); ++i) {
162  _dp.push_back((_pmax[i])/(static_cast<T>(_Np[i])));
163  }
164 
165  for (size_t i(0); i < _Npx.size(); ++i) {
166  _px.push_back( Axis<T>( static_cast<T>(-1.0) * _pmax[i], _pmax[i], _Npx[i]));
167  }
168  for (size_t i(0); i < _Npy.size(); ++i) {
169  _py.push_back( Axis<T>( static_cast<T>(-1.0) * _pmax[i], _pmax[i], _Npy[i]));
170  }
171  for (size_t i(0); i < _Npz.size(); ++i) {
172  _pz.push_back( Axis<T>( static_cast<T>(-1.0) * _pmax[i], _pmax[i], _Npz[i]));
173  }
174  }
175 
176  AxisBundle(const AxisBundle& a){
177  for (size_t i(0); i < a.xdim(); ++i) {
178  _x.push_back( CAxis<T>(a.xmin(i), a.xmax(i), a.Nx(i)) );
179  }
180  for (size_t i(0); i < a.xgdim(); ++i) {
181  _xg.push_back( CAxis<T>(a.xgmin(i), a.xgmax(i), a.Nxg(i)) );
182  }
183  for (size_t i(0); i < a.xdim(); ++i) {
184  _dx.push_back( (a.xgmax(i)-a.xgmin(i))/a.Nxg(i) );
185  }
186  // for (size_t i(0); i < a.pdim(); ++i) {
187  // _p.push_back( Axis<T>(a.pmax(i)/static_cast<T>(2 * a.Np(i)),a.pmax(i),a.Np(i)) );
188  // }
189  for (size_t i(0); i < a.pdim(); ++i) {
190  _p.push_back( CAxis<T>(a.pmin(i),a.pmax(i),a.Np(i)) );
191  }
192  for (size_t i(0); i < a.pdim(); ++i) {
193  _dp.push_back( a.pmax(i)/a.Np(i) );
194  }
195 
196  for (size_t i(0); i < a.pxdim(); ++i) {
197  _px.push_back( Axis<T>( a.pxmin(i), a.pxmin(i), a.Npx(i)) );
198  }
199  for (size_t i(0); i < a.pydim(); ++i) {
200  _py.push_back( Axis<T>( a.pymin(i), a.pymin(i), a.Npy(i)) );
201  }
202  for (size_t i(0); i < a.pzdim(); ++i) {
203  _pz.push_back( Axis<T>( a.pzmin(i), a.pzmin(i), a.Npz(i)) );
204  }
205  }
206 
207  ~AxisBundle(){ }
208 
209  private:
210  // vector< Axis<T> > _x, _xg, _p, _px;
211  //
212  vector< Axis<T> > _px, _py, _pz;
213  vector< CAxis<T> > _x, _xg, _p;
214  vector<T> _dx, _dp;
215 
216 
217 
218 
219  };
220 //--------------------------------------------------------------
221 
222 //--------------------------------------------------------------
223 // LEGENDRE POLYNOMIALS
224 // Calculate the legendre polynomials using the recurrance relations
225 // For m0 = 0
226  template<class T>
227  valarray<T> Legendre(const T x, const size_t Nl){ // where Nl is to denote l0+1
228 
229  valarray<T> P_Legendre(0.0, Nl);
230  if ( abs(x) > 1 ) return P_Legendre;
231 
232  T r1, sqrtx = sqrt(1.0-x*x);
233 
234 // Initialization
235  P_Legendre[0] = 1.0;
236  P_Legendre[1] = x;
237 
238  for (size_t l(1); l < Nl - 1; ++l){
239  r1 = 1.0 / double(l + 1);
240  P_Legendre[l+1] = P_Legendre[l] * (x*(2.0*l+1.0) * r1) -
241  P_Legendre[l-1]*(double(l) * r1);
242  }
243 
244  return P_Legendre;
245  }
246 
247 // For m0 > 0
248  template<class T>
249  Array2D<T> Legendre(const T x, const size_t Nl, const size_t Nm){
250 
251 // Local variables
252  if (Nl < Nm) {
253  cout << "ERROR: " << "l0+1 = " << Nl << " < " << Nm << " = m0+1\n";
254  exit(1);
255  }
256 
257  Array2D<T> P_Legendre(Nl+1,Nm+1);
258 
259  if ( abs(x) > 1 ) {
260 // cout << "ERROR: " << "|" << x << "| > 1 is not a valid cosine\n";
261 // cout << "setting values to 0.0";
262  P_Legendre = 0.0;
263 
264  }
265  else
266  {
267 
268  T r1, sqrtx = sqrt(1.0-x*x), fact = 1.0;
269 
270  // Initialization
271  P_Legendre = 0.0;
272  P_Legendre(0,0) = 1.0;
273 
274  for (size_t l = 1; l < Nm+1; ++l){
275  // std::cout << "\n l=" << l;
276  P_Legendre(l,l) = - P_Legendre(l-1,l-1)*(fact*sqrtx);
277  fact += 2.0;
278  }
279  // std::cout << "\n ";
280 
281  for (size_t l = 0; l < ((Nm < Nl) ? Nm+1 : Nl+1); ++l){
282  // std::cout << "\n l=" << l;
283  P_Legendre(l+1,l) = P_Legendre(l,l)*(x*(2.0*l+1.0));
284  }
285  // std::cout << "\n ";
286 
287  for (size_t m = 0; m < Nm+1; ++m){
288  for (size_t l = m+1; l < Nl ; ++l){
289  // std::cout << "\n l=" << l << ", m=" << m;
290  r1 = 1.0 / double(l - m + 1);
291  P_Legendre(l+1,m) = P_Legendre(l,m) * (x*(2.0*l+1.0) * r1) -
292  P_Legendre(l-1,m)*(double(l+m) * r1);
293  }
294  }
295  }
296 
297  return P_Legendre;
298  }
299 //--------------------------------------------------------------
300 
301 //--------------------------------------------------------------
302 // MOMENTS
303 // p-th moment of a quantity q(x)
304  template<class T>
305  T moment(const vector<T> q, const vector<T> x, const int p){
306 
307  T integral(0.0);
308 // TODO: Integral values up to the zeroth cell and above the last cell
309 // integral += q[0] * pow(x[0], p) // += Q_0*x_0^p * (x_1-x_0)
310 // * (x[1]-x[0]);
311 // for (size_t i(1); i < q.size()-1; ++i){ // += Q_i*x_i^p * (x_{i+1}-x_{i-1})
312 // integral += q[i] * pow(x[i], p)
313 // * (x[i+1] - x[i-1]);
314 // }
315 // integral += q[q.size()-1] * pow(x[q.size()-1], p) // += Q_n*x_n^p * (x_{n}-x_{n-1})
316 // * (x[q.size()-1]-x[q.size()-2]);
317 // return integral*0.5;
318 
319  integral += q[0] * pow(x[0], p) // += Q_0*x_0^p * (x_1-x_0)
320  * (x[1]-x[0]);
321  for (size_t i(1); i < q.size()-1; ++i){ // += Q_i*x_i^p * (x_{i+1}-x_{i-1})
322  integral += q[i] * pow(x[i], p)
323  * 0.5*(x[i+1] - x[i-1]);
324  }
325  integral += q[q.size()-1] * pow(x[q.size()-1], p) // += Q_n*x_n^p * (x_{n}-x_{n-1})
326  * (x[q.size()-1]-x[q.size()-2]);
327  return integral;
328  }
329 
330  template<class T>
331  T moment(const vector<T> q, const valarray<T> x, const int p){
332 
333  T integral(0.0);
334 // TODO: Integral values up to the zeroth cell and above the last cell
335 // integral += q[0] * pow(x[0], p) // += Q_0*x_0^p * (x_1-x_0)
336 // * (x[1]-x[0]);
337 // for (size_t i(1); i < q.size()-1; ++i){ // += Q_i*x_i^p * (x_{i+1}-x_{i-1})
338 // integral += q[i] * pow(x[i], p)
339 // * (x[i+1] - x[i-1]);
340 // }
341 // integral += q[q.size()-1] * pow(x[q.size()-1], p) // += Q_n*x_n^p * (x_{n}-x_{n-1})
342 // * (x[q.size()-1]-x[q.size()-2]);
343 // return integral*0.5;
344  integral += q[0] * pow(x[0], p) // += Q_0*x_0^p * (x_1-x_0)
345  * (x[1]-x[0]);
346  for (size_t i(1); i < q.size()-1; ++i){ // += Q_i*x_i^p * (x_{i+1}-x_{i-1})
347  integral += q[i] * pow(x[i], p)
348  * 0.5*(x[i+1] - x[i-1]);
349  }
350  integral += q[q.size()-1] * pow(x[q.size()-1], p) // += Q_n*x_n^p * (x_{n}-x_{n-1})
351  * (x[q.size()-1]-x[q.size()-2]);
352  return integral;
353  }
354  template<class T>
355  T moment(const valarray<T> q, const valarray<T> x, const int p){
356 
357  T integral(0.0);
358 // TODO: Integral values up to the zeroth cell and above the last cell
359 // integral += q[0] * pow(x[0], p) // += Q_0*x_0^p * (x_1-x_0)
360 // * (x[1]-x[0]);
361 // for (size_t i(1); i < q.size()-1; ++i){ // += Q_i*x_i^p * (x_{i+1}-x_{i-1})
362 // integral += q[i] * pow(x[i], p)
363 // * (x[i+1] - x[i-1]);
364 // }
365 // integral += q[q.size()-1] * pow(x[q.size()-1], p) // += Q_n*x_n^p * (x_{n}-x_{n-1})
366 // * (x[q.size()-1]-x[q.size()-2]);
367 // return integral*0.5;
368  integral += q[0] * pow(x[0], p) // += Q_0*x_0^p * (x_1-x_0)
369  * (x[1]-x[0]);
370  for (size_t i(1); i < q.size()-1; ++i){ // += Q_i*x_i^p * (x_{i+1}-x_{i-1})
371  integral += q[i] * pow(x[i], p)
372  * 0.5*(x[i+1] - x[i-1]);
373  }
374  integral += q[q.size()-1] * pow(x[q.size()-1], p) // += Q_n*x_n^p * (x_{n}-x_{n-1})
375  * (x[q.size()-1]-x[q.size()-2]);
376  return integral;
377  }
378 //--------------------------------------------------------------
379 // relativistic moments for inverse gamma and gamma
380 //--------------------------------------------------------------
381  template<class T>
382  T relativistic_invg_moment(const vector<T> q, const valarray<T> x, const int p){
383 
384  T integral(0.0);
385 // TODO: Integral values up to the zeroth cell and above the last cell
386  // integral += q[0] * pow(x[0], p) // += Q_0*x_0^p * (x_1-x_0)
387  // * (x[1]-x[0])
388  // / sqrt(1.0+x[0]*x[0]);
389  // for (size_t i(1); i < q.size()-1; ++i){ // += Q_i*x_i^p * (x_{i+1}-x_{i-1})
390  // integral += q[i] * pow(x[i], p)
391  // * (x[i+1] - x[i-1])
392  // / sqrt(1.0+x[i]*x[i]);
393  // }
394  // integral += q[q.size()-1] * pow(x[q.size()-1], p) // += Q_n*x_n^p * (x_{n}-x_{n-1})
395  // * (x[q.size()-1]-x[q.size()-2])
396  // / sqrt(1.0+x[q.size()-1]*x[q.size()-1]);
397  // return integral*0.5;
398 
399  integral += q[0] * pow(x[0], p) // += Q_0*x_0^p * (x_1-x_0)
400  * (x[1]-x[0])
401  / sqrt(1.0+x[0]*x[0]);
402  for (size_t i(1); i < q.size()-1; ++i){ // += Q_i*x_i^p * (x_{i+1}-x_{i-1})
403  integral += q[i] * pow(x[i], p)
404  * 0.5*(x[i+1] - x[i-1])
405  / sqrt(1.0+x[i]*x[i]);
406  }
407  integral += q[q.size()-1] * pow(x[q.size()-1], p) // += Q_n*x_n^p * (x_{n}-x_{n-1})
408  * (x[q.size()-1]-x[q.size()-2])
409  / sqrt(1.0+x[q.size()-1]*x[q.size()-1]);
410  return integral;
411  }
412 
413  template<class T>
414  T relativistic_invg_moment(const valarray<T> q, const valarray<T> x, const int p){
415 
416  T integral(0.0);
417 // TODO: Integral values up to the zeroth cell and above the last cell
418  // integral += q[0] * pow(x[0], p) // += Q_0*x_0^p * (x_1-x_0)
419  // * (x[1]-x[0])
420  // / sqrt(1.0+x[0]*x[0]);
421  // for (size_t i(1); i < q.size()-1; ++i){ // += Q_i*x_i^p * (x_{i+1}-x_{i-1})
422  // integral += q[i] * pow(x[i], p)
423  // * (x[i+1] - x[i-1])
424  // / sqrt(1.0+x[i]*x[i]);
425  // }
426  // integral += q[q.size()-1] * pow(x[q.size()-1], p) // += Q_n*x_n^p * (x_{n}-x_{n-1})
427  // * (x[q.size()-1]-x[q.size()-2])
428  // / sqrt(1.0+x[q.size-1]*x[q.size-1]);
429  // return integral*0.5;
430  integral += q[0] * pow(x[0], p) // += Q_0*x_0^p * (x_1-x_0)
431  * (x[1]-x[0])
432  / sqrt(1.0+x[0]*x[0]);
433  for (size_t i(1); i < q.size()-1; ++i){ // += Q_i*x_i^p * (x_{i+1}-x_{i-1})
434  integral += q[i] * pow(x[i], p)
435  * 0.5*(x[i+1] - x[i-1])
436  / sqrt(1.0+x[i]*x[i]);
437  }
438  integral += q[q.size()-1] * pow(x[q.size()-1], p) // += Q_n*x_n^p * (x_{n}-x_{n-1})
439  * (x[q.size()-1]-x[q.size()-2])
440  / sqrt(1.0+x[q.size-1]*x[q.size-1]);
441  return integral;
442  }
443 //--------------------------------------------------------------
444  template<class T>
445  T relativistic_gamma_moment(const valarray<T> q, const valarray<T> x, const int p){
446 
447  T integral(0.0);
448 // TODO: Integral values up to the zeroth cell and above the last cell
449 // integral += q[0] * pow(x[0], p) // += Q_0*x_0^p * (x_1-x_0)
450 // * (x[1]-x[0])
451 // * sqrt(1.0+x[0]*x[0]);
452 // for (size_t i(1); i < q.size()-1; ++i){ // += Q_i*x_i^p * (x_{i+1}-x_{i-1})
453 // integral += q[i] * pow(x[i], p)
454 // * (x[i+1] - x[i-1])
455 // * sqrt(1.0+x[i]*x[i]);
456 // }
457 // integral += q[q.size()-1] * pow(x[q.size()-1], p) // += Q_n*x_n^p * (x_{n}-x_{n-1})
458 // * (x[q.size()-1]-x[q.size()-2])
459 // * sqrt(1.0+x[q.size-1]*x[q.size-1]);
460 // return integral*0.5;
461  integral += q[0] * pow(x[0], p) // += Q_0*x_0^p * (x_1-x_0)
462  * (x[1]-x[0])
463  * sqrt(1.0+x[0]*x[0]);
464  for (size_t i(1); i < q.size()-1; ++i){ // += Q_i*x_i^p * (x_{i+1}-x_{i-1})
465  integral += q[i] * pow(x[i], p)
466  * 0.5*(x[i+1] - x[i-1])
467  * sqrt(1.0+x[i]*x[i]);
468  }
469  integral += q[q.size()-1] * pow(x[q.size()-1], p) // += Q_n*x_n^p * (x_{n}-x_{n-1})
470  * (x[q.size()-1]-x[q.size()-2])
471  * sqrt(1.0+x[q.size-1]*x[q.size-1]);
472  return integral;
473  }
474 //--------------------------------------------------------------
475 
476 
477 //--------------------------------------------------------------
478 // RUNGE-KUTTA METHODS
479 //--------------------------------------------------------------
480  template<typename T> class AbstFunctor {
481 // abstract functor interface
482  public:
483  virtual void operator()(const T& Yin, T& Yslope)=0; // call using operator
484  virtual void operator()(const T& Yin, T& Yslope, size_t dir)=0; // call using operator
485  };
486 
487 
488 // RK2
489  template<class T> class RK2 {
490  public:
491 // Constructor
492  RK2(T& Yin): Y0(Yin), Yh(Yin) { }
493 
494 // Main function
495  T& operator()(T& Y, double h, AbstFunctor<T>* F);
496  T& operator()(T& Y, double h, AbstFunctor<T>* F, size_t dir);
497 
498  private:
499 // R-K copies for the data
500  T Y0, Yh;
501  };
502 
503  template<class T> T& RK2<T>::operator()
504  (T& Y, double h, AbstFunctor<T>* F) {
505 // Take a step using RK2
506 
507 // Initialization
508  Y0 = Y;
509 
510 // Step 1
511  (*F)(Y0,Yh); Yh *= h; // Yh = h*F(Y0)
512  Y0 += Yh; // Y0 = Y0 + h*Yh
513  Yh *= 0.5; Y += Yh; // Y = Y + (h/2)*F(Y0)
514 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
515 
516 // Step 2
517  (*F)(Y0,Yh); Yh *= 0.5*h; // Yh = (h/2)*F(Y0)
518  Y += Yh; // Y = Y + (h/2)*F(Y0)
519 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
520 
521  return Y;
522  }
523  template<class T> T& RK2<T>::operator()
524  (T& Y, double h, AbstFunctor<T>* F, size_t dir) {
525 // Take a step using RK2
526 
527 // Initialization
528  Y0 = Y;
529 
530 // Step 1
531  (*F)(Y0,Yh,dir); Yh *= h; // Yh = h*F(Y0)
532  Y0 += Yh; // Y0 = Y0 + h*Yh
533  Yh *= 0.5; Y += Yh; // Y = Y + (h/2)*F(Y0)
534 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
535 
536 // Step 2
537  (*F)(Y0,Yh,dir); Yh *= 0.5*h; // Yh = (h/2)*F(Y0)
538  Y += Yh; // Y = Y + (h/2)*F(Y0)
539 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
540 
541  return Y;
542  }
543 //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
544 
545 // RK3 (Osher&Shu Method)
546  template<class T> class RK3 {
547  public:
548 // Constructor
549  RK3(T& Yin): Y0(Yin), Yh(Yin) { }
550 
551 // Main function
552  T& operator()(T& Y, double h, AbstFunctor<T>* F);
553  T& operator()(T& Y, double h, AbstFunctor<T>* F, size_t dir);
554 
555 
556  private:
557 // R-K copies for the data
558  T Y0, Yh;
559  };
560 
561  template<class T> T& RK3<T>::operator()
562  (T& Y, double h, AbstFunctor<T>* F) {
563 // Take a step using RK3
564 
565 // Initialization
566  Y0 = Y;
567 
568 // Step 1
569  (*F)(Y0,Yh); Yh *= h; // Yh = h*F(Y0)
570  Y0 += Yh;
571 // Y0 = Y0 + h*Yh
572 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
573 
574 // Step 2
575  (*F)(Y0,Yh); Yh *= h; // Yh = h*F(Y0)
576  Y0 += Yh; Y *= 3.0; Y0 += Y; Y0 *= 0.25; // Changed Y to 3*Y!
577 // Y0 = 1/4 * ( 3*Y + (Y0 + h*Yh) )
578 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
579 
580 // Step 3
581  (*F)(Y0,Yh); Yh *= h; // Yh = h*F(Y0)
582  Y0 += Yh; Y0 *= 6.0; Y += Y0; Y *= (1.0/9.0);
583 // Y = 1/3 * ( Y + 2 * (Y0+h*Yh) )
584 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
585 
586  return Y;
587  }
588 //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
589  template<class T> T& RK3<T>::operator()
590  (T& Y, double h, AbstFunctor<T>* F, size_t dir) {
591 // Take a step using RK3
592 
593 // Initialization
594  Y0 = Y;
595 
596 // Step 1
597  (*F)(Y0,Yh,dir); Yh *= h; // Yh = h*F(Y0)
598  Y0 += Yh;
599 // Y0 = Y0 + h*Yh
600 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
601 
602 // Step 2
603  (*F)(Y0,Yh,dir); Yh *= h; // Yh = h*F(Y0)
604  Y0 += Yh; Y *= 3.0; Y0 += Y; Y0 *= 0.25; // Changed Y to 3*Y!
605 // Y0 = 1/4 * ( 3*Y + (Y0 + h*Yh) )
606 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
607 
608 // Step 3
609  (*F)(Y0,Yh,dir); Yh *= h; // Yh = h*F(Y0)
610  Y0 += Yh; Y0 *= 6.0; Y += Y0; Y *= (1.0/9.0);
611 // Y = 1/3 * ( Y + 2 * (Y0+h*Yh) )
612 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
613 
614  return Y;
615  }
616 //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
617 
618 
619 // RK4
620  template<class T> class RK4 {
621  public:
622 // Constructor
623  RK4(T& Yin): Y0(Yin), Y1(Yin), Yh(Yin) { }
624 
625 // Main function
626  T& operator()(T& Y, double h, AbstFunctor<T>* F);
627  T& operator()(T& Y, double h, AbstFunctor<T>* F, size_t dir);
628 
629  private:
630 // R-K copies for the data
631  T Y0, Y1, Yh;
632  };
633 
634  template<class T> T& RK4<T>::operator()
635  (T& Y, double h, AbstFunctor<T>* F) {
636 // Take a step using RK4
637 
638 // Initialization
639  Y0 = Y; Y1 = Y;
640 
641 // Step 1
642  (*F)(Y1,Yh); // slope in the beginning
643  Yh *= (0.5*h); Y1 += Yh; // Y1 = Y1 + (h/2)*Yh
644  Yh *= (1.0/3.0); Y += Yh; // Y = Y + (h/6)*Yh
645 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
646 
647 // Step 2
648  (*F)(Y1,Yh); Y1 = Y0; // slope in the middle
649  Yh *= (0.5*h); Y1 += Yh; // Y1 = Y0 + (h/2)*Yh
650  Yh *= (2.0/3.0); Y += Yh; // Y = Y + (h/3)*Yh
651 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
652 
653 // Step 3
654  (*F)(Y1,Yh); // slope in the middle again
655  Yh *= h; Y0 += Yh; // Y0 = Y0 + h*Yh
656  Yh *= (1.0/3.0); Y += Yh; // Y = Y + (h/3)*Yh
657 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
658 
659 // Step 4
660  (*F)(Y0,Yh); // slope at the end
661  Yh *= (h/6.0); Y += Yh; // Y = Y + (h/6)*Yh
662 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
663 
664  return Y;
665  }
666  template<class T> T& RK4<T>::operator()
667  (T& Y, double h, AbstFunctor<T>* F, size_t dir) {
668 // Take a step using RK4
669 
670 // Initialization
671  Y0 = Y; Y1 = Y;
672 
673 // Step 1
674  (*F)(Y1,Yh,dir); // slope in the beginning
675  Yh *= (0.5*h); Y1 += Yh; // Y1 = Y1 + (h/2)*Yh
676  Yh *= (1.0/3.0); Y += Yh; // Y = Y + (h/6)*Yh
677 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
678 
679 // Step 2
680  (*F)(Y1,Yh,dir); Y1 = Y0; // slope in the middle
681  Yh *= (0.5*h); Y1 += Yh; // Y1 = Y0 + (h/2)*Yh
682  Yh *= (2.0/3.0); Y += Yh; // Y = Y + (h/3)*Yh
683 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
684 
685 // Step 3
686  (*F)(Y1,Yh,dir); // slope in the middle again
687  Yh *= h; Y0 += Yh; // Y0 = Y0 + h*Yh
688  Yh *= (1.0/3.0); Y += Yh; // Y = Y + (h/3)*Yh
689 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
690 
691 // Step 4
692  (*F)(Y0,Yh,dir); // slope at the end
693  Yh *= (h/6.0); Y += Yh; // Y = Y + (h/6)*Yh
694 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
695 
696  return Y;
697  }
698 
699 
700 //--------------------------------------------------------------
701 
702  // Leapfrog space (Position verlet)
703  template<class T> class LEAPs {
704  public:
705 // Constructor
706  LEAPs(T& Yin): Y0(Yin), Yh(Yin) { }
707 
708 // Main function
709  T& operator()(T& Y, double h, AbstFunctor<T>* F_space, AbstFunctor<T>* F_momentum);
710  T& operator()(T& Y, double h, AbstFunctor<T>* F_space, AbstFunctor<T>* F_momentum, size_t dir);
711 
712  private:
713 // R-K copies for the data
714  T Y0, Yh;
715  };
716 
717  template<class T> T& LEAPs<T>::operator()
718  (T& Y, double h, AbstFunctor<T>* F_space, AbstFunctor<T>* F_momentum) {
719 // Take a step using LEAPspace
720 
721 // Initialization
722  Y0 = Y;
723 
724  (*F_space)(Y0,Yh); Yh *= 0.5*h; // x = h * F_space(Y0)
725  Y0 += Yh;
726  (*F_momentum)(Y0,Yh); Yh *= h; // p = h * F_momentum(Y0)
727  Y0 += Yh; // Y0 = Y0 + h*Yh
728  (*F_space)(Y0,Yh); Yh *= 0.5*h; // x = h * F_space(Y0)
729  Y0 += Yh;
730 
731  Y = Y0;
732 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
733 
734  return Y;
735  }
736  template<class T> T& LEAPs<T>::operator()
737  (T& Y, double h, AbstFunctor<T>* F_space, AbstFunctor<T>* F_momentum, size_t dir) {
738 // Take a step using LEAPspace
739 
740 // Initialization
741  Y0 = Y;
742 
743 // First step momentum
744  (*F_space)(Y0,Yh); Yh *= 0.5*h; // x = h * F_space(Y0)
745  Y0 += Yh;
746  (*F_momentum)(Y0,Yh); Yh *= h; // p = h * F_momentum(Y0)
747  Y0 += Yh; // Y0 = Y0 + h*Yh
748  (*F_space)(Y0,Yh); Yh *= 0.5*h; // x = h * F_space(Y0)
749  Y0 += Yh;
750 
751  Y = Y0;
752 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
753 
754  return Y;
755  }
756 //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
757 //--------------------------------------------------------------
758 
759  // Leapfrog space (velocity verlet)
760  template<class T> class LEAPv {
761  public:
762 // Constructor
763  LEAPv(T& Yin): Y0(Yin), Yh(Yin) { }
764 
765 // Main function
766  T& operator()(T& Y, double h, AbstFunctor<T>* F_space, AbstFunctor<T>* F_momentum);
767  T& operator()(T& Y, double h, AbstFunctor<T>* F_space, AbstFunctor<T>* F_momentum, size_t dir);
768 
769  private:
770 // R-K copies for the data
771  T Y0, Yh;
772  };
773 
774  template<class T> T& LEAPv<T>::operator()
775  (T& Y, double h, AbstFunctor<T>* F_space, AbstFunctor<T>* F_momentum) {
776 // Take a step using LEAPmomentum
777 
778 // Initialization
779  Y0 = Y;
780 
781  (*F_momentum)(Y0,Yh); Yh *= 0.5*h; // x = h * F_space(Y0)
782  Y0 += Yh;
783  (*F_space)(Y0,Yh); Yh *= h; // p = h * F_momentum(Y0)
784  Y0 += Yh; // Y0 = Y0 + h*Yh
785  (*F_momentum)(Y0,Yh); Yh *= 0.5*h; // x = h * F_space(Y0)
786  Y0 += Yh;
787 
788  Y = Y0;
789 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
790 
791  return Y;
792  }
793  template<class T> T& LEAPv<T>::operator()
794  (T& Y, double h, AbstFunctor<T>* F_space, AbstFunctor<T>* F_momentum, size_t dir) {
795 // Take a step using LEAPmomentum
796 
797 // Initialization
798  Y0 = Y;
799 
800  (*F_momentum)(Y0,Yh); Yh *= 0.5*h; // x = h * F_space(Y0)
801  Y0 += Yh;
802  (*F_space)(Y0,Yh); Yh *= h; // p = h * F_momentum(Y0)
803  Y0 += Yh; // Y0 = Y0 + h*Yh
804  (*F_momentum)(Y0,Yh); Yh *= 0.5*h; // x = h * F_space(Y0)
805  Y0 += Yh;
806 
807  Y = Y0;
808 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
809 
810 
811  return Y;
812  }
813 //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
814 
815  // Position Extended Forest-Ruth Like
816  template<class T> class PEFRL {
817  public:
818 // Constructor
819  PEFRL(T& Yin): Y0(Yin), Yh(Yin),
820  xsi(0.1786178958448091),
821  lambda(-0.2123418310626054),
822  chi(-0.06626458266981849)
823  {}
824 
825 // Main function
826  T& operator()(T& Y, double h, AbstFunctor<T>* F_space, AbstFunctor<T>* F_momentum);
827  T& operator()(T& Y, double h, AbstFunctor<T>* F_space, AbstFunctor<T>* F_momentum, size_t dir);
828 
829  private:
830 // R-K copies for the data
831  T Y0, Yh;
832  double xsi, lambda, chi;
833  };
834 
835  template<class T> T& PEFRL<T>::operator()
836  (T& Y, double h, AbstFunctor<T>* F_space, AbstFunctor<T>* F_momentum) {
837 // Take a step using PEFRL
838 
839 // Initialization
840  Y0 = Y;
841 
842 // First step space
843  (*F_space)(Y0,Yh); Yh *= xsi*h; // x = h * F_space(Y0)
844  Y0 += Yh;
845 
846  (*F_momentum)(Y0,Yh); Yh *= (1.0-2.0*lambda)*0.5*h; // p = h * F_momentum(Y0)
847  Y0 += Yh; // Y0 = Y0 + h*Yh
848 
849  (*F_space)(Y0,Yh); Yh *= chi*h; // x = h * F_space(Y0)
850  Y0 += Yh;
851 
852  (*F_momentum)(Y0,Yh); Yh *= lambda*h; // p = h * F_momentum(Y0)
853  Y0 += Yh; // Y0 = Y0 + h*Yh
854 
855  (*F_space)(Y0,Yh); Yh *= (1.0-2.0*(chi+xsi))*h; // x = h * F_space(Y0)
856  Y0 += Yh;
857 
858  (*F_momentum)(Y0,Yh); Yh *= lambda*h; // p = h * F_momentum(Y0)
859  Y0 += Yh; // Y0 = Y0 + h*Yh
860 
861  (*F_space)(Y0,Yh); Yh *= chi*h; // x = h * F_space(Y0)
862  Y0 += Yh;
863 
864  (*F_momentum)(Y0,Yh); Yh *= (1.0-2.0*lambda)*0.5*h; // p = h * F_momentum(Y0)
865  Y0 += Yh;
866 
867  (*F_space)(Y0,Yh); Yh *= xsi*h; // x = h * F_space(Y0)
868  Y0 += Yh;
869 
870  Y = Y0;
871 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
872 
873  return Y;
874  }
875  template<class T> T& PEFRL<T>::operator()
876  (T& Y, double h, AbstFunctor<T>* F_space, AbstFunctor<T>* F_momentum, size_t dir) {
877 // Take a step using PEFRL
878 
879 // Initialization
880  Y0 = Y;
881 
882 // First step space
883  (*F_space)(Y0,Yh); Yh *= xsi*h; // x = h * F_space(Y0)
884  Y0 += Yh;
885 
886  (*F_momentum)(Y0,Yh); Yh *= (1.0-2.0*lambda)*0.5*h; // p = h * F_momentum(Y0)
887  Y0 += Yh; // Y0 = Y0 + h*Yh
888 
889  (*F_space)(Y0,Yh); Yh *= chi*h; // x = h * F_space(Y0)
890  Y0 += Yh;
891 
892  (*F_momentum)(Y0,Yh); Yh *= lambda*h; // p = h * F_momentum(Y0)
893  Y0 += Yh; // Y0 = Y0 + h*Yh
894 
895  (*F_space)(Y0,Yh); Yh *= (1.0-2.0*(chi+xsi))*h; // x = h * F_space(Y0)
896  Y0 += Yh;
897 
898  (*F_momentum)(Y0,Yh); Yh *= lambda*h; // p = h * F_momentum(Y0)
899  Y0 += Yh; // Y0 = Y0 + h*Yh
900 
901  (*F_space)(Y0,Yh); Yh *= chi*h; // x = h * F_space(Y0)
902  Y0 += Yh;
903 
904  (*F_momentum)(Y0,Yh); Yh *= (1.0-2.0*lambda)*0.5*h; // p = h * F_momentum(Y0)
905  Y0 += Yh;
906 
907  (*F_space)(Y0,Yh); Yh *= xsi*h; // x = h * F_space(Y0)
908  Y0 += Yh;
909 
910  Y = Y0;
911 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
912  return Y;
913  }
914 //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
915 
916 }
917 //**************************************************************
918 
919 
920 
921 //**************************************************************
922 // This Clock controls an iteration loop, given an initial
923 // time tout_start*dt_out it evaluates the number of time-
924 // steps it takes to get to (tout_start+1)*dt_out and can
925 // be incremented;
926 
927 class Clock {
928 public:
929 // Constructor
930  Clock(int tout_start, double dt_out, double CFL) {
931  _hstep = 0;
932  _t_start = double(tout_start)*dt_out;
933  _numh = size_t(static_cast<int>(dt_out/CFL))+1;
934  _h = dt_out/static_cast<double>(_numh);
935  }
936 
937 // Clock readings
938  double h() const {return _h;} // Step size
939  size_t numh() const {return _numh;} // # steps
940  size_t tick() const {return _hstep;} // Current step
941  double time() const {return tick()*h() + _t_start;} // Current time
942 
943 // Increment time
944  Clock& operator++() { ++_hstep; return *this;}
945 
946 private:
947  double _t_start; // Initial output timestep
948  size_t _numh; // # of steps derived from this CFL
949  double _h; // resulting time-step from CFL
950  size_t _hstep;
951 };
952 //--------------------------------------------------------------
953 //**************************************************************
954 
955 
956 #endif
Clock::_h
double _h
Definition: lib-algorithms.h:949
Algorithms::AxisBundle::py
valarray< T > py(size_t i) const
Definition: lib-algorithms.h:127
Clock::Clock
Clock(int tout_start, double dt_out, double CFL)
Definition: lib-algorithms.h:930
Algorithms::AxisBundle::AxisBundle
AxisBundle(const vector< T > _xmin, const vector< T > _xmax, const vector< size_t > _Nx, const vector< T > _xgmin, const vector< T > _xgmax, const vector< size_t > _Nxg, const vector< T > _pmax, const vector< size_t > _Np, const vector< size_t > _Npx, const vector< size_t > _Npy, const vector< size_t > _Npz)
Definition: lib-algorithms.h:141
Algorithms::AxisBundle::_xg
vector< CAxis< T > > _xg
Definition: lib-algorithms.h:213
Algorithms::AxisBundle::_pz
vector< Axis< T > > _pz
Definition: lib-algorithms.h:212
Algorithms::AxisBundle::pmin
T pmin(size_t i) const
Definition: lib-algorithms.h:114
Algorithms::Legendre
valarray< T > Legendre(const T x, const size_t Nl)
Definition: lib-algorithms.h:227
Algorithms::AxisBundle::dp
T dp(size_t i) const
Definition: lib-algorithms.h:117
Algorithms::AxisBundle::px
valarray< T > px(size_t i) const
Definition: lib-algorithms.h:120
Algorithms::AxisBundle::dx
T dx(size_t i) const
Definition: lib-algorithms.h:102
Algorithms::Axis::Axis
Axis(const T min=0, const T max=0, const size_t N=1)
Definition: lib-algorithms.h:51
Algorithms::AxisBundle::Npy
size_t Npy(size_t i) const
Definition: lib-algorithms.h:128
Algorithms::CAxis::v
valarray< T > * v
Definition: lib-algorithms.h:72
Algorithms::AxisBundle::xg
valarray< T > xg(size_t i) const
Definition: lib-algorithms.h:105
Algorithms::AxisBundle::pxmin
T pxmin(size_t i) const
Definition: lib-algorithms.h:122
Algorithms::AxisBundle::pymax
T pymax(size_t i) const
Definition: lib-algorithms.h:130
Algorithms::AxisBundle::xgmin
T xgmin(size_t i) const
Definition: lib-algorithms.h:107
Algorithms::AxisBundle::pymin
T pymin(size_t i) const
Definition: lib-algorithms.h:129
Clock::tick
size_t tick() const
Definition: lib-algorithms.h:940
Algorithms::AxisBundle::AxisBundle
AxisBundle(const AxisBundle &a)
Definition: lib-algorithms.h:176
Algorithms::AxisBundle::pzmin
T pzmin(size_t i) const
Definition: lib-algorithms.h:136
Clock::numh
size_t numh() const
Definition: lib-algorithms.h:939
Algorithms::AxisBundle::xgmax
T xgmax(size_t i) const
Definition: lib-algorithms.h:108
Algorithms::Axis::Axis
Axis(const Axis &other)
Definition: lib-algorithms.h:59
Algorithms::relativistic_invg_moment
T relativistic_invg_moment(const vector< T > q, const valarray< T > x, const int p)
Definition: lib-algorithms.h:382
Clock::time
double time() const
Definition: lib-algorithms.h:941
Algorithms::AxisBundle::pxmax
T pxmax(size_t i) const
Definition: lib-algorithms.h:123
Clock::h
double h() const
Definition: lib-algorithms.h:938
Algorithms::AxisBundle::Np
size_t Np(size_t i) const
Definition: lib-algorithms.h:113
Algorithms::MakeAxis
valarray< T > MakeAxis(const T min, const T max, const size_t N)
Definition: lib-algorithms.h:25
Clock::_numh
size_t _numh
Definition: lib-algorithms.h:948
Algorithms::AxisBundle::xmax
T xmax(size_t i) const
Definition: lib-algorithms.h:100
Algorithms::Axis::v
valarray< T > * v
Definition: lib-algorithms.h:48
Clock::_hstep
size_t _hstep
Definition: lib-algorithms.h:950
Algorithms::CAxis::CAxis
CAxis(const T min=0, const T max=0, const size_t N=1)
Definition: lib-algorithms.h:75
Algorithms::AxisBundle::pzmax
T pzmax(size_t i) const
Definition: lib-algorithms.h:137
Algorithms::relativistic_gamma_moment
T relativistic_gamma_moment(const valarray< T > q, const valarray< T > x, const int p)
Definition: lib-algorithms.h:445
Algorithms::AxisBundle::p
valarray< T > p(size_t i) const
Definition: lib-algorithms.h:112
Algorithms::CAxis::CAxis
CAxis(const CAxis &other)
Definition: lib-algorithms.h:83
Algorithms::AxisBundle::Npx
size_t Npx(size_t i) const
Definition: lib-algorithms.h:121
Algorithms::MakeCAxis
valarray< T > MakeCAxis(const T min, const T max, const size_t N)
Definition: lib-algorithms.h:36
Algorithms::AxisBundle::pmax
T pmax(size_t i) const
Definition: lib-algorithms.h:115
Algorithms::AxisBundle::Npz
size_t Npz(size_t i) const
Definition: lib-algorithms.h:135
Clock::operator++
Clock & operator++()
Definition: lib-algorithms.h:944
Algorithms::AxisBundle::Nx
size_t Nx(size_t i) const
Definition: lib-algorithms.h:98
Algorithms::AxisBundle::Nxg
size_t Nxg(size_t i) const
Definition: lib-algorithms.h:106
Algorithms::AxisBundle::xmin
T xmin(size_t i) const
Definition: lib-algorithms.h:99
Algorithms::AxisBundle::x
valarray< T > x(size_t i) const
Definition: lib-algorithms.h:97
Algorithms::AxisBundle::pz
valarray< T > pz(size_t i) const
Definition: lib-algorithms.h:134
Algorithms::moment
T moment(const vector< T > q, const vector< T > x, const int p)
Definition: lib-algorithms.h:305
Clock::_t_start
double _t_start
Definition: lib-algorithms.h:947