Line data Source code
1 : #include "crpropa/PhotonBackground.h"
2 : #include "crpropa/Units.h"
3 : #include "crpropa/Random.h"
4 :
5 : #include "kiss/logger.h"
6 :
7 : #include <fstream>
8 : #include <stdexcept>
9 : #include <limits>
10 : #include <cmath>
11 :
12 : namespace crpropa {
13 :
14 94 : TabularPhotonField::TabularPhotonField(std::string fieldName, bool isRedshiftDependent) {
15 94 : this->fieldName = fieldName;
16 94 : this->isRedshiftDependent = isRedshiftDependent;
17 :
18 282 : readPhotonEnergy(getDataPath("") + "Scaling/" + this->fieldName + "_photonEnergy.txt");
19 282 : readPhotonDensity(getDataPath("") + "Scaling/" + this->fieldName + "_photonDensity.txt");
20 94 : if (this->isRedshiftDependent)
21 292 : readRedshift(getDataPath("") + "Scaling/" + this->fieldName + "_redshift.txt");
22 :
23 94 : checkInputData();
24 :
25 94 : if (this->isRedshiftDependent)
26 73 : initRedshiftScaling();
27 94 : }
28 :
29 :
30 9539396 : double TabularPhotonField::getPhotonDensity(double Ephoton, double z) const {
31 9539396 : if ((this->isRedshiftDependent)) {
32 : // fix behaviour for future redshift. See issue #414
33 : // with redshift < 0 the photon density is set to 0 in interpolate2d.
34 : // 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.
35 9539396 : double zMin = this->redshifts[0];
36 9539396 : if(z < zMin){
37 0 : if(z < -1) {
38 0 : KISS_LOG_WARNING << "Photon Field " << fieldName << " uses FutureRedshift with z < -1. The photon density is set to n(Ephoton, z=0). \n";
39 : }
40 0 : return getPhotonDensity(Ephoton, zMin);
41 : } else {
42 9539396 : return interpolate2d(Ephoton, z, this->photonEnergies, this->redshifts, this->photonDensity);
43 : }
44 : } else {
45 0 : return interpolate(Ephoton, this->photonEnergies, this->photonDensity);
46 : }
47 : }
48 :
49 :
50 1787 : double TabularPhotonField::getRedshiftScaling(double z) const {
51 1787 : if (!this->isRedshiftDependent)
52 : return 1.;
53 :
54 1603 : if (z < this->redshifts.front())
55 : return 1.;
56 :
57 1603 : if (z > this->redshifts.back())
58 : return 0.;
59 :
60 1603 : return interpolate(z, this->redshifts, this->redshiftScalings);
61 : }
62 :
63 1 : double TabularPhotonField::getMinimumPhotonEnergy(double z) const{
64 1 : return photonEnergies[0];
65 : }
66 :
67 1 : double TabularPhotonField::getMaximumPhotonEnergy(double z) const{
68 1 : return photonEnergies[photonEnergies.size() -1];
69 : }
70 :
71 94 : void TabularPhotonField::readPhotonEnergy(std::string filePath) {
72 94 : std::ifstream infile(filePath.c_str());
73 94 : if (!infile.good())
74 0 : throw std::runtime_error("TabularPhotonField::readPhotonEnergy: could not open " + filePath);
75 :
76 : std::string line;
77 15320 : while (std::getline(infile, line)) {
78 15226 : if ((line.size() > 0) & (line[0] != '#') )
79 14944 : this->photonEnergies.push_back(std::stod(line));
80 : }
81 94 : infile.close();
82 94 : }
83 :
84 94 : void TabularPhotonField::readPhotonDensity(std::string filePath) {
85 94 : std::ifstream infile(filePath.c_str());
86 94 : if (!infile.good())
87 0 : throw std::runtime_error("TabularPhotonField::readPhotonDensity: could not open " + filePath);
88 :
89 : std::string line;
90 4772519 : while (std::getline(infile, line)) {
91 4772425 : if ((line.size() > 0) & (line[0] != '#') )
92 4772143 : this->photonDensity.push_back(std::stod(line));
93 : }
94 94 : infile.close();
95 94 : }
96 :
97 73 : void TabularPhotonField::readRedshift(std::string filePath) {
98 73 : std::ifstream infile(filePath.c_str());
99 73 : if (!infile.good())
100 0 : throw std::runtime_error("TabularPhotonField::initRedshift: could not open " + filePath);
101 :
102 : std::string line;
103 14599 : while (std::getline(infile, line)) {
104 14526 : if ((line.size() > 0) & (line[0] != '#') )
105 14307 : this->redshifts.push_back(std::stod(line));
106 : }
107 73 : infile.close();
108 73 : }
109 :
110 73 : void TabularPhotonField::initRedshiftScaling() {
111 : double n0 = 0.;
112 14380 : for (int i = 0; i < this->redshifts.size(); ++i) {
113 14307 : double z = this->redshifts[i];
114 : double n = 0.;
115 4770022 : for (int j = 0; j < this->photonEnergies.size()-1; ++j) {
116 4755715 : double e_j = this->photonEnergies[j];
117 4755715 : double e_j1 = this->photonEnergies[j+1];
118 4755715 : double deltaLogE = std::log10(e_j1) - std::log10(e_j);
119 4755715 : if (z == 0.)
120 12750 : n0 += (getPhotonDensity(e_j, 0) + getPhotonDensity(e_j1, 0)) / 2. * deltaLogE;
121 4755715 : n += (getPhotonDensity(e_j, z) + getPhotonDensity(e_j1, z)) / 2. * deltaLogE;
122 : }
123 14307 : this->redshiftScalings.push_back(n / n0);
124 : }
125 73 : }
126 :
127 94 : void TabularPhotonField::checkInputData() const {
128 94 : if (this->isRedshiftDependent) {
129 73 : if (this->photonDensity.size() != this->photonEnergies.size() * this-> redshifts.size())
130 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");
131 : } else {
132 21 : if (this->photonEnergies.size() != this->photonDensity.size())
133 0 : throw std::runtime_error("TabularPhotonField::checkInputData: length of photon energy input is unequal to length of photon density input");
134 : }
135 :
136 15038 : for (int i = 0; i < this->photonEnergies.size(); ++i) {
137 : double ePrevious = 0.;
138 14944 : double e = this->photonEnergies[i];
139 14944 : if (e <= 0.)
140 0 : throw std::runtime_error("TabularPhotonField::checkInputData: a value in the photon energy input is not positive");
141 : if (e <= ePrevious)
142 : throw std::runtime_error("TabularPhotonField::checkInputData: photon energy values are not strictly increasing");
143 : ePrevious = e;
144 : }
145 :
146 4772237 : for (int i = 0; i < this->photonDensity.size(); ++i) {
147 4772143 : if (this->photonDensity[i] < 0.)
148 0 : throw std::runtime_error("TabularPhotonField::checkInputData: a value in the photon density input is negative");
149 : }
150 :
151 94 : if (this->isRedshiftDependent) {
152 73 : if (this->redshifts[0] != 0.)
153 0 : throw std::runtime_error("TabularPhotonField::checkInputData: redshift input must start with zero");
154 :
155 14380 : for (int i = 0; i < this->redshifts.size(); ++i) {
156 : double zPrevious = -1.;
157 14307 : double z = this->redshifts[i];
158 14307 : if (z < 0.)
159 0 : throw std::runtime_error("TabularPhotonField::checkInputData: a value in the redshift input is negative");
160 : if (z <= zPrevious)
161 : throw std::runtime_error("TabularPhotonField::checkInputData: redshift values are not strictly increasing");
162 : zPrevious = z;
163 : }
164 :
165 73 : for (int i = 0; i < this->redshiftScalings.size(); ++i) {
166 0 : double scalingFactor = this->redshiftScalings[i];
167 0 : if (scalingFactor <= 0.)
168 0 : throw std::runtime_error("TabularPhotonField::checkInputData: initRedshiftScaling has created a non-positive scaling factor");
169 : }
170 : }
171 94 : }
172 :
173 40 : BlackbodyPhotonField::BlackbodyPhotonField(std::string fieldName, double blackbodyTemperature) {
174 40 : this->fieldName = fieldName;
175 40 : this->blackbodyTemperature = blackbodyTemperature;
176 40 : this->quantile = 0.0001; // tested to be sufficient, only used for extreme values of primary energy or temperature
177 40 : }
178 :
179 9532 : double BlackbodyPhotonField::getPhotonDensity(double Ephoton, double z) const {
180 9532 : return 8 * M_PI * pow_integer<3>(Ephoton / (h_planck * c_light)) / std::expm1(Ephoton / (k_boltzmann * this->blackbodyTemperature));
181 : }
182 :
183 29 : double BlackbodyPhotonField::getMinimumPhotonEnergy(double z) const {
184 : double A;
185 29 : int quantile_int = 10000 * quantile;
186 29 : switch (quantile_int)
187 : {
188 : case 1: // 0.01 % percentil
189 : A = 1.093586e-5 * eV / kelvin;
190 : break;
191 0 : case 10: // 0.1 % percentil
192 : A = 2.402189e-5 * eV / kelvin;
193 0 : break;
194 0 : case 100: // 1 % percentil
195 : A = 5.417942e-5 * eV / kelvin;
196 0 : break;
197 0 : default:
198 0 : throw std::runtime_error("Quantile not understood. Please use 0.01 (1%), 0.001 (0.1%) or 0.0001 (0.01%) \n");
199 : break;
200 : }
201 29 : return A * this -> blackbodyTemperature;
202 : }
203 :
204 29 : double BlackbodyPhotonField::getMaximumPhotonEnergy(double z) const {
205 29 : double factor = std::max(1., blackbodyTemperature / 2.73);
206 29 : return 0.1 * factor * eV; // T dependent scaling, starting at 0.1 eV as suitable for CMB
207 : }
208 :
209 0 : void BlackbodyPhotonField::setQuantile(double q) {
210 0 : if(not ((q == 0.0001) or (q == 0.001) or (q == 0.01)))
211 0 : throw std::runtime_error("Quantile not understood. Please use 0.01 (1%), 0.001 (0.1%) or 0.0001 (0.01%) \n");
212 0 : this -> quantile = q;
213 0 : }
214 :
215 : } // namespace crpropa
|
