Line data Source code
1 : #include "crpropa/PhotonBackground.h"
2 : #include "crpropa/Units.h"
3 : #include "crpropa/Random.h"
4 : #include "crpropa/Common.h"
5 :
6 : #include "kiss/logger.h"
7 : #include "kiss/path.h"
8 :
9 : #include <fstream>
10 :
11 : namespace crpropa {
12 :
13 95 : TabularPhotonField::TabularPhotonField(std::string fieldName, bool isRedshiftDependent) {
14 95 : this->fieldName = fieldName;
15 95 : this->isRedshiftDependent = isRedshiftDependent;
16 95 : this->isPositionDependent = false;
17 : this->surface = NULL;
18 :
19 285 : readPhotonEnergy(getDataPath("") + "Scaling/" + this->fieldName + "_photonEnergy.txt");
20 285 : readPhotonDensity(getDataPath("") + "Scaling/" + this->fieldName + "_photonDensity.txt");
21 :
22 95 : if (this->isRedshiftDependent)
23 296 : readRedshift(getDataPath("") + "Scaling/" + this->fieldName + "_redshift.txt");
24 :
25 95 : checkInputData();
26 :
27 95 : if (this->isRedshiftDependent)
28 74 : initRedshiftScaling();
29 95 : }
30 :
31 :
32 9541130 : double TabularPhotonField::getPhotonDensity(double Ephoton, double z, const Vector3d &pos) const {
33 9541130 : if (this->isRedshiftDependent) {
34 : // fix behaviour for future redshift. See issue #414
35 : // with redshift < 0 the photon density is set to 0 in interpolate2d.
36 : // Therefore it is assumed that the photon density does not change from values at z = 0. This is only valid for small changes in redshift.
37 9541130 : double zMin = this->redshifts[0];
38 :
39 9541130 : if(z < zMin){
40 0 : if(z < -1) {
41 0 : KISS_LOG_WARNING << "Photon Field " << fieldName << " uses FutureRedshift with z < -1. The photon density is set to n(Ephoton, z=0). \n";
42 : }
43 0 : return getPhotonDensity(Ephoton, zMin);
44 : } else {
45 9541130 : return interpolate2d(Ephoton, z, this->photonEnergies, this->redshifts, this->photonDensity);
46 : }
47 : } else {
48 0 : return interpolate(Ephoton, this->photonEnergies, this->photonDensity);
49 : }
50 : }
51 :
52 3988 : double TabularPhotonField::getRedshiftScaling(double z) const {
53 3988 : if (!this->isRedshiftDependent)
54 : return 1.;
55 3708 : if (z < this->redshifts.front())
56 : return 1.;
57 3708 : if (z > this->redshifts.back())
58 : return 0.;
59 3708 : return interpolate(z, this->redshifts, this->redshiftScalings);
60 : }
61 :
62 0 : double TabularPhotonField::getMinimumPhotonEnergy(double z, const Vector3d &pos) const{
63 0 : return photonEnergies[0];
64 : }
65 :
66 0 : double TabularPhotonField::getMaximumPhotonEnergy(double z, const Vector3d &pos) const{
67 0 : return photonEnergies[photonEnergies.size() -1];
68 : }
69 :
70 95 : void TabularPhotonField::readPhotonEnergy(std::string filePath) {
71 95 : std::ifstream infile(filePath.c_str());
72 :
73 95 : if (!infile.good())
74 0 : throw std::runtime_error("TabularPhotonField::readPhotonEnergy: could not open " + filePath);
75 :
76 : std::string line;
77 15425 : while (std::getline(infile, line)) {
78 15330 : if ((line.size() > 0) & (line[0] != '#') )
79 15045 : this->photonEnergies.push_back(std::stod(line));
80 : }
81 95 : infile.close();
82 95 : }
83 :
84 95 : void TabularPhotonField::readPhotonDensity(std::string filePath) {
85 95 : std::ifstream infile(filePath.c_str());
86 :
87 95 : if (!infile.good())
88 0 : throw std::runtime_error("TabularPhotonField::readPhotonDensity: could not open " + filePath);
89 :
90 : std::string line;
91 4774543 : while (std::getline(infile, line)) {
92 4774448 : if ((line.size() > 0) & (line[0] != '#') )
93 4774163 : this->photonDensity.push_back(std::stod(line));
94 : }
95 95 : infile.close();
96 95 : }
97 :
98 74 : void TabularPhotonField::readRedshift(std::string filePath) {
99 74 : std::ifstream infile(filePath.c_str());
100 :
101 74 : if (!infile.good())
102 0 : throw std::runtime_error("TabularPhotonField::initRedshift: could not open " + filePath);
103 :
104 : std::string line;
105 14623 : while (std::getline(infile, line)) {
106 14549 : if ((line.size() > 0) & (line[0] != '#') )
107 14327 : this->redshifts.push_back(std::stod(line));
108 : }
109 74 : infile.close();
110 74 : }
111 :
112 74 : void TabularPhotonField::initRedshiftScaling() {
113 : double n0 = 0.;
114 :
115 14401 : for (int i = 0; i < this->redshifts.size(); ++i) {
116 :
117 14327 : double z = this->redshifts[i];
118 : double n = 0.;
119 :
120 4772042 : for (int j = 0; j < this->photonEnergies.size()-1; ++j) {
121 4757715 : double e_j = this->photonEnergies[j];
122 4757715 : double e_j1 = this->photonEnergies[j+1];
123 4757715 : double deltaLogE = std::log10(e_j1) - std::log10(e_j);
124 4757715 : if (z == 0.)
125 12850 : n0 += (getPhotonDensity(e_j, 0) + getPhotonDensity(e_j1, 0)) / 2. * deltaLogE;
126 4757715 : n += (getPhotonDensity(e_j, z) + getPhotonDensity(e_j1, z)) / 2. * deltaLogE;
127 : }
128 14327 : this->redshiftScalings.push_back(n / n0);
129 : }
130 74 : }
131 :
132 95 : void TabularPhotonField::checkInputData() const {
133 95 : if (this->isRedshiftDependent) {
134 :
135 74 : if (this->photonDensity.size() != this->photonEnergies.size() * this-> redshifts.size())
136 0 : throw std::runtime_error("TabularPhotonField::checkInputData: length of photon density input is unequal to length of photon energy input times length of redshift input");
137 : } else {
138 21 : if (this->photonEnergies.size() != this->photonDensity.size())
139 0 : throw std::runtime_error("TabularPhotonField::checkInputData: length of photon energy input is unequal to length of photon density input");
140 : }
141 :
142 15140 : for (int i = 0; i < this->photonEnergies.size(); ++i) {
143 : double ePrevious = 0.;
144 15045 : double e = this->photonEnergies[i];
145 :
146 15045 : if (e <= 0.)
147 0 : throw std::runtime_error("TabularPhotonField::checkInputData: a value in the photon energy input is not positive");
148 : if (e <= ePrevious)
149 : throw std::runtime_error("TabularPhotonField::checkInputData: photon energy values are not strictly increasing");
150 : ePrevious = e;
151 : }
152 :
153 4774258 : for (int i = 0; i < this->photonDensity.size(); ++i) {
154 4774163 : if (this->photonDensity[i] < 0.)
155 0 : throw std::runtime_error("TabularPhotonField::checkInputData: a value in the photon density input is negative");
156 : }
157 :
158 95 : if (this->isRedshiftDependent) {
159 :
160 74 : if (this->redshifts[0] != 0.)
161 0 : throw std::runtime_error("TabularPhotonField::checkInputData: redshift input must start with zero");
162 :
163 14401 : for (int i = 0; i < this->redshifts.size(); ++i) {
164 : double zPrevious = -1.;
165 14327 : double z = this->redshifts[i];
166 :
167 14327 : if (z < 0.)
168 0 : throw std::runtime_error("TabularPhotonField::checkInputData: a value in the redshift input is negative");
169 : if (z <= zPrevious)
170 : throw std::runtime_error("TabularPhotonField::checkInputData: redshift values are not strictly increasing");
171 : zPrevious = z;
172 : }
173 :
174 74 : for (int i = 0; i < this->redshiftScalings.size(); ++i) {
175 0 : double scalingFactor = this->redshiftScalings[i];
176 0 : if (scalingFactor <= 0.)
177 0 : throw std::runtime_error("TabularPhotonField::checkInputData: initRedshiftScaling has created a non-positive scaling factor");
178 : }
179 : }
180 95 : }
181 :
182 0 : TabularSpatialPhotonField::TabularSpatialPhotonField(std::string fieldName, ref_ptr<Surface> surface) {
183 :
184 0 : this->fieldName = fieldName;
185 0 : this->isRedshiftDependent = isRedshiftDependent;
186 0 : this->isPositionDependent = true;
187 0 : this->surface = surface;
188 :
189 0 : std::string dirE = getDataPath("") + "Scaling/" + this->fieldName + "/photonEnergy/";
190 0 : if (!is_directory(dirE)) {
191 : std::cout << "Photon tables not found in " << dirE << std::endl;
192 : return;
193 : }
194 :
195 : std::unordered_map<int, Vector3d> photonDict;
196 0 : int iFile = 0;
197 :
198 : std::vector<std::string> dirsE;
199 0 : if(!list_directory(dirE, dirsE))
200 0 : throw std::runtime_error("Could not find any files in " + dirE + "!\n");
201 :
202 0 : for (auto const& dir_entry : dirsE) {
203 0 : std::vector<double> vE = readPhotonEnergy(concat_path(dirE, dir_entry));
204 :
205 0 : this->photonEnergies = vE;
206 : break;
207 0 : }
208 :
209 0 : std::string dirD = getDataPath("") + "Scaling/" + this->fieldName + "/photonDensity/";
210 0 : if (!is_directory(dirD)) {
211 : std::cout << "Photon tables not found in " << dirD << std::endl;
212 : return;
213 : }
214 :
215 : std::vector<std::string> dirsD;
216 0 : if(!list_directory(dirD, dirsD))
217 0 : throw std::runtime_error("Could not find any files in " + dirD + "!\n");
218 :
219 0 : for (auto const& dir_entry : dirsD) {
220 : double x, y, z;
221 : std::string str;
222 0 : std::stringstream ss;
223 :
224 0 : std::string filename = splitFilename(dir_entry);
225 : ss << filename;
226 :
227 : //Getline function to take and store the x, y, z coordinates of each node
228 : int iLine = 0;
229 : // it ensures the double numbers are of the type 1.00329, with the . for the decimal part
230 0 : std::locale::global(std::locale("C"));
231 :
232 0 : while (getline(ss, str, '_')) {
233 0 : if (iLine == 2) {
234 0 : x = stod(str) * kpc;
235 : }
236 0 : if (iLine == 3) {
237 0 : y = stod(str) * kpc;
238 : }
239 0 : if (iLine == 4) {
240 0 : z = stod(str) * kpc;
241 : }
242 0 : iLine = iLine + 1;
243 : }
244 :
245 : Vector3d vPos(x, y, z);
246 :
247 : // Continue when not "inside" surface
248 0 : if (getSurface() && getSurface()->distance(vPos)>=0)
249 : continue;
250 :
251 : photonDict[iFile] = vPos;
252 0 : iFile = iFile + 1;
253 :
254 0 : std::vector<double> vD = readPhotonDensity(concat_path(dirD, dir_entry));
255 0 : this->photonDensity.push_back(vD);
256 0 : }
257 :
258 0 : if (this->photonDensity.empty())
259 0 : throw std::runtime_error("Tabular spatial photon field for " + fieldName + " empty! Check if the surface is properly set.");
260 :
261 : this->photonDict = photonDict;
262 0 : checkInputData();
263 0 : }
264 :
265 0 : void TabularSpatialPhotonField::setSurface(ref_ptr<Surface> surface) {
266 0 : this->surface = surface;
267 0 : }
268 :
269 0 : double TabularSpatialPhotonField::getPhotonDensity(const double ePhoton, double z, const Vector3d &pos) const {
270 : double dMin = 1000.;
271 : int iMin = -1;
272 :
273 0 : for (const auto& el : this->photonDict) {
274 : Vector3d posNode = el.second;
275 : double d;
276 0 : d = sqrt((- posNode.x / kpc - pos.x / kpc) * (- posNode.x / kpc - pos.x / kpc) + (posNode.y / kpc - pos.y / kpc) * (posNode.y / kpc - pos.y / kpc ) + (posNode.z / kpc - pos.z / kpc) * (posNode.z / kpc - pos.z / kpc));
277 :
278 0 : if (d<dMin) {
279 : dMin = d;
280 0 : iMin = el.first;
281 : }
282 : }
283 :
284 0 : if (iMin == -1) {
285 : return -1.;
286 : } else {
287 0 : if ((ePhoton < photonEnergies[0]) || (ePhoton > photonEnergies[photonEnergies.size() - 1])) {
288 : return 0;
289 : } else {
290 0 : std::vector<double> rowE = this->photonEnergies; // assuming all the nodes have the same energy binning
291 0 : std::vector<double> rowD = this->photonDensity[iMin];
292 0 : return interpolate(ePhoton, rowE, rowD);
293 0 : }
294 : }
295 : }
296 :
297 0 : double TabularSpatialPhotonField::getMinimumPhotonEnergy(double z, const Vector3d &pos) const {
298 0 : return photonEnergies[0]; // assuming all the nodes have the same energy bins
299 : }
300 :
301 0 : double TabularSpatialPhotonField::getMaximumPhotonEnergy(double z, const Vector3d &pos) const {
302 0 : return photonEnergies[photonEnergies.size() - 1]; // assuming all the nodes have the same energy bins
303 : }
304 :
305 0 : std::vector<double> TabularSpatialPhotonField::readPhotonEnergy(std::string filePath) {
306 0 : std::ifstream infile(filePath.c_str());
307 :
308 0 : if (!infile.good())
309 0 : throw std::runtime_error("TabularPhotonField::readPhotonEnergy: could not open " + filePath);
310 :
311 : std::string line;
312 : std::vector<double> vE;
313 :
314 0 : while (std::getline(infile, line)) {
315 0 : if ((line.size() > 0) & (line[0] != '#') ) {
316 0 : vE.insert(vE.begin(),std::stod(line));
317 : }
318 : }
319 :
320 0 : infile.close();
321 0 : return vE;
322 0 : }
323 :
324 0 : std::vector<double> TabularSpatialPhotonField::readPhotonDensity(std::string filePath) {
325 0 : std::ifstream infile(filePath.c_str());
326 :
327 0 : if (!infile.good())
328 0 : throw std::runtime_error("TabularPhotonField::readPhotonDensity: could not open " + filePath);
329 :
330 : std::string line;
331 : std::vector<double> vD;
332 :
333 0 : while (std::getline(infile, line)) {
334 0 : if ((line.size() > 0) & (line[0] != '#') )
335 0 : vD.insert(vD.begin(),std::stod(line));
336 : }
337 :
338 0 : infile.close();
339 0 : return vD;
340 0 : }
341 :
342 0 : void TabularSpatialPhotonField::checkInputData() const {
343 :
344 : std::size_t numRowsDens = this->photonEnergies.size();
345 : std::size_t numRowsEn = this->photonDensity.size();
346 :
347 0 : for (int j = 0; j < this->photonDensity.size(); ++j) { //take the proper row size!
348 :
349 0 : if (this->photonEnergies.size() != this->photonDensity[j].size())
350 0 : throw std::runtime_error("TabularPhotonField::checkInputData: length of photon energy input is unequal to length of photon density input");
351 :
352 0 : for (int i = 0; i < this->photonEnergies.size(); ++i) {
353 : double ePrevious = 0.;
354 0 : double e = this->photonEnergies[i];
355 :
356 0 : if (e <= 0.)
357 0 : throw std::runtime_error("TabularSpatialPhotonField::checkInputData: a value in the photon energy input is not positive");
358 : if (e <= ePrevious)
359 : throw std::runtime_error("TabularSpatialPhotonField::checkInputData: photon energy values are not strictly increasing");
360 : ePrevious = e;
361 : }
362 :
363 0 : for (int i = 0; i < this->photonDensity[j].size(); ++i) {
364 :
365 0 : if (this->photonDensity[j][i] < 0.)
366 0 : throw std::runtime_error("TabularSpatialPhotonField::checkInputData: a value in the photon density input is negative");
367 : }
368 : }
369 0 : }
370 :
371 47 : BlackbodyPhotonField::BlackbodyPhotonField(std::string fieldName, double blackbodyTemperature) {
372 47 : this->fieldName = fieldName;
373 47 : this->blackbodyTemperature = blackbodyTemperature;
374 47 : this->quantile = 0.0001; // tested to be sufficient, only used for extreme values of primary energy or temperature
375 47 : }
376 :
377 9626 : double BlackbodyPhotonField::getPhotonDensity(double Ephoton, double z, const Vector3d &pos) const {
378 9626 : return 8 * M_PI * pow_integer<3>(Ephoton / (h_planck * c_light)) / std::expm1(Ephoton / (k_boltzmann * this->blackbodyTemperature));
379 : }
380 :
381 26 : double BlackbodyPhotonField::getMinimumPhotonEnergy(double z, const Vector3d &pos) const {
382 : double A;
383 26 : int quantile_int = 10000 * quantile;
384 26 : switch (quantile_int)
385 : {
386 : case 1: // 0.01 % percentil
387 : A = 1.093586e-5 * eV / kelvin;
388 : break;
389 0 : case 10: // 0.1 % percentil
390 : A = 2.402189e-5 * eV / kelvin;
391 0 : break;
392 0 : case 100: // 1 % percentil
393 : A = 5.417942e-5 * eV / kelvin;
394 0 : break;
395 0 : default:
396 0 : throw std::runtime_error("Quantile not understood. Please use 0.01 (1%), 0.001 (0.1%) or 0.0001 (0.01%) \n");
397 : break;
398 : }
399 26 : return A * this -> blackbodyTemperature;
400 : }
401 :
402 51 : double BlackbodyPhotonField::getMaximumPhotonEnergy(double z, const Vector3d &pos) const {
403 51 : double factor = std::max(1., blackbodyTemperature / 2.73);
404 51 : return 0.1 * factor * eV; // T dependent scaling, starting at 0.1 eV as suitable for CMB
405 : }
406 :
407 0 : void BlackbodyPhotonField::setQuantile(double q) {
408 0 : if(not ((q == 0.0001) or (q == 0.001) or (q == 0.01)))
409 0 : throw std::runtime_error("Quantile not understood. Please use 0.01 (1%), 0.001 (0.1%) or 0.0001 (0.01%) \n");
410 0 : this -> quantile = q;
411 0 : }
412 :
413 : } // namespace crpropa
|
