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          436 : void Redshift::process(Candidate *c) const {
      10          436 :         double z = c->getRedshift();
      11              : 
      12              :         // check if z = 0
      13          436 :         if (z <= std::numeric_limits<double>::min())
      14            0 :                 return;
      15              : 
      16              :         // use small step approximation:  dz = H(z) / c * ds
      17          436 :         double dz = hubbleRate(z) / c_light * c->getCurrentStep();
      18              : 
      19              :         // prevent dz > z
      20          436 :         dz = std::min(dz, z);
      21              : 
      22              :         // update redshift
      23          436 :         c->setRedshift(z - dz);
      24              : 
      25              :         // adiabatic energy loss: dE / dz = E / (1 + z)
      26          436 :         double E = c->current.getEnergy();
      27          436 :         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