Source Code for Module pylal.cosmography

 1  # Copyright (C) 2011 Nickolas Fotopoulos 
 2  # 
 3  # This program is free software; you can redistribute it and/or modify it 
 4  # under the terms of the GNU General Public License as published by the 
 5  # Free Software Foundation; either version 2 of the License, or (at your 
 6  # option) any later version. 
 7  # 
 8  # This program is distributed in the hope that it will be useful, but 
 9  # WITHOUT ANY WARRANTY; without even the implied warranty of 
10  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General 
11  # Public License for more details. 
12  # 
13  # You should have received a copy of the GNU General Public License along 
14  # with this program; if not, write to the Free Software Foundation, Inc., 
15  # 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA. 
16  """ 
17  Cosmography based on Hogg 1998 (astro-ph/9905116) 
18  Cosmological parameters taken from WMAP 7 year BAO + H0 ML value (arXiv:1001.4538) 
19  """ 
20   
21  from numpy import pi, sinh, arcsinh 
22  from scipy.integrate import quad 
23  from scipy.optimize import newton 
24   
25  from lal import PC_SI as LAL_PC_SI 
26  from lal import C_SI as LAL_C_SI 
27   
28  h0 = 0.703  # H0 = h0 * 100 km / s / Mpc 
29  H0_SI = h0 * 100000 / (1000000 * LAL_PC_SI)  # Hubble constant in inverse seconds 
30  DH_SI = LAL_C_SI / H0_SI  # Hubble distance in meters 
31  OmegaL = 0.729 
32  OmegaM = 0.271 # sum of baryonic and cold dark matter Omega 
33  OmegaR = 0 # the remainder is consistent with zero given the error bars 
34   
35 -def Einv(z): 
36      """ 
37      Integrand used in many cosmography integrals. 1 / E(z) in Hogg's notation. 
38      """ 
39      return (OmegaM * (1 + z)**3 + OmegaR * (1 + z)**2 + OmegaL)**-0.5 
40   
41 -def DM(z): 
42      """ 
43      Comoving distance (transverse) at z. Hard-coded for OmegaR = 0. 
44      """ 
45      return DH_SI * quad(Einv, 0, z)[0] 
46   
47 -def DL(z): 
48      """ 
49      Luminosity distance 
50      """ 
51      return DM(z) * (1 + z) 
52   
53 -def compute_redshift(DL0): 
54      """ 
55      Compute the redshift from a luminosity distance. 
56      Use the Newton-Raphson refinement method on a first-order guess. 
57      """ 
58      return newton(lambda z: DL(z) - DL0, DL0 / DH_SI) 
59   
60 -def dVdz(z): 
61      """ 
62      Different volume element per unit redshift. Hard-coded for OmegaR = 0. 
63      """ 
64      return DH_SI * DM(z)**2 * Einv(z) * 4 * pi 
65   
66 -def V(z): 
67      """ 
68      Analytical integration of dVdz. Hard-coded for OmegaR = 0. 
69      Double precision craps out below about 100 kpc (z=2.3e-6). 
70      """ 
71      tmpDM = DM(z) 
72      return 4 * pi / 3 * tmpDM * tmpDM * tmpDM 
73   
74 -def V_from_DL_z(D, z=None): 
75      """ 
76      Analytical integration of dVdz. Hard-coded for OmegaR = 0. 
77      Sped up for the case that you have D rather than (or in addition to) z. 
78      Double precision craps out below about 100 kpc (z=2.3e-6). 
79      """ 
80      tmpDM = D / (1 + (z or compute_redshift(D))) 
81      return 4 * pi / 3 * tmpDM * tmpDM * tmpDM 
82