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 567 : std::string removeNullCharacter(std::string path) {
17 : // check for null character in data path and remove it
18 567 : if (path.find('\x00') != std::string::npos)
19 0 : path.erase(std::remove(path.begin(), path.end(), '\x00'), path.end());
20 567 : return path;
21 : }
22 :
23 554 : std::string getDataPath(std::string filename) {
24 554 : static std::string dataPath;
25 :
26 554 : if (dataPath.size())
27 1082 : return removeNullCharacter(concat_path(dataPath, filename));
28 :
29 13 : const char *env_path = getenv("CRPROPA_DATA_PATH");
30 13 : if (env_path) {
31 0 : if (is_directory(env_path)) {
32 0 : dataPath = removeNullCharacter(env_path);
33 0 : KISS_LOG_INFO << "getDataPath: use environment variable, "
34 0 : << dataPath << std::endl;
35 0 : return concat_path(dataPath, filename);
36 : }
37 : }
38 :
39 : #ifdef CRPROPA_INSTALL_PREFIX
40 : {
41 13 : std::string _path = CRPROPA_INSTALL_PREFIX "/share/crpropa";
42 13 : if (is_directory(_path)) {
43 0 : dataPath = removeNullCharacter(_path);
44 0 : KISS_LOG_INFO
45 0 : << "getDataPath: use install prefix, " << dataPath << std::endl;
46 0 : return concat_path(dataPath, filename);
47 : }
48 : }
49 : #endif
50 :
51 : {
52 26 : std::string _path = removeNullCharacter(executable_path() + "../data");
53 13 : if (is_directory(_path)) {
54 0 : dataPath = removeNullCharacter(_path);
55 0 : KISS_LOG_INFO << "getDataPath: use executable path, " << dataPath
56 0 : << std::endl;
57 0 : return concat_path(dataPath, filename);
58 : }
59 : }
60 :
61 : #ifdef CRPROPA_BINARY_DIR
62 : {
63 13 : std::string _path = CRPROPA_BINARY_DIR "/data";
64 13 : if (is_directory(_path)) {
65 13 : dataPath = removeNullCharacter(_path);
66 13 : KISS_LOG_INFO
67 0 : << "getDataPath: use binary path, " << dataPath << std::endl;
68 13 : return concat_path(dataPath, filename);
69 : }
70 : }
71 : #endif
72 :
73 : dataPath = "data";
74 0 : KISS_LOG_INFO << "getDataPath: use default, " << dataPath << std::endl;
75 0 : return removeNullCharacter(concat_path(dataPath, filename));
76 : }
77 :
78 :
79 0 : std::string getInstallPrefix()
80 : {
81 0 : std::string _path = "";
82 : #ifdef CRPROPA_INSTALL_PREFIX
83 : _path += CRPROPA_INSTALL_PREFIX;
84 : #endif
85 0 : return _path;
86 : };
87 :
88 114267 : double interpolate(double x, const std::vector<double> &X,
89 : const std::vector<double> &Y) {
90 114267 : std::vector<double>::const_iterator it = std::upper_bound(X.begin(),
91 : X.end(), x);
92 114267 : if (it == X.begin())
93 1 : return Y.front();
94 114266 : if (it == X.end())
95 1 : return Y.back();
96 :
97 : size_t i = it - X.begin() - 1;
98 114265 : return Y[i] + (x - X[i]) * (Y[i + 1] - Y[i]) / (X[i + 1] - X[i]);
99 : }
100 :
101 9539396 : double interpolate2d(double x, double y, const std::vector<double> &X,
102 : const std::vector<double> &Y, const std::vector<double> &Z) {
103 :
104 9539396 : std::vector<double>::const_iterator itx = std::upper_bound(X.begin(), X.end(), x);
105 9539396 : std::vector<double>::const_iterator ity = std::upper_bound(Y.begin(), Y.end(), y);
106 :
107 9539396 : if (x > X.back() || x < X.front())
108 : return 0;
109 9539395 : if (y > Y.back() || y < Y.front())
110 : return 0;
111 :
112 9539395 : if (itx == X.begin() && ity == Y.begin())
113 0 : return Z.front();
114 9539395 : if (itx == X.end() && ity == Y.end())
115 73 : return Z.back();
116 :
117 9539322 : size_t i = itx - X.begin() - 1;
118 9539322 : size_t j = ity - Y.begin() - 1;
119 :
120 9539322 : double Q11 = Z[index(i,j)];
121 9539322 : double Q12 = Z[index(i,j+1)];
122 9539322 : double Q21 = Z[index(i+1,j)];
123 9539322 : double Q22 = Z[index(i+1,j+1)];
124 :
125 9539322 : double R1 = ((X[i+1]-x)/(X[i+1]-X[i]))*Q11+((x-X[i])/(X[i+1]-X[i]))*Q21;
126 9539322 : double R2 = ((X[i+1]-x)/(X[i+1]-X[i]))*Q12+((x-X[i])/(X[i+1]-X[i]))*Q22;
127 :
128 9539322 : return ((Y[j+1]-y)/(Y[j+1]-Y[j]))*R1+((y-Y[j])/(Y[j+1]-Y[j]))*R2;
129 : }
130 :
131 698 : double interpolateEquidistant(double x, double lo, double hi,
132 : const std::vector<double> &Y) {
133 698 : if (x <= lo)
134 1 : return Y.front();
135 697 : if (x >= hi)
136 1 : return Y.back();
137 :
138 696 : double dx = (hi - lo) / (Y.size() - 1);
139 696 : double p = (x - lo) / dx;
140 696 : size_t i = floor(p);
141 696 : return Y[i] + (p - i) * (Y[i + 1] - Y[i]);
142 : }
143 :
144 563 : size_t closestIndex(double x, const std::vector<double> &X) {
145 563 : size_t i1 = std::lower_bound(X.begin(), X.end(), x) - X.begin();
146 563 : if (i1 == 0)
147 : return i1;
148 563 : size_t i0 = i1 - 1;
149 563 : if (std::fabs(X[i0] - x) < std::fabs(X[i1] - x))
150 : return i0;
151 : else
152 306 : return i1;
153 : }
154 :
155 : } // namespace crpropa
156 :
|
