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        1304 :         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      524299 :         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      166078 :   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       63808 : double gaussInt(Integrand&& integrand, double A, double B) {
      75       63808 :         const double XM = 0.5 * (B + A);
      76       63808 :         const double XR = 0.5 * (B - A);
      77             :         double SS = 0.;
      78      574281 :         for (int i = 0; i < 8; ++i) {
      79      510472 :                 double DX = XR * X[i];
      80      510472 :                 SS += W[i] * (integrand(XM + DX) + integrand(XM - DX));
      81             :         }
      82       63808 :         return XR * SS;
      83             : }
      84             : 
      85             : } // namespace crpropa
      86             : 
      87             : #endif // CRPROPA_COMMON_H