Line data    Source code

       1             : #include "crpropa/Common.h"
       2             : 
       3             : #include "kiss/path.h"
       4             : #include "kiss/logger.h"
       5             : 
       6             : #include <cstdlib>
       7             : #include <fstream>
       8             : #include <string>
       9             : #include <cmath>
      10             : #include <algorithm>
      11             : 
      12             : #define index(i,j) ((j)+(i)*Y.size())
      13             : 
      14             : namespace crpropa {
      15             : 
      16         547 : std::string getDataPath(std::string filename) {
      17         547 :         static std::string dataPath;
      18         547 :         if (dataPath.size())
      19         534 :                 return concat_path(dataPath, filename);
      20             : 
      21          13 :         const char *env_path = getenv("CRPROPA_DATA_PATH");
      22          13 :         if (env_path) {
      23           0 :                 if (is_directory(env_path)) {
      24             :                         dataPath = env_path;
      25           0 :                         KISS_LOG_INFO << "getDataPath: use environment variable, "
      26           0 :                                         << dataPath << std::endl;
      27           0 :                         return concat_path(dataPath, filename);
      28             :                 }
      29             :         }
      30             : 
      31             : #ifdef CRPROPA_INSTALL_PREFIX
      32             :         {
      33          13 :                 std::string _path = CRPROPA_INSTALL_PREFIX "/share/crpropa";
      34          13 :                 if (is_directory(_path)) {
      35             :                         dataPath = _path;
      36           0 :                         KISS_LOG_INFO
      37           0 :                         << "getDataPath: use install prefix, " << dataPath << std::endl;
      38           0 :                         return concat_path(dataPath, filename);
      39             :                 }
      40             :         }
      41             : #endif
      42             : 
      43             :         {
      44          13 :                 std::string _path = executable_path() + "../data";
      45          13 :                 if (is_directory(_path)) {
      46             :                         dataPath = _path;
      47           0 :                         KISS_LOG_INFO << "getDataPath: use executable path, " << dataPath
      48           0 :                                         << std::endl;
      49           0 :                         return concat_path(dataPath, filename);
      50             :                 }
      51             :         }
      52             : 
      53             :         dataPath = "data";
      54          13 :         KISS_LOG_INFO << "getDataPath: use default, " << dataPath << std::endl;
      55          13 :         return concat_path(dataPath, filename);
      56             : }
      57             : 
      58             : 
      59           0 : std::string getInstallPrefix()
      60             : {
      61           0 :   std::string _path = "";
      62             :   #ifdef CRPROPA_INSTALL_PREFIX
      63             :     _path += CRPROPA_INSTALL_PREFIX;
      64             :   #endif
      65           0 :   return _path;
      66             : };
      67             : 
      68      307611 : double interpolate(double x, const std::vector<double> &X,
      69             :                 const std::vector<double> &Y) {
      70             :         std::vector<double>::const_iterator it = std::upper_bound(X.begin(),
      71             :                         X.end(), x);
      72      307611 :         if (it == X.begin())
      73           1 :                 return Y.front();
      74      307610 :         if (it == X.end())
      75           1 :                 return Y.back();
      76             : 
      77             :         size_t i = it - X.begin() - 1;
      78      307609 :         return Y[i] + (x - X[i]) * (Y[i + 1] - Y[i]) / (X[i + 1] - X[i]);
      79             : }
      80             : 
      81     9586406 : double interpolate2d(double x, double y, const std::vector<double> &X,
      82             :                 const std::vector<double> &Y, const std::vector<double> &Z) {
      83             : 
      84             :         std::vector<double>::const_iterator itx = std::upper_bound(X.begin(), X.end(), x);
      85             :         std::vector<double>::const_iterator ity = std::upper_bound(Y.begin(), Y.end(), y);
      86             : 
      87     9586406 :         if (x > X.back() || x < X.front())
      88             :                 return 0;
      89     9586384 :         if (y > Y.back() || y < Y.front())
      90             :                 return 0;
      91             : 
      92     9586384 :         if (itx == X.begin() && ity == Y.begin())
      93           0 :                 return Z.front();
      94     9586384 :         if (itx == X.end() && ity == Y.end())
      95          73 :                 return Z.back();
      96             : 
      97     9586311 :         size_t i = itx - X.begin() - 1;
      98     9586311 :         size_t j = ity - Y.begin() - 1;
      99             : 
     100     9586311 :         double Q11 = Z[index(i,j)];
     101     9586311 :         double Q12 = Z[index(i,j+1)];
     102     9586311 :         double Q21 = Z[index(i+1,j)];
     103     9586311 :         double Q22 = Z[index(i+1,j+1)];
     104             : 
     105     9586311 :         double R1 = ((X[i+1]-x)/(X[i+1]-X[i]))*Q11+((x-X[i])/(X[i+1]-X[i]))*Q21;
     106     9586311 :         double R2 = ((X[i+1]-x)/(X[i+1]-X[i]))*Q12+((x-X[i])/(X[i+1]-X[i]))*Q22;
     107             : 
     108     9586311 :         return ((Y[j+1]-y)/(Y[j+1]-Y[j]))*R1+((y-Y[j])/(Y[j+1]-Y[j]))*R2;
     109             : }
     110             : 
     111        1203 : double interpolateEquidistant(double x, double lo, double hi,
     112             :                 const std::vector<double> &Y) {
     113        1203 :         if (x <= lo)
     114           1 :                 return Y.front();
     115        1202 :         if (x >= hi)
     116           1 :                 return Y.back();
     117             : 
     118        1201 :         double dx = (hi - lo) / (Y.size() - 1);
     119        1201 :         double p = (x - lo) / dx;
     120        1201 :         size_t i = floor(p);
     121        1201 :         return Y[i] + (p - i) * (Y[i + 1] - Y[i]);
     122             : }
     123             : 
     124         561 : size_t closestIndex(double x, const std::vector<double> &X) {
     125         561 :         size_t i1 = std::lower_bound(X.begin(), X.end(), x) - X.begin();
     126         561 :         if (i1 == 0)
     127             :                 return i1;
     128         561 :         size_t i0 = i1 - 1;
     129         561 :         if (std::fabs(X[i0] - x) < std::fabs(X[i1] - x))
     130             :                 return i0;
     131             :         else
     132         304 :                 return i1;
     133             : }
     134             : 
     135             : } // namespace crpropa
     136             :