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              : // Maximum nuclear charge number, neutron number, and flat-array stride
      18              : // for all nuclear interaction tables (PhotoDisintegration, NuclearDecay, ParticleMass).
      19              : // Coverage matches PD_Talys1.9: Z up to Pb (82), N up to Pb-214 (132).
      20              : static const int NUCLEAR_ZMAX    = 82;
      21              : static const int NUCLEAR_NMAX    = 132;
      22              : static const int NUCLEAR_NSTRIDE = 133; // = NUCLEAR_NMAX + 1
      23              : 
      24              : // Returns the full path to a CRPropa data file
      25              : std::string getDataPath(std::string filename);
      26              : 
      27              : // Returns the install prefix
      28              : std::string getInstallPrefix();
      29              : 
      30              : // Returns a certain digit from a given integer
      31              : inline int digit(const int& value, const int& d) {
      32        55566 :         return (value % (d * 10)) / d;
      33              : }
      34              : 
      35              : // Return value xclip which is the closest to x, so that lower <= xclip <= upper
      36              : template <typename T>
      37              : T clip(const T& x, const T& lower, const T& upper) {
      38       492352 :         return std::max(lower, std::min(x, upper));
      39              : }
      40              : 
      41              : // Perform linear interpolation on a set of n tabulated data points X[0 .. n-1] -> Y[0 .. n-1]
      42              : // Returns Y[0] if x < X[0] and Y[n-1] if x > X[n-1]
      43              : double interpolate(double x, const std::vector<double>& X,
      44              :                 const std::vector<double>& Y);
      45              : 
      46              : 
      47              : // 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]
      48              : // Returns 0 if x < X[0] or x > X[n-1] or y < Y[0] or y > Y[m-1]
      49              : double interpolate2d(double x, double y, const std::vector<double>& X,
      50              :                 const std::vector<double>& Y, const std::vector<double>& Z);
      51              : 
      52              : // Perform linear interpolation on equidistant tabulated data
      53              : // Returns Y[0] if x < lo and Y[n-1] if x > hi
      54              : double interpolateEquidistant(double x, double lo, double hi,
      55              :                 const std::vector<double>& Y);
      56              : 
      57              : // Find index of value in a sorted vector X that is closest to x
      58              : size_t closestIndex(double x, const std::vector<double> &X);
      59              : 
      60              : // Takes the filename from a full data path. Used in EM* modules.
      61              : std::string splitFilename(const std::string str);
      62              : /** @}*/
      63              : 
      64              : // pow implementation as template for integer exponents pow_integer<2>(x)
      65              : // evaluates to x*x
      66              : template <unsigned int exponent>
      67              : inline double pow_integer(double base)
      68              : {
      69        35430 :   return pow_integer<(exponent >> 1)>(base*base) * (((exponent & 1) > 0) ? base : 1);
      70              : }
      71              : 
      72              : template <>
      73              : inline double pow_integer<0>(double base)
      74              : {
      75              :   return 1;
      76              : }
      77              : 
      78              : // - input:  function over which to integrate, integration limits A and B
      79              : // - output: 8-points Gauß-Legendre integral
      80              : static const double X[8] = {.0950125098, .2816035507, .4580167776, .6178762444, .7554044083, .8656312023, .9445750230, .9894009349};
      81              : static const double W[8] = {.1894506104, .1826034150, .1691565193, .1495959888, .1246289712, .0951585116, .0622535239, .0271524594};
      82              : template<typename Integrand>
      83         9601 : double gaussInt(Integrand&& integrand, double A, double B) {
      84         9601 :         const double XM = 0.5 * (B + A);
      85         9601 :         const double XR = 0.5 * (B - A);
      86              :         double SS = 0.;
      87        86418 :         for (int i = 0; i < 8; ++i) {
      88        76816 :                 double DX = XR * X[i];
      89        76816 :                 SS += W[i] * (integrand(XM + DX) + integrand(XM - DX));
      90              :         }
      91         9601 :         return XR * SS;
      92              : }
      93              : 
      94              : } // namespace crpropa
      95              : 
      96              : #endif // CRPROPA_COMMON_H