Line data    Source code

       1             : #include "crpropa/module/Redshift.h"
       2             : #include "crpropa/Units.h"
       3             : #include "crpropa/Cosmology.h"
       4             : 
       5             : #include <limits>
       6             : 
       7             : namespace crpropa {
       8             : 
       9       32602 : void Redshift::process(Candidate *c) const {
      10       32602 :         double z = c->getRedshift();
      11             : 
      12             :         // check if z = 0
      13       32602 :         if (z <= std::numeric_limits<double>::min())
      14           4 :                 return;
      15             : 
      16             :         // use small step approximation:  dz = H(z) / c * ds
      17       32598 :         double dz = hubbleRate(z) / c_light * c->getCurrentStep();
      18             : 
      19             :         // prevent dz > z
      20       32598 :         dz = std::min(dz, z);
      21             : 
      22             :         // update redshift
      23       32598 :         c->setRedshift(z - dz);
      24             : 
      25             :         // adiabatic energy loss: dE / dz = E / (1 + z)
      26       32598 :         double E = c->current.getEnergy();
      27       32598 :         c->current.setEnergy(E * (1 - dz / (1 + z)));
      28             : }
      29             : 
      30           0 : std::string Redshift::getDescription() const {
      31           0 :         std::stringstream s;
      32           0 :         s << "Redshift: h0 = " << hubbleRate() / 1e5 * Mpc << ", omegaL = "
      33           0 :                         << omegaL() << ", omegaM = " << omegaM();
      34           0 :         return s.str();
      35           0 : }
      36             : 
      37           0 : void FutureRedshift::process(Candidate *c) const {
      38           0 :         double z = c->getRedshift();
      39             : 
      40             :         // check if z = -1
      41           0 :         if (z <= -1)
      42             :                 return;
      43             : 
      44             :         // use small step approximation:  dz = H(z) / c * ds
      45           0 :         double dz = hubbleRate(z) / c_light * c->getCurrentStep();
      46             : 
      47             :         // update redshift
      48           0 :         c->setRedshift(z - dz);
      49             : 
      50             :         // adiabatic energy loss: dE / dz = E / (1 + z)
      51           0 :         double E = c->current.getEnergy();
      52           0 :         c->current.setEnergy(E * (1 - dz / (1 + z)));
      53             : }
      54             : 
      55           0 : std::string FutureRedshift::getDescription() const {
      56           0 :         std::stringstream s;
      57           0 :         s << "FutureRedshift: h0 = " << hubbleRate() / 1e5 * Mpc << ", omegaL = "
      58           0 :                         << omegaL() << ", omegaM = " << omegaM();
      59           0 :         return s.str();
      60           0 : }
      61             : 
      62             : } // namespace crpropa