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         188 :         readPhotonEnergy(getDataPath("") + "Scaling/" + this->fieldName + "_photonEnergy.txt");
      19         188 :         readPhotonDensity(getDataPath("") + "Scaling/" + this->fieldName + "_photonDensity.txt");
      20          94 :         if (this->isRedshiftDependent)
      21         146 :                 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     9586406 : double TabularPhotonField::getPhotonDensity(double Ephoton, double z) const {   
      31     9586406 :         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     9586406 :                 double zMin = this->redshifts[0];
      36     9586406 :                 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     9586406 :                         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       66390 : double TabularPhotonField::getRedshiftScaling(double z) const {
      51       66390 :         if (!this->isRedshiftDependent)
      52             :                 return 1.;
      53             :  
      54       66206 :         if (z < this->redshifts.front())
      55             :                 return 1.;
      56             :  
      57       66206 :         if (z > this->redshifts.back())
      58             :                 return 0.;
      59             :  
      60       66206 :         return interpolate(z, this->redshifts, this->redshiftScalings);
      61             : }
      62             : 
      63          22 : double TabularPhotonField::getMinimumPhotonEnergy(double z) const{
      64          22 :         return photonEnergies[0];
      65             : }
      66             : 
      67          22 : double TabularPhotonField::getMaximumPhotonEnergy(double z) const{
      68          22 :         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       14307 :                         if (z <= zPrevious)
     161           0 :                                 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       14419 : double BlackbodyPhotonField::getPhotonDensity(double Ephoton, double z) const {
     180       14419 :         return 8 * M_PI * pow_integer<3>(Ephoton / (h_planck * c_light)) / std::expm1(Ephoton / (k_boltzmann * this->blackbodyTemperature));
     181             : }
     182             : 
     183          44 : double BlackbodyPhotonField::getMinimumPhotonEnergy(double z) const {
     184             :         double A;
     185          44 :         int quantile_int = 10000 * quantile;
     186          44 :         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          44 :         return A * this -> blackbodyTemperature;
     202             : }
     203             : 
     204          44 : double BlackbodyPhotonField::getMaximumPhotonEnergy(double z) const {
     205          44 :         double factor = std::max(1., blackbodyTemperature / 2.73);
     206          44 :         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