1616`projected_P_BBPS_real` interpolates over a table of comoving distances to
1717obtain a more precise answer.
1818'''
19- import astropy .constants as c
20- import astropy .units as u
2119import numpy as np
2220from scipy .integrate import quad
2321
2422
2523def _rho_crit (z , omega_m , omega_lambda ):
2624 '''
27- The critical density of the universe, in units of $Msun*Mpc^{-3}*h^{-2} $.
25+ The critical density of the universe, in units of $Msun*Mpc^{-3}*h^2 $.
2826 '''
2927 # The constant is 3 * (100 km / s / Mpc)**2 / (8 * pi * G)
3028 # in units of Msun h^2 Mpc^{-3}
29+ # (source: astropy's constants module and unit conversions)
3130 return 2.77536627e+11 * (omega_m * (1 + z )** 3 + omega_lambda )
3231
3332
@@ -41,7 +40,9 @@ def P_delta(M, z, omega_b, omega_m, omega_lambda, delta=200):
4140 See BBPS, section 4.1 for details.
4241 Units: Msun s^{-2} Mpc^{-1}
4342 '''
44- return c .G * M * delta * _rho_crit (z , omega_m , omega_lambda ) * \
43+ # G = 4.51710305e-48 Mpc^3 Msun^{-1} s^{-2}
44+ # (source: astropy's constants module and unit conversions)
45+ return 4.51710305e-48 * M * delta * _rho_crit (z , omega_m , omega_lambda ) * \
4546 omega_b / omega_m / 2 / R_delta (M , omega_m , omega_lambda , z , delta )
4647
4748
@@ -151,7 +152,7 @@ def projected_P_BBPS(r, M, z, omega_b, omega_m, omega_lambda,
151152 R_del = R_delta (M , z , omega_m , omega_lambda )
152153 return quad (lambda x : P_BBPS (np .sqrt (x * x + r * r ), M , z ,
153154 omega_b , omega_m ,
154- omega_lambda ). value ,
155+ omega_lambda ),
155156 - dist * R_del , dist * R_del ,
156157 epsrel = 1e-3 )[0 ] / (1 + z )
157158
@@ -181,8 +182,7 @@ def projected_P_BBPS_real(r, M, z, omega_b, omega_m, omega_lambda, chis, zs,
181182 return quad (lambda x : P_BBPS (np .sqrt ((x - chi_cluster )** 2 + r * r ),
182183 M , z ,
183184 omega_b , omega_m ,
184- omega_lambda ).value
185- / (1 + np .interp (x , chis , zs )),
185+ omega_lambda ) / (1 + np .interp (x , chis , zs )),
186186 chi_cluster - dist * R_del ,
187187 chi_cluster + dist * R_del ,
188188 epsrel = 1e-3 )[0 ]
0 commit comments