Line data    Source code

       1             : #include "crpropa/Cosmology.h"
       2             : #include "crpropa/Units.h"
       3             : #include "crpropa/Common.h"
       4             : 
       5             : #include <vector>
       6             : #include <cmath>
       7             : #include <stdexcept>
       8             : 
       9             : namespace crpropa {
      10             : 
      11             : /**
      12             :  @class Cosmology
      13             :  @brief Cosmology calculations
      14             :  */
      15             : struct Cosmology {
      16             :         double H0; // Hubble parameter at z=0
      17             :         double omegaM; // matter density parameter
      18             :         double omegaL; // vacuum energy parameter
      19             : 
      20             :         static const int n;
      21             :         static const double zmin;
      22             :         static const double zmax;
      23             : 
      24             :         std::vector<double> Z;  // redshift
      25             :         std::vector<double> Dc; // comoving distance [m]
      26             :         std::vector<double> Dl; // luminosity distance [m]
      27             :         std::vector<double> Dt; // light travel distance [m]
      28             : 
      29          19 :         void update() {
      30          19 :                 double dH = c_light / H0; // Hubble distance
      31             : 
      32          19 :                 std::vector<double> E(n);
      33          19 :                 E[0] = 1;
      34             : 
      35             :                 // Relation between comoving distance r and redshift z (cf. J.A. Peacock, Cosmological physics, p. 89 eq. 3.76)
      36             :                 // dr = c / H(z) dz, integration using midpoint rule
      37             :                 double dlz = log10(zmax) - log10(zmin);
      38       19000 :                 for (int i = 1; i < n; i++) {
      39       18981 :                         Z[i] = zmin * pow(10, i * dlz / (n - 1)); // logarithmic even spacing
      40       18981 :                         double dz = (Z[i] - Z[i - 1]); // redshift step
      41       18981 :                         E[i] = sqrt(omegaL + omegaM * pow_integer<3>(1 + Z[i]));
      42       18981 :                         Dc[i] = Dc[i - 1] + dH * dz * (1 / E[i] + 1 / E[i - 1]) / 2;
      43       18981 :                         Dl[i] = (1 + Z[i]) * Dc[i];
      44       18981 :                         Dt[i] = Dt[i - 1]
      45       18981 :                                         + dH * dz
      46       18981 :                                                         * (1 / ((1 + Z[i]) * E[i])
      47       18981 :                                                                         + 1 / ((1 + Z[i - 1]) * E[i - 1])) / 2;
      48             :                 }
      49          19 :         }
      50             : 
      51          19 :         Cosmology() {
      52             :         // Cosmological parameters (K.A. Olive et al. (Particle Data Group), Chin. Phys. C, 38, 090001 (2014))
      53          19 :                 H0 = 67.3 * 1000 * meter / second / Mpc; // default values
      54          19 :                 omegaM = 0.315;
      55          19 :                 omegaL = 1 - omegaM;
      56             : 
      57          19 :                 Z.resize(n);
      58          19 :                 Dc.resize(n);
      59          19 :                 Dl.resize(n);
      60          19 :                 Dt.resize(n);
      61             : 
      62          19 :                 Z[0] = 0;
      63          19 :                 Dc[0] = 0;
      64          19 :                 Dl[0] = 0;
      65          19 :                 Dt[0] = 0;
      66             : 
      67          19 :                 update();
      68          19 :         }
      69             : 
      70             :         void setParameters(double h, double oM) {
      71           0 :                 H0 = h * 1e5 / Mpc;
      72           0 :                 omegaM = oM;
      73           0 :                 omegaL = 1 - oM;
      74           0 :                 update();
      75             :         }
      76             : };
      77             : 
      78             : const int Cosmology::n = 1000;
      79             : const double Cosmology::zmin = 0.0001;
      80             : const double Cosmology::zmax = 100;
      81             : 
      82             : static Cosmology cosmology; // instance is created at runtime
      83             : 
      84           0 : void setCosmologyParameters(double h, double oM) {
      85             :         cosmology.setParameters(h, oM);
      86           0 : }
      87             : 
      88       32598 : double hubbleRate(double z) {
      89       32598 :         return cosmology.H0
      90       32598 :                         * sqrt(cosmology.omegaL + cosmology.omegaM * pow(1 + z, 3));
      91             : }
      92             : 
      93           0 : double omegaL() {
      94           0 :         return cosmology.omegaL;
      95             : }
      96             : 
      97           0 : double omegaM() {
      98           0 :         return cosmology.omegaM;
      99             : }
     100             : 
     101           0 : double H0() {
     102           0 :         return cosmology.H0;
     103             : }
     104             : 
     105           1 : double comovingDistance2Redshift(double d) {
     106           1 :         if (d < 0)
     107           0 :                 throw std::runtime_error("Cosmology: d < 0");
     108           1 :         if (d > cosmology.Dc.back())
     109           0 :                 throw std::runtime_error("Cosmology: d > dmax");
     110           1 :         return interpolate(d, cosmology.Dc, cosmology.Z);
     111             : }
     112             : 
     113           0 : double redshift2ComovingDistance(double z) {
     114           0 :         if (z < 0)
     115           0 :                 throw std::runtime_error("Cosmology: z < 0");
     116           0 :         if (z > cosmology.zmax)
     117           0 :                 throw std::runtime_error("Cosmology: z > zmax");
     118           0 :         return interpolate(z, cosmology.Z, cosmology.Dc);
     119             : }
     120             : 
     121           0 : double luminosityDistance2Redshift(double d) {
     122           0 :         if (d < 0)
     123           0 :                 throw std::runtime_error("Cosmology: d < 0");
     124           0 :         if (d > cosmology.Dl.back())
     125           0 :                 throw std::runtime_error("Cosmology: d > dmax");
     126           0 :         return interpolate(d, cosmology.Dl, cosmology.Z);
     127             : }
     128             : 
     129           0 : double redshift2LuminosityDistance(double z) {
     130           0 :         if (z < 0)
     131           0 :                 throw std::runtime_error("Cosmology: z < 0");
     132           0 :         if (z > cosmology.zmax)
     133           0 :                 throw std::runtime_error("Cosmology: z > zmax");
     134           0 :         return interpolate(z, cosmology.Z, cosmology.Dl);
     135             : }
     136             : 
     137           0 : double lightTravelDistance2Redshift(double d) {
     138           0 :         if (d < 0)
     139           0 :                 throw std::runtime_error("Cosmology: d < 0");
     140           0 :         if (d > cosmology.Dt.back())
     141           0 :                 throw std::runtime_error("Cosmology: d > dmax");
     142           0 :         return interpolate(d, cosmology.Dt, cosmology.Z);
     143             : }
     144             : 
     145           0 : double redshift2LightTravelDistance(double z) {
     146           0 :         if (z < 0)
     147           0 :                 throw std::runtime_error("Cosmology: z < 0");
     148           0 :         if (z > cosmology.zmax)
     149           0 :                 throw std::runtime_error("Cosmology: z > zmax");
     150           0 :         return interpolate(z, cosmology.Z, cosmology.Dt);
     151             : }
     152             : 
     153           2 : double comoving2LightTravelDistance(double d) {
     154           2 :         if (d < 0)
     155           0 :                 throw std::runtime_error("Cosmology: d < 0");
     156           2 :         if (d > cosmology.Dc.back())
     157           0 :                 throw std::runtime_error("Cosmology: d > dmax");
     158           2 :         return interpolate(d, cosmology.Dc, cosmology.Dt);
     159             : }
     160             : 
     161           1 : double lightTravel2ComovingDistance(double d) {
     162           1 :         if (d < 0)
     163           0 :                 throw std::runtime_error("Cosmology: d < 0");
     164           1 :         if (d > cosmology.Dt.back())
     165           0 :                 throw std::runtime_error("Cosmology: d > dmax");
     166           1 :         return interpolate(d, cosmology.Dt, cosmology.Dc);
     167             : }
     168             : 
     169             : } // namespace crpropa