Line data    Source code

       1              : #ifndef CRPROPA_COMMON_H
       2              : #define CRPROPA_COMMON_H
       3              : 
       4              : #include <string>
       5              : #include <vector>
       6              : /**
       7              :  @file
       8              :  @brief Common helper functions
       9              :  */
      10              : 
      11              : namespace crpropa {
      12              : /**
      13              :  * \addtogroup Core
      14              :  * @{
      15              :  */
      16              : 
      17              : // Returns the full path to a CRPropa data file
      18              : std::string getDataPath(std::string filename);
      19              : 
      20              : // Returns the install prefix
      21              : std::string getInstallPrefix();
      22              : 
      23              : // Returns a certain digit from a given integer
      24              : inline int digit(const int& value, const int& d) {
      25         1199 :         return (value % (d * 10)) / d;
      26              : }
      27              : 
      28              : // Return value xclip which is the closest to x, so that lower <= xclip <= upper
      29              : template <typename T>
      30              : T clip(const T& x, const T& lower, const T& upper) {
      31       492133 :         return std::max(lower, std::min(x, upper));
      32              : }
      33              : 
      34              : // Perform linear interpolation on a set of n tabulated data points X[0 .. n-1] -> Y[0 .. n-1]
      35              : // Returns Y[0] if x < X[0] and Y[n-1] if x > X[n-1]
      36              : double interpolate(double x, const std::vector<double>& X,
      37              :                 const std::vector<double>& Y);
      38              : 
      39              : 
      40              : // Perform bilinear interpolation on a set of (n,m) tabulated data points X[0 .. n-1], Y[0 .. m-1] -> Z[0.. n-1*m-1]
      41              : // Returns 0 if x < X[0] or x > X[n-1] or y < Y[0] or y > Y[m-1]
      42              : double interpolate2d(double x, double y, const std::vector<double>& X,
      43              :                 const std::vector<double>& Y, const std::vector<double>& Z);
      44              : 
      45              : // Perform linear interpolation on equidistant tabulated data
      46              : // Returns Y[0] if x < lo and Y[n-1] if x > hi
      47              : double interpolateEquidistant(double x, double lo, double hi,
      48              :                 const std::vector<double>& Y);
      49              : 
      50              : // Find index of value in a sorted vector X that is closest to x
      51              : size_t closestIndex(double x, const std::vector<double> &X);
      52              : /** @}*/
      53              : 
      54              : 
      55              : // pow implementation as template for integer exponents pow_integer<2>(x)
      56              : // evaluates to x*x
      57              : template <unsigned int exponent>
      58              : inline double pow_integer(double base)
      59              : {
      60        30946 :   return pow_integer<(exponent >> 1)>(base*base) * (((exponent & 1) > 0) ? base : 1);
      61              : }
      62              : 
      63              : template <>
      64              : inline double pow_integer<0>(double base)
      65              : {
      66              :   return 1;
      67              : }
      68              : 
      69              : // - input:  function over which to integrate, integration limits A and B
      70              : // - output: 8-points Gauß-Legendre integral
      71              : static const double X[8] = {.0950125098, .2816035507, .4580167776, .6178762444, .7554044083, .8656312023, .9445750230, .9894009349};
      72              : static const double W[8] = {.1894506104, .1826034150, .1691565193, .1495959888, .1246289712, .0951585116, .0622535239, .0271524594};
      73              : template<typename Integrand>
      74        11968 : double gaussInt(Integrand&& integrand, double A, double B) {
      75        11968 :         const double XM = 0.5 * (B + A);
      76        11968 :         const double XR = 0.5 * (B - A);
      77              :         double SS = 0.;
      78       107721 :         for (int i = 0; i < 8; ++i) {
      79        95752 :                 double DX = XR * X[i];
      80        95752 :                 SS += W[i] * (integrand(XM + DX) + integrand(XM - DX));
      81              :         }
      82        11968 :         return XR * SS;
      83              : }
      84              : 
      85              : } // namespace crpropa
      86              : 
      87              : #endif // CRPROPA_COMMON_H