20 #include <gsl/gsl_errno.h>
21 #include <gsl/gsl_spline.h>
23 #include "../NumtoolsLib/NumtoolsLib.h"
31 const vector<double>& vec_interpretation,
35 _vec_interpretation(vec_interpretation)
43 if( v > _par_ode._par_pop._theta )
45 if( v < _par_ode._V_min )
47 Number n = _vec_interpretation.size();
48 if ( v > _vec_interpretation.back() )
52 if ( v < _vec_interpretation[
i] )
59 int sg = (dv > 0) - (dv < 0);
60 int n =
static_cast<int>(_vec_interpretation.size());
61 for (
int i = 0;
i < n;
i++){
63 if (_vec_interpretation[j] <= v && (j == static_cast<MPILib::Index>(n-1) || _vec_interpretation[j+1] > v) )
73 if ( v < _par_ode._V_min)
75 int n =
static_cast<int>(_vec_interpretation.size());
76 if ( v > _par_ode._par_pop._theta)
82 int i_ret = Search(i,v, dv);
93 int n = _vec_interpretation.size();
94 if (i_tr_low < 0 || i_tr_low > n-1)
97 Potential low_bin = _vec_interpretation[i_tr_low];
98 Potential high_bin = (i_tr_low !=
static_cast<int>(_vec_interpretation.size()) - 1)? _vec_interpretation[i_tr_low + 1] : _par_ode._par_pop._theta;
99 assert( low_bin <= v && high_bin >= v);
100 return (high_bin -v )/(high_bin - low_bin);
104 int n = _vec_interpretation.size();
105 if (i_tr_high < 0 || i_tr_high > n-1)
108 Potential low_bin = _vec_interpretation[i_tr_high];
109 Potential high_bin = (i_tr_high !=
static_cast<int>(_vec_interpretation.size()) - 1)? _vec_interpretation[i_tr_high + 1] : _par_ode._par_pop._theta;
110 assert( low_bin <= v && high_bin >= v);
111 return (v - low_bin)/(high_bin - low_bin);
115 assert(i < _vec_interpretation.size());
118 Potential high = (i != _vec_interpretation.size() - 1) ? _vec_interpretation[i+1] : _par_ode._par_pop._theta;
120 Potential trans_low = Translate(low, delta_v);
121 Potential trans_high = Translate(high, delta_v);
124 int i_tr_low = this->SearchBin(i,trans_low, delta_v);
125 int i_tr_high = this->SearchBin(i,trans_high,delta_v);
127 double alpha_low = BinLowFraction (trans_low, i_tr_low);
128 double alpha_high = BinHighFraction(trans_high, i_tr_high);
130 pair_ret.first._alpha = alpha_low;
131 pair_ret.first._index = i_tr_low;
132 pair_ret.second._alpha = alpha_high;
133 pair_ret.second._index = i_tr_high;
BinEstimator(const vector< Potential > &, const OdeParameter &)
Constructor.
Contains the parameters necessary to configure a concrete OdeSystem instance. See AbstractOdeSystem a...
double BinHighFraction(Potential, int) const
CoverPair CalculateBinCover(MPILib::Index i, Potential h) const
Cover pair looks at the boundaries of the current bins and calculates where they would end up if a po...
Base class for all exceptions thrown in GeomLib.
MPILib::Index SearchBin(Potential) const
Use for testing purposes, for any normal usage CalculateBinCover should be used. This function return...
int Search(MPILib::Index, Potential, Potential) const
std::pair< BinCover, BinCover > CoverPair
Type representing begin an end bins of a translated bins and the respective covering fraction...
Potential Translate(Potential, Potential) const
The objective is find a numerical solution for this equation This requires a numerical representation of the density We will work in the state space of a two dimensional and define a mesh there We first give two examples and then define the general procedure and given in Table for given fixed Delta g see Fig and we will denote coordinates in this dimension by a small letter $v The second dimension can be used to represent parameters as varied as and will represented by $w A strip is constructed by choosing two neighbouring points in state e g and integrating the vector field for a time $T that is assumed to be an integer multiple of a period of time Delta which we assume to be a defining characteristic of the grid Let then the set of points the set of points which is quadrilateral in shape The quadrilateral should be but not necessarily as long as they are but it is convenient to number them in order of creation In the we will assume that strip numbers created by the integration procedure start and are so that the numbers i in each identify a unique strip Strip no is reserved for stationary points There may or more cells in strip The number of cells in strip $i denoted by with $i the strip number and $j the cell as the i
double BinLowFraction(Potential, int) const