blob: 4c96dc8656b3279043f198edb394c11efae2db9f (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
|
// calc_distance.cc
// calculate commonly used cosmology distances
// Ref: https://en.wikipedia.org/wiki/Distance_measures_(cosmology)
//
// Junhua Gu
//
// ChangeLogs:
// 2012/08/12, LIweitiaNux
// fix a bug in `calc_angular_distance()'
// add `calc_transcomv_distance()' and
// account `omega_k >0, <0' cases
// add `calc_ligtrav_distance()'
//
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <cstddef>
#include <cassert>
// integrator
#include "adapt_trapezoid.h"
// extern variables in `calc_distance.h'
extern double c;
extern double km;
extern double s;
extern double Mpc;
// cosmology parameters
extern double H0; // units: [km/s/Mpc]
extern double omega_m;
extern double omega_l;
extern double omega_k;
// precision to determine float number equality
extern double eps;
using namespace std;
// auxiliary functions
// for `calc_comoving_distance'
double f_comvdist(double z);
// for `calc_ligtrav_distance'
double f_ligdist(double z);
double f_age(double z);
/////////////////////////////////////////////////
// main functions
// //////////////////////////////////////////////
// dimensionless Hubble parameter
double E(double z)
{
return sqrt(omega_m*(1+z)*(1+z)*(1+z) + omega_k*(1+z)*(1+z) + omega_l);
}
// Hubble distance
double calc_hubble_distance(void)
{
//
// cerr << "*** H0 = " << H0 << " ***" << endl;
//
double _H0_ = H0 * km/s/Mpc;
return c / _H0_;
}
// calc `comoving distance'
double calc_comoving_distance(double z)
{
double d_H = calc_hubble_distance();
return d_H * adapt_trapezoid(f_comvdist, 0.0, z, 1e-4);
}
// calc `transverse comoving distance'
double calc_transcomv_distance(double z)
{
if (fabs(omega_k) <= eps) {
// omega_k = 0
// flat, d_M(z) = d_C(z)
return calc_comoving_distance(z);
}
else if (omega_k > 0.0) {
// omega_k > 0
// d_M(z) = d_H/sqrt(omega_k) * sinh(sqrt(omega_k) * d_C(z) / d_H)
// ^ ^
double d_H = calc_hubble_distance();
double d_C = calc_comoving_distance(z);
return (d_H / sqrt(omega_k) * sinh(sqrt(omega_k) * d_C / d_H));
}
else {
// omega_k < 0
// d_M(z) = d_H/sqrt(-omega_k) * sin(sqrt(-omega_k) * d_C(z) / d_H)
// ^ ^ ^
double d_H = calc_hubble_distance();
double d_C = calc_comoving_distance(z);
return (d_H / sqrt(-omega_k) * sin(sqrt(-omega_k) * d_C / d_H));
}
}
// calc `angular diameter distance'
// d_A(z) = d_M(z)/(1+z)
double calc_angdia_distance(double z)
{
return (calc_transcomv_distance(z) / (1+z));
}
// calc `luminoisity distance'
// d_L(z) = (1+z)*d_M(z)
double calc_luminosity_distance(double z)
{
return (calc_transcomv_distance(z) * (1+z));
}
// calc `light-travel distance'
// d_T(z) = d_H \int_0^z 1/((1+z)*E(z)) dz
double calc_ligtrav_distance(double z)
{
double d_H = calc_hubble_distance();
return d_H * adapt_trapezoid(f_ligdist, 0.0, z, 1e-4);
}
// auxiliary functions
// for `calc_comoving_distance'
double f_comvdist(double z)
{
return 1.0/E(z);
}
// for `calc_ligtrav_distance'
double f_ligdist(double z)
{
return 1.0 / ((1+z)*E(z));
}
double f_age(double z)
{
return f_comvdist(1.0/z) / (z*z);
}
// EOF
/* vim: set ts=8 sw=4 tw=0 fenc=utf-8 ft=cpp: */
|
