Source Code for Module cosmolopy.density

  1  """Cosmological densities like matter density, baryon density, etc. 
  2  """ 
  3  import math 
  4   
  5  import numpy 
  6  import scipy 
  7  import scipy.special 
  8  import scipy.integrate as si 
  9   
 10  import constants as cc 
 11  import distance as cd 
 12  from distance import get_omega_k_0, set_omega_k_0 
 13   
 14 -def omega_M_z(z, **cosmo): 
 15      """Matter density omega_M as a function of redshift z. 
 16   
 17      Notes 
 18      ----- 
 19   
 20      From Lahav et al. (1991, MNRAS 251, 128) equations 11b-c. This is 
 21      equivalent to equation 10 of Eisenstein & Hu (1999 ApJ 511 5). 
 22   
 23      """ 
 24      if get_omega_k_0(**cosmo) == 0: 
 25          return 1.0 / (1. + (1. - cosmo['omega_M_0'])/ 
 26                        (cosmo['omega_M_0'] * (1. + z)**3.)) 
 27      else: 
 28          return (cosmo['omega_M_0'] * (1. + z)**3. /  
 29                  cd.e_z(z, **cosmo)**2.) 
 30   
 31 -def cosmo_densities(**cosmo): 
 32      """The critical and mean densities of the universe. 
 33   
 34      Returns 
 35      ------- 
 36      rho_crit and rho_0 in solar masses per cubic Megaparsec. 
 37   
 38      """ 
 39   
 40      omega_M_0 = cosmo['omega_M_0'] 
 41      h = cosmo['h'] 
 42   
 43      rho_crit = 3. * (h * cc.H100_s)**2. / (8. * math.pi * cc.G_const_Mpc_Msun_s) 
 44      rho_0 = omega_M_0 * rho_crit 
 45       
 46      #print " Critical density rho_crit = %.3g Msun/Mpc^3" % rho_crit 
 47      #print " Matter density      rho_0 = %.3g Msun/Mpc^3" % rho_0 
 48   
 49      return rho_crit, rho_0 
 50   
 51 -def get_X_Y(**cosmo): 
 52      """The fraction of baryonic mass in hydrogen and helium. 
 53   
 54      Assumes X_H + Y_He = 1. 
 55   
 56      You must specify either 'X_H', or 'Y_He', or both. 
 57      """ 
 58      if 'X_H' in cosmo and 'Y_He' not in cosmo: 
 59          X_H = cosmo['X_H'] 
 60          Y_He = 1. - X_H 
 61      elif 'Y_He' in cosmo and 'X_H' not in cosmo: 
 62          Y_He = cosmo['Y_He'] 
 63          X_H = 1. - Y_He 
 64      else: 
 65          X_H = cosmo['X_H'] 
 66          Y_He = cosmo['Y_He'] 
 67      return X_H, Y_He 
 68       
 69 -def baryon_densities(**cosmo): 
 70      """Hydrogen number density at z=0. 
 71   
 72      Parameters 
 73      ---------- 
 74   
 75      cosmo: cosmological parameters 
 76   
 77      parameters used: 'omega_b_0', 'X_H' and/or 'Y_He', plus those 
 78      needed by cosmo_densities. 
 79          
 80   
 81      Returns 
 82      ------- 
 83   
 84      rho_crit, rho_0, n_He_0, n_H_0 
 85   
 86      The first two are in units of solar masses per cubic 
 87      Megaparsec. The later two are in number per cubic Megaparsec. 
 88       
 89      """ 
 90       
 91      X_H, Y_He = get_X_Y(**cosmo) 
 92   
 93      rho_crit, rho_0 = cosmo_densities(**cosmo) 
 94   
 95      n_H_0 = (rho_crit * cosmo['omega_b_0'] * X_H * cc.M_sun_g /  
 96               cc.m_H_g) 
 97      n_He_0 = (rho_crit * cosmo['omega_b_0'] * Y_He * cc.M_sun_g /  
 98                cc.m_He_g) 
 99      #    print " Hydrogen number density n_H_0 = %.4g (Mpc^-3)" % n_H_0 
100      return rho_crit, rho_0, n_He_0, n_H_0 
101