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 18 : void update() {
30 18 : double dH = c_light / H0; // Hubble distance
31 :
32 18 : std::vector<double> E(n);
33 18 : 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 18000 : for (int i = 1; i < n; i++) {
39 17982 : Z[i] = zmin * pow(10, i * dlz / (n - 1)); // logarithmic even spacing
40 17982 : double dz = (Z[i] - Z[i - 1]); // redshift step
41 17982 : E[i] = sqrt(omegaL + omegaM * pow_integer<3>(1 + Z[i]));
42 17982 : Dc[i] = Dc[i - 1] + dH * dz * (1 / E[i] + 1 / E[i - 1]) / 2;
43 17982 : Dl[i] = (1 + Z[i]) * Dc[i];
44 17982 : Dt[i] = Dt[i - 1]
45 17982 : + dH * dz
46 17982 : * (1 / ((1 + Z[i]) * E[i])
47 17982 : + 1 / ((1 + Z[i - 1]) * E[i - 1])) / 2;
48 : }
49 18 : }
50 :
51 18 : Cosmology() {
52 : // Cosmological parameters (K.A. Olive et al. (Particle Data Group), Chin. Phys. C, 38, 090001 (2014))
53 18 : H0 = 67.3 * 1000 * meter / second / Mpc; // default values
54 18 : omegaM = 0.315;
55 18 : omegaL = 1 - omegaM;
56 :
57 18 : Z.resize(n);
58 18 : Dc.resize(n);
59 18 : Dl.resize(n);
60 18 : Dt.resize(n);
61 :
62 18 : Z[0] = 0;
63 18 : Dc[0] = 0;
64 18 : Dl[0] = 0;
65 18 : Dt[0] = 0;
66 :
67 18 : update();
68 18 : }
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 436 : double hubbleRate(double z) {
89 436 : return cosmology.H0
90 436 : * 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
|
