-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathdistance_functions.py
183 lines (133 loc) · 4.09 KB
/
distance_functions.py
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
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
"""
basic cosmology utility functions
"""
from __future__ import absolute_import, division, print_function, unicode_literals
import numpy as np
from astropy import units as u
from astropy import constants as const
import scipy.integrate as integrate
from .default_cosmo import default_cosmo # define a default cosology for utilities
__all__=('hubble_distance', 'comoving_distance', 'transverse_comoving_distance',
'angular_diameter_distance', 'luminosity_distance', 'comoving_volume',)
__author__=('Duncan Campbell')
def hubble_distance(H0):
"""
Calculate the Hubble distance
parameters
----------
H0: float
Hubble constant in km/s/Mpc
returns
-------
DH: float
Hubble distance in Mpc
"""
return (const.c.to('km/s')/H0).value
def comoving_distance(z,cosmo=None):
"""
Calculate the line of sight comoving distance
parameters
----------
z: float
redshift
cosmo: astropy.cosmology object, optional
cosmology object specifying cosmology. If None, FlatLambdaCDM(H0=70,Om0=0.3)
returns
-------
DC: float
Comoving line of sight distance in Mpc
"""
if cosmo==None:
cosmo = default_cosmo
f = lambda zz: 1.0/_Ez(zz, cosmo.Om0, cosmo.Ok0, cosmo.Ode0)
DC = integrate.quadrature(f,0.0,z)[0]
return hubble_distance(cosmo.H0.value)*DC
def transverse_comoving_distance(z,cosmo=None):
"""
Calculate the transverse comoving distance
parameters
----------
z: float
redshift
cosmo: astropy.cosmology object, optional
cosmology object specifying cosmology. If None, FlatLambdaCDM(H0=70,Om0=0.3)
returns
-------
DM: float
Comoving transverse distance in Mpc
"""
if cosmo==None:
cosmo = default_cosmo
if cosmo.Ok0==0.0:
return comoving_distance(z,cosmo)
elif cosmo.Ok0>0:
DC = comoving_distance(z,cosmo)
DH = hubble_distance(cosmo.H0.value)
return DH*1.0/np.sqrt(cosmo.Ok0)*np.sinh(np.sqrt(cosmo.Ok0)*DC/DH)
elif cosmo.Ok0<0:
DC = comoving_distance(z,cosmo)
DH = hubble_distance(cosmo.H0.value)
return DH*1.0/np.sqrt(np.fabs(cosmo.Ok0))*np.sin(np.sqrt(np.fabs(cosmo.Ok0))*DC/DH)
else:
raise ValueError("omega curavture value not specified.")
def angular_diameter_distance(z,cosmo=None):
"""
Calculate the angular diameter distance
parameters
----------
z: float
redshift
cosmo: astropy.cosmology object, optional
cosmology object specifying cosmology. If None, FlatLambdaCDM(H0=70,Om0=0.3)
returns
-------
DA: float
Angular diameter distance in Mpc
"""
if cosmo==None:
cosmo = default_cosmo
return transverse_comoving_distance(z,cosmo)/(1.0+z)
def luminosity_distance(z,cosmo=None):
"""
Calculate the luminosity distance.
parameters
----------
z: float
redshift
cosmo: astropy.cosmology object, optional
cosmology object specifying cosmology. If None, FlatLambdaCDM(H0=70,Om0=0.3)
returns
-------
DL: float
Luminosity distance in Mpc
"""
if cosmo==None:
cosmo = default_cosmo
return transverse_comoving_distance(z,cosmo)*(1.0+z)
def comoving_volume(z,dw,cosmo=None):
"""
Calculate comoving volume
parameters
----------
z: float
redshift
dw: float
solid angle
cosmo: astropy.cosmology object, optional
cosmology object specifying cosmology. If None, FlatLambdaCDM(H0=70,Om0=0.3)
returns
-------
VC: float
comoving volume in Mpc^3
"""
if cosmo==None:
cosmo = default_cosmo
DH = hubble_distance(cosmo.H0.value)
f = lambda zz: DH*((1.0+zz)**2.0*angular_diameter_distance(zz,cosmo)**2.0)/(_Ez(zz, cosmo.Om0, cosmo.Ok0, cosmo.Ode0))
VC = integrate.quadrature(f, 0.0, z, vec_func=False)[0]*dw
return VC
def _Ez(z, omega_m, omega_k, omega_l):
"""
internal function used for distance calculations
"""
return np.sqrt(omega_m*(1.0+z)**3.0+omega_k*(1.0+z)**2.0+omega_l)