r""" A general FLRW cosmology module. Assumes a spatially flat universe. """
import numpy as np
from scipy.integrate import quad
from scipy.interpolate import interp1d
from rabacus.constants import physical
from rabacus.constants import units
from rabacus.utils import utils
__all__ = ['Cosmology']
class Error(Exception):
r"""Base class for exceptions in this module."""
pass
class InitError(Error):
r"""Exception raised for not having required keys in input. """
def __init__(self, message ):
message = message + '\n' + \
'cpdict must contain the following keys ... \n' + \
'"omegam" -> current matter density in units of critical \n' + \
'"omegal" -> current lambda density in units of critical \n' + \
'"omegab" -> current baryon density in units of critical \n' + \
'"h" -> Hubble parameter H0 = 100 h km/s/Mpc \n' + \
'"sigma8 -> density variance in spheres w/ R = 8 Mpc/h \n' + \
'"ns" -> slope of primordial power spectrum \n' + \
'"Yp" -> primordial mass fraction of helium \n'
Exception.__init__( self, message )
[docs]class Cosmology(object):
r""" General Cosmology Class. Assumes a spatialy flat universe.
Args:
`cpdict` (dict) A dictionary of cosmological parameters.
For example, see
:class:`~rabacus.cosmology.parameters.planck.load.PlanckParameters`.
The dictionray `cpdict` must include the following keys,
- ``omegam`` -> current matter density in units of critical today
- ``omegal`` -> current lambda density in units of critical today
- ``omegab`` -> current baryon density in units of critical today
- ``h`` -> Hubble parameter H0 = 100 h km/s/Mpc
- ``sigma8`` -> amplitude of fluctuations in spheres w/ R = 8 Mpc/h
- ``ns`` -> slope of primordial power spectrum
- ``Yp`` -> primordial mass fraction of helium
Attributes:
H0 (real): hubble parameter now, :math:`H_0`.
OmegaB (real): baryon density / critical density now,
:math:`\Omega_b`
OmegaC (real): cold dark matter density / critical density now,
:math:`\Omega_c`
OmegaL (real): dark energy density / critical density now,
:math:`\Omega_{\Lambda}`
OmegaM (real): matter density / critical density now,
:math:`\Omega_m`
Yp (real): primoridial helium mass fraction,
:math:`Y_p`
cu (:class:`~rabacus.constants.units.CosmoUnits`): cosmological
units which are aware of the hubble parameter.
dH0 (real): hubble distance now,
:math:`d_{\rm H_0} = c / H_0`
tH0 (real): hubble time now,
:math:`t_{\rm H_0} = 1 / H_0`
rho_crit0 (real): critical density now,
:math:`\rho_{c,0} = (3 H_0^2) / (8 \pi G)`
eps_crit0 (real): critical energy density now,
:math:`\epsilon_{c,0} = \rho_{c,0} \, c^2`
nH_crit0 (real): critical hydrogen number density now,
:math:`n_{\rm H,c,0} = \rho_{c,0} \, \Omega_b \, (1-Y_p) \, / \,
m_{\rm H}`
"""
def __init__( self, cpdict, verbose=False, zlo=0.0, zhi=200.0, Nz=500 ):
# check input dictionary
#--------------------------------------------------------
need_keys = ['omegam', 'omegal', 'omegab', 'h', 'sigma8', 'ns', 'Yp']
for key in need_keys:
if not cpdict.has_key( key ):
raise InitError, ' key missing from cpdict \n'
# transfer input parametrs
#--------------------------------------------------------
self.OmegaM = cpdict['omegam']
self.OmegaL = cpdict['omegal']
self.OmegaB = cpdict['omegab']
self.h = cpdict['h']
self.sigma8 = cpdict['sigma8']
self.ns = cpdict['ns']
self.Yp = cpdict['Yp']
self.OmegaC = self.OmegaM - self.OmegaB
if cpdict.has_key('omegar'):
self.OmegaR = cpdict['omegar']
else:
# if not provided we use OmegaM and zeq from Planck13
# column one of table 2 in http://arxiv.org/abs/1303.5076
# OmegaR = OmegaM / (1+zeq), OmegaM=0.3175, zeq = 3402
self.OmegaR = 0.3175 / (1.0+3402.)
self.cpdict = cpdict
# make sure flat universe is specified
#--------------------------------------------------------
err = ( self.OmegaR + self.OmegaM + self.OmegaL ) - 1.0
if err > 1.0e-3:
raise InitError, ' OmegaR + OmegaM + OmegaL must be 1.0 \n' + \
' OmegaR: ' + str(self.OmegaR) + '\n' + \
' OmegaM: ' + str(self.OmegaM) + '\n' + \
' OmegaL: ' + str(self.OmegaL) + '\n' + \
' sum: ' + str(self.OmegaR + self.OmegaM + self.OmegaL) + '\n'
# create physical constants and cosmological units
#--------------------------------------------------------
self.pc = physical.PhysicalConstants()
self.cu = units.CosmoUnits( self.h, a=1.0 )
# derive constants
#--------------------------------------------------------
# Hubble parameter now
self.H0 = 1.0e2 * self.cu.h * self.cu.km / self.cu.s / self.cu.Mpc
self.H0.units = self.cu.km / self.cu.s / (self.cu.Mpc / self.cu.h)
# Hubble time now
self.tH0 = 1.0 / self.H0
self.tH0.units = self.cu.Myr/self.cu.h
# Hubble distance now
self.dH0 = self.pc.c / self.H0
self.dH0.units = self.cu.Mpc/self.cu.h
# Critical mass density now
self.rho_crit0 = ( 3.0 * self.H0**2 ) / (8.0 * np.pi * self.pc.G)
self.rho_crit0.units = (self.cu.Msun/self.cu.h) / \
(self.cu.Mpc/self.cu.h)**3
# Critical energy density now
self.eps_crit0 = self.rho_crit0 * self.pc.c**2
self.eps_crit0.units = self.cu.erg / self.cu.cm**3
# Critical hydrogen number density now
H_rho = self.rho_crit0 * self.OmegaB * (1.0-self.Yp)
self.nH_crit0 = H_rho / ( self.pc.m_p + self.pc.m_e )
self.nH_crit0.units = self.cu.cm**(-3)
# normalize the growth functions such that ...
# D1a(a=1) = 1.0
# D1z(z=0) = 1.0
#---------------------------------------------------------
self._D1a_Norm = 1.0 / quad( self._D1a_integrand, 0.0, 1.0 )[0]
self._D1z_Norm = 1.0 / quad( self._D1z_integrand, 0.0, np.inf )[0]
# make interpolating functions for ....
# redshift and comoving distance,
# redshift and lookback time
# for redshift we use log(1+z)
#--------------------------------------------------------
# first tabulate values with respect to log(1+z)
# then create interpolating functions
self._lopz = np.linspace( np.log10(1+zlo),
np.log10(1+zhi),
Nz )
self._z = 10**self._lopz - 1
self._a = 1.0 / (1.0 + self._z)
self._Dc = self.Dcz( self._z )
self._Dc.units = 'Mpc/hh'
self._Dc2lopz = interp1d( self._Dc.magnitude, self._lopz )
self._tL = self.tLz( self._z )
self._tL.units = 'Myr/hh'
self._tL2lopz = interp1d( self._tL.magnitude, self._lopz )
[docs] def z_at_Dc(self, Dc):
r""" Redshift at a given comoving distance, z(Dc). Inverts the
function :func:`Dcz` by interpolating between tabulated values.
.. math::
z(D_C) = {\rm Inverse}[ D_C(z) ]
"""
if hasattr(Dc,'units'):
Dc.units = 'Mpc/hh'
else:
msg = '\n Input variable Dc must have units of length. \n'
raise ValueError(msg)
lopz = self._Dc2lopz( Dc )
z = (10**lopz-1)
return z
[docs] def z_at_tL(self, tL):
r""" Redshift at a given lookback time, z(tL). Inverts the
function :func:`tLz` by interpolating between tabulated values.
.. math::
z(t_L) = {\rm Inverse}[ t_L(z) ]
"""
if hasattr(tL,'units'):
tL.units = 'Myr/hh'
else:
msg = '\n Input variable tL must have units of time. \n'
raise ValueError(msg)
lopz = self._tL2lopz( tL )
z = (10**lopz-1)
return z
[docs] def tL_at_Dc(self, Dc):
r""" Lookback time at a given comoving distance, tL(Dc). First
finds a redshift `z` by inverting the function :func:`Dcz` using
tabulated values. Second, calls :func:`tLz` to get a lookback
time from `z`.
.. math::
z = {\rm Inverse}[ D_C(z) ] \\
t_L(z) = \int_0^z \frac{dz'}{(1+z') H(z')}
"""
if hasattr(Dc,'units'):
Dc.units = 'Mpc/hh'
else:
msg = '\n Input variable Dc must have units of length. \n'
raise ValueError(msg)
lopz = self._Dc2lopz( Dc )
z = (10**lopz-1)
tL = self.tLz(z)
return tL
[docs] def Dc_at_tL(self, tL):
r""" Comoving distance at a given lookback time, Dc(tL). First
finds a redshift `z` by inverting the function :func:`tLz` using
tabulated values. Second, calls :func:`Dcz` to get a comoving
distance from `z`.
.. math::
z = {\rm Inverse}[ t_L(z) ] \\
D_C(z) = d_{\rm H_0} \int_{0}^{z} \frac{dz'}{E(z')}
"""
if hasattr(tL,'units'):
tL.units = 'Myr/hh'
else:
raise ValueError('\n Input variable tL must have units \n')
lopz = self._tL2lopz( tL )
z = (10**lopz-1)
Dc = self.Dcz(z)
return Dc
[docs] def Ea(self, a):
r""" Helper function H(a) = H0 * E(a)
.. math::
E(a) = ( \Omega_R a^{-4} + \Omega_M a^{-3} + \Omega_{\Lambda} )^{1/2}
"""
E = np.sqrt( self.OmegaR*a**(-4) + self.OmegaM*a**(-3) + self.OmegaL )
return E
[docs] def Ez(self, z):
r""" Helper function H(z) = H0 * E(z)
.. math::
E(z) = [ \Omega_R (1+z)^{4} + \Omega_M (1+z)^{3} +
\Omega_{\Lambda} ]^{1/2}
"""
E = np.sqrt( self.OmegaR*(1+z)**4 + self.OmegaM*(1+z)**3 + self.OmegaL )
return E
[docs] def iEa(self, a):
r""" Reciprocal of :func:`Ea`.
.. math::
E(a) = ( \Omega_R a^{-4} + \Omega_M a^{-3} +
\Omega_{\Lambda} )^{-1/2}
"""
iE = 1.0 / self.Ea(a)
return iE
[docs] def iEz(self, z):
r""" Reciprocal of :func:`Ez`.
.. math::
E(z) = [ \Omega_R (1+z)^{4} + \Omega_M (1+z)^{3} +
\Omega_{\Lambda} ]^{-1/2}
"""
iE = 1.0 / self.Ez(z)
return iE
[docs] def Ha(self, a):
r""" Hubble parameter at scale factor `a`, H(a) = H0 * E(a)
.. math::
H(a) = H_0 E(a)
"""
H = self.H0 * self.Ea(a)
return H
[docs] def Hz(self, z):
r""" Hubble parameter at redshift `z`, H(z) = H0 * E(z)
.. math::
H(z) = H_0 E(z)
"""
H = self.H0 * self.Ez(z)
return H
[docs] def tHa(self, a):
r""" Hubble time at scale factor `a`, tH(a)
.. math::
t_{\rm H(a)} = 1 / H(a)
"""
t = 1.0 / self.Ha(a)
t.units = 'Myr/hh'
return t
[docs] def tHz(self, z):
r""" Hubble time at redshift `z`, tH(z)
.. math::
t_{\rm H(z)} = 1 / H(z)
"""
t = 1.0 / self.Hz(z)
t.units = 'Myr/hh'
return t
[docs] def dHa(self, a):
r""" Hubble distance at scale factor `a`, dH(a)
.. math::
d_{\rm H(a)} = c / H(a)
"""
d = self.pc.c / self.Ha(a)
d.units = 'Mpc/hh'
return d
[docs] def dHz(self, z):
r""" Hubble distance at redshift `z`, dH(z)
.. math::
d_{\rm H(z)} = c / H(z)
"""
d = self.pc.c / self.Hz(z)
d.units = 'Mpc/hh'
return d
[docs] def rho_crita(self, a):
r""" Critical mass density at scale factor `a`, rho_crit(a)
.. math::
\rho_{c,a} = \rho_{c,0} \, E(a)^2
"""
rho = self.rho_crit0 * self.Ea(a)**2
return rho
[docs] def rho_critz(self, z):
r""" Critical mass density at redshift `z`, rho_crit(z)
.. math::
\rho_{c,z} = \rho_{c,0} \, E(z)^2
"""
rho = self.rho_crit0 * self.Ez(z)**2
return rho
[docs] def eps_crita(self, a):
r""" Critical energy density at scale factor `a`, eps_crit(a)
.. math::
\epsilon_{c,a} = \epsilon_{c,0} \, E(a)^2
"""
eps = self.eps_crit0 * self.Ea(a)**2
return eps
[docs] def eps_critz(self, z):
r""" Critical energy density at redshift `z`, eps_crit(z)
.. math::
\epsilon_{c,z} = \epsilon_{c,0} \, E(z)^2
"""
eps = self.eps_crit0 * self.Ez(z)**2
return eps
[docs] def nH_crita(self, a):
r""" Critical hydrogen number density at scale factor `a`, nH_crit(a)
.. math::
n_{\rm H,c,a} = n_{\rm H,c,0} \, E(a)^2
"""
nH = self.nH_crit0 * self.Ea(a)**2
return nH
[docs] def nH_critz(self, z):
r""" Critical hydrogen number density at redshift `z`, nH_crit(z)
.. math::
n_{\rm H,c,z} = n_{\rm H,c,0} \, E(z)^2
"""
nH = self.nH_crit0 * self.Ez(z)**2
return nH
[docs] def ta( self, a):
r""" Time since a=0 or age of the Universe, t(a).
.. math::
t(a) = \int_0^a \frac{da'}{a' H(a')}
"""
# quad works but takes a long time.
#-------------------------------------------------------------------
# if utils.isiterable( a ):
# t = np.array( [quad( self._dt_da, 0.0, aa )[0] for aa in a] )
# else:
# t = quad( self._dt_da, 0.0, a )[0]
# t = t * self.cu.s
# print 't = ', t.rescale("Myr/hh")
# trap rule includes OmegaR and is fast and accurate to 6 decimals
#-------------------------------------------------------------------
N = 5000
a_min = 1.0e-10
if utils.isiterable( a ):
t = np.zeros( a.size ) * self.cu.s
for ii,aa in enumerate(a):
xx = np.linspace( a_min, aa, N )
yy = 1.0 / ( xx * self.Ha(xx) )
t[ii] = utils.trap( xx, yy )
else:
xx = np.linspace( a_min, a, N )
yy = 1.0 / ( xx * self.Ha(xx) )
t = utils.trap( xx, yy )
t.units = 'Myr/hh'
# this is analytic but does not include OmegaR
#-------------------------------------------------------------------
# pre = 2.0 / ( 3.0 * np.sqrt( self.OmegaL ) )
# aeq = (self.OmegaM/self.OmegaL)**(1./3.)
# arg = (a/aeq)**(3./2.) + np.sqrt(1 + (a/aeq)**3)
# t = pre * np.log(arg) * self.tH0
return t
def _dt_da( self, a ):
r""" Differential dt/da """
Ha = self.Ha(a).rescale('1/s')
Ha = Ha.magnitude
dtda = 1.0 / (a * Ha)
return dtda
[docs] def tz( self, z):
r""" Time since z=inf or age of the Universe t(z).
.. math::
t(z) = \int_z^\infty \frac{dz'}{(1+z') H(z')}
"""
a = 1.0/(1.0+z)
t = self.ta(a)
return t
[docs] def tLa( self, a):
r""" Lookback from a = 1 to a, tL(a).
.. math::
t_{\rm L}(a) = \int_a^1 \frac{da'}{a' H(a')}
"""
t1 = self.ta(1.0)
tL = t1 - self.ta(a)
return tL
[docs] def tLz( self, z):
r""" Lookback from z = 0 to z, tL(z).
.. math::
t_{\rm L}(z) = \int_0^z \frac{dz'}{(1+z') H(z')}
"""
t1 = self.tz(0.0)
tL = t1 - self.tz(z)
return tL
[docs] def Dcz(self,z):
r""" Comoving distance between z=0 and z, Dc(z)
.. math::
D_C(z) = d_{\rm H_0} \int_{0}^{z} \frac{dz'}{E(z')}
"""
if utils.isiterable( z ):
Dc = np.array( [quad( self.iEz, 0.0, zz )[0] for zz in z] )
Dc = Dc * self.dH0
else:
Dc = quad( self.iEz, 0.0, z )[0] * self.dH0
Dc.units = self.cu.Mpc/self.cu.h
return Dc
[docs] def Dca(self,a):
r""" Comoving distance between a=1 and a, Dc(a). The scale factor
`a` is converted to redshift `z` and then :func:`Dcz` is called.
.. math::
D_C(a) = d_{\rm H_0} \int_{a}^{1} \frac{da'}{a'^2 E(a')}
"""
z = 1.0/a - 1.0
Dc = self.Dcz(z)
return Dc
[docs] def dz2dDc(self,z1,z2):
r""" Comoving distance between `z1` and `z2`. Two calls to :func:`Dcz`
are made and the results subtracted. `z1` < `z2` produces positive
comoving distance.
.. math::
D_C(z1,z2) = d_{\rm H_0} \int_{z1}^{z2} \frac{dz'}{E(z')}
"""
if utils.isiterable(z1) and utils.isiterable(z2):
if not len(z1) == len(z2):
raise utils.InputError, "len(z1) /= len(z2)"
Dc1 = self.Dcz(z1)
Dc2 = self.Dcz(z2)
Dc = Dc2 - Dc1
return Dc
[docs] def da2dDc(self,a1,a2):
r""" Comoving distance between `a1` and `a2`. The scale factors
`a1` and `a2` are converted to redshifts `z1` and `z2` and then
:func:`dz2dDc` is called. `a1` < `a2` produces positive comoving
distance.
.. math::
D_C(a1,a2) = d_{\rm H_0} \int_{a1}^{a2} \frac{da'}{a'^2 E(a')}
"""
z2 = 1.0/a1 - 1.0
z1 = 1.0/a2 - 1.0
Dc = self.dz2dDc(z1,z2)
return Dc
[docs] def da2dtL(self,a1,a2):
r""" Lookback time between `a1` and `a2`. Two calls to :func:`tLa`
are made and the results subtracted. `a1` < `a2` produces positive
lookback time.
.. math::
t_{\rm L}(a1,a2) = \int_{a1}^{a2} \frac{da'}{a' H(a')}
"""
t1 = self.ta( a1 )
t2 = self.ta( a2 )
dt = t2 - t1
return dt
[docs] def dz2dtL(self,z1,z2):
r""" Lookback time between `z1` and `z2`. The redshifts `z1` and
`z2` are converted to scale factors and then :func:`da2dtL` is called.
`z1` < `z2` produces positive lookback time.
.. math::
t_{\rm L}(z1,z2) = \int_{z1}^{z2} \frac{dz'}{(1+z') H(z')}
"""
a2 = 1.0 / (1.0 + z1)
a1 = 1.0 / (1.0 + z2)
dt = self.da2dtL( a1, a2 )
return dt
[docs] def X(self,z1,z2):
r""" Absorption distance between `z1` and `z2`,
.. math::
X(z1,z2) = \int_{z1}^{z2} dz' (1+z')^2 / E(z')
"""
X = quad( self._X_integrand, z1, z2 )
return X[0]
def _X_integrand( self, z ):
r""" function to interate in order to find absorption distance. """
dXdz = (1.0+z)*(1.0+z) * self.iEz(z)
return dXdz
# Cosmological growth functions
#----------------------------------------------------------------
[docs] def D1z( self, z ):
r""" Linear growth function D1(z)
.. math::
D_1(z) \propto H(z) \int_z^\infty \frac{1+z'}{H(z')^3} dz'
The normalization is such that D1z(0) = 1
"""
if utils.isiterable(z):
int_u = [quad( self._D1z_integrand, zz, np.inf )[0] for zz in z]
else:
int_u = quad( self._D1z_integrand, z, np.inf )[0]
D1 = self.Hz(z) / self.H0 * int_u * self._D1z_Norm
return D1
def _D1z_integrand( self, z ):
r""" Integrand for D1z """
return (1.0 + z) / self.Hz(z)**3
[docs] def D1a(self, a):
r""" Linear growth function D1(a).
.. math::
D_1(a) \propto H(a) \int_0^a \frac{da'}{a'^3 H(a')^3}
The normalization is such that D1a(1) = 1
"""
if utils.isiterable(a):
integral = [quad( self._D1a_integrand, 0.0, aa )[0] for aa in a]
else:
integral = quad( self._D1a_integrand, 0.0, a )[0]
D1 = self.Ha(a) / self.H0 * integral * self._D1a_Norm
return D1
def _D1a_integrand( self, a ):
r""" Integrand for D1a """
return 1.0 / ( a**3 * self.Ha(a)**3 )
# Integrands
#----------------------------------------------------------------
# plotting
#----------------------------------------------------------------
# def _plot_D1a(self):
# N = 100
# aa = np.linspace( 1.0e-4, 1.0, N )
# zz = 1.0/aa - 1.0
# D1z = self.D1z(zz)
# D1a = self.D1a(aa)
#
# fig = plt.figure()
# plt.plot( aa, D1a, color='green' )
# plt.plot( aa, D1z, ls='--', color='red' )
#
# plt.xlabel( 'a', fontsize=20 )
# plt.ylabel( 'D1(a)', fontsize=20 )
# def _plot_Dc_vs_z(self):
# fig = plt.figure()
# Dc = self._Dc[::8].copy()
# lopz = self._lopz[::8].copy()
# plt.plot( lopz, Dc, 'ro' )
# z = self.zDc( Dc )
# lopz = np.log10( 1 + z )
# plt.plot( lopz, Dc, 'g-' )
# plt.xlabel( r'$\log(1+z)$', fontsize=20 )
# plt.ylabel( r'$D_C \, [\rm Mpc/h]$', fontsize=20 )
def __str__(self):
print 'Cosmology instance'
print ' OmegaR: %5.3f' % self.OmegaR
print ' OmegaM: %5.3f' % self.OmegaM
print ' OmegaB: %5.3f' % self.OmegaB
print ' OmegaC: %5.3f' % self.OmegaC
print ' OmegaL: %5.3f' % self.OmegaL
print ' h: %5.3f' % self.h
print ' sigma_8: %5.3f' % self.sigma8
print ' n_s: %5.3f' % self.ns
print ' Y_P: %5.3f' % self.Yp
print
print ' H0: ', self.H0
print ' dH: ', self.dH0
print ' tH: ', self.tH0
print ' rho_crit_0: ', self.rho_crit0
print ' eps_crit_0: ', self.eps_crit0
return ''