#!/Users/weinberg.21/anaconda3/bin/python """ cosmodist2.py -- compute cosmology quantities, w0-wa model cosmodist2.py Omega_{m,0} Omega_{k,0} H0 w0 wa [zoption] Omega_{m,0} = present day value of Omega_m Omega_{k,0} = present day value of Omega_k, set to zero for flat universe H0 = present day Hubble constant in km/s/Mpc [set to zero for theta_* norm] w0, wa specify DE equation of state as w(a) = w0 + wa*(1-a) [zoption] = "smooth" "desi" or "des" (default is "smooth") for which redshifts are output --- For a cosmological constant, w0=-1, wa=0 --- Omega_{DE,0} will be calculated as 1 - Omega_{m,0}-Omega_{k,0}-Omega{r,0} Returns z Omega_m(z) t_loobkack(z) (Gyr) d_c(z) = \int_0^z H_0/H(z') dz' (dimensionless) D_C(z) = (c/H_0) d_c(z) (Mpc, for specified H0) D_M(z) = comoving ang. diam. distance (Mpc) D_H(z) = c/H(z) (Mpc) D_M/rd (dimensionless) D_H/rd (dimensionless) D_V/rd (dimensionless) Gfac (dimensionless) The angle-averaged D_V = (z D_H D_M^2)^{1/3} (eq. 2.6 of 2403.03002) rd can be fixed or scaled from other parameters, as seen below. Gfac uses eq. 16 of Weinberg++2013 as an approxmation for G(z)/G(z=0) If H0 is <=0, the value of H0 will be chosen to give D_M(z=1090)/rd = 94.23, which is the value found for flat LCDM with Omega_m=0.3145 and H0=67.36. For reference on parameter values, note that the cosmoprimo package implements Planck 2018 as omega_b = 0.02237 omega_cdm = 0.1200 h=0.6736 implying Omega_bc = (omega_cdm+omega_b)/h**2 = 0.3138 and Omega_m = 0.3145 when adding 0.06 eV neutrinos """ import numpy as np # may be used in integrand import sys # needed for command line reads below from scipy.integrate import quad # quadrature integration # functionality of a printf statement to stdout def printf(format, *args): sys.stdout.write(format %args) if (len(sys.argv)<6): print("\n") print("cosmodist2.py Omega_{m,0) Omega_{k,0} H0 w0 wa [zoption]\n") sys.exit() # read parameters omega_m=float(sys.argv[1]) omega_k=float(sys.argv[2]) h00=float(sys.argv[3]) # if <=0, h0 will be set by D_M(z=1090)/rd w0=float(sys.argv[4]) wa=float(sys.argv[5]) if (len(sys.argv)==6): zoption="smooth" else: zoption=sys.argv[6] if (zoption=="desi"): zout=np.array([0.295,0.510,0.706,0.934,1.321,1.484,2.330]) elif (zoption=="des"): zout=np.array([0.049,0.199,0.284,0.365,0.531,0.693,0.914]) else: zout=np.array([0.02,0.05,0.1,0.2,0.3,0.4,0.5,0.7,0.9,1.2,1.5,1.8,2.1,2.4,2.7]) # compute radiation density parameter omega_r, using values from Huterer tab 3.1 if (h00<=0): h0=67 # initial value, will be refined else: h0=h00 hh=h0/100. omega_gamma=2.47e-5/hh**2 omega_nu=1.68e-5/hh**2 # note that this is for massless neutrinos omega_r=omega_gamma+omega_nu # compute omega_de omega_de=1.-(omega_m+omega_k+omega_r) # compute sound horizon distance rd using Eq. 2 of DESI DR2 paper obh2=0.02236 # scaled value based on Planck neff=3.04 # in computing omega_bc, subtract the small contribution from 0.06 eV neutrinos obch2=omega_m*(hh**2)-0.06/93.14 rd=147.05*((obch2/0.1432)**-0.23)*((neff/3.04)**-0.1)*((obh2/0.02236)**-0.13) # or fix to fiducial value from CMB constraints (quoted error is +/- 0.26) # rd=147.05 # ratio of DE density at redshift z to present day density def de_ratio(z): if (wa==0.0): deratio=(1+z)**(3.*(1+w0)) else: az=1./(1+z) deratio=az**(-3.*(1+w0+wa))*np.exp(-3.*wa*(1-az)) return(deratio) # comoving distance integrand def dc_integrand(z): return(1./np.sqrt(omega_m*(1+z)**3+ omega_k*(1+z)**2+ omega_de*de_ratio(z)+ omega_r*(1+z)**4)) # lookback time integrand def tlookback_integrand(z): return(1./((1+z)*np.sqrt(omega_m*(1+z)**3+ omega_k*(1+z)**2+ omega_de*de_ratio(z)+ omega_r*(1+z)**4))) # time integrand forward from t=0, in a instead of z def tforward_integrand(a): znow=1./a-1 return(1./(a*np.sqrt(omega_m/(a**3)+ omega_k/(a**2)+ omega_de*de_ratio(znow)+ omega_r/(a**4) ))) # growth factor integrand (see Weinberg++2013, Eqs. 16 & 17) def growth_integrand(z): omega_mz=omega_m*(1+z)**3/(omega_m*(1+z)**3+ omega_k*(1+z)**2+ omega_de*de_ratio(z)+ omega_r*(1+z)**4 ) gamma=0.55+0.05*(1+w0+0.5*wa) return(((omega_mz)**gamma)/(1+z)) # dimensionless comoving angular diameter distance (without dh0 factor) def dmfunc(zz,ok): dc,err=quad(dc_integrand,0.0,zz) # quad returns integral and error if (ok==0): dmz=dc elif (ok>0): dmz=np.sinh(np.sqrt(ok)*dc)/np.sqrt(ok) else: dmz=np.sin(np.sqrt(-ok)*dc)/np.sqrt(-ok) return(dmz) dh0=2.997925e5/h0 # c/H0 in Mpc, if H0 is in km/s/Mpc th=9.7778*100./h0 # 1/H0 in Gyr, if H0 is in km/s/Mpc # compute ratio of D_M(z=1090)/rd, connected to angular scale of CMB peaks z=1090 ratio_zrec = dmfunc(z,omega_k)*dh0/rd ############################################################################## # THIS SECTION IS SKIPPED IF H0 > 0 IS SPECIFIED AS INPUT # if h00 <= 0, iterate a fixed number of times to match CMB-LCDM value of # D_M(z=1090)/rd if (h00<=0): # print(h0,rd,ratio_zrec) target=94.23 # D_M(z=1090)/rd for Omega_m=0.3138, H0=67.36 tolerance=1.e-4 # will try to match to this fractional tolerance nitermax=4 # maximum allowed iterations niter=0 while ((np.abs(ratio_zrec/target-1.)>tolerance) and (niter