From 454e8842885272f8b8a13a4530c69a792375137c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ga=C3=A9tan=20Facchinetti?= Date: Wed, 31 Jan 2024 12:09:19 +0100 Subject: [PATCH] removing useless packages --- .github/workflows/CI.yml | 2 +- Project.toml | 6 - README.md | 1 + src/Background.jl | 286 +++++++++++++++++++++++++++++++++++++++ src/BackgroundCosmo.jl | 271 ------------------------------------- src/CosmoTools.jl | 6 +- src/MassFunction.jl | 1 + src/PowerSpectrum.jl | 19 ++- src/TransferFunction.jl | 39 +++--- 9 files changed, 324 insertions(+), 307 deletions(-) create mode 100644 src/Background.jl delete mode 100644 src/BackgroundCosmo.jl diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index ccb6017..a9c62b6 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -23,7 +23,7 @@ jobs: fail-fast: false matrix: version: - - '1.8' + - '1.5' - '1.9' - 'nightly' os: diff --git a/Project.toml b/Project.toml index ee4848b..5cdf77f 100644 --- a/Project.toml +++ b/Project.toml @@ -5,18 +5,12 @@ version = "1.0.0-DEV" [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" -Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" -UnitfulAstro = "6112ee07-acf9-5e0f-b108-d242c714bf9f" [compat] -PhysicalConstants = "0.2" Roots = "2" -Unitful = "1" QuadGK = "2" -UnitfulAstro = "1" Documenter = "1" julia = "1.8" diff --git a/README.md b/README.md index 438362f..56f7645 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,4 @@ # CosmoTools [![Build Status](https://github.com/gaetanfacchinetti/CosmoTools.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/gaetanfacchinetti/CosmoTools.jl/actions/workflows/CI.yml?query=branch%3Amain) + diff --git a/src/Background.jl b/src/Background.jl new file mode 100644 index 0000000..ce862fb --- /dev/null +++ b/src/Background.jl @@ -0,0 +1,286 @@ +######################################################## +# Subfile of the module Cosmojuly.jl +# +# Contains functions related to the transfer functions +# +# author: Gaetan Facchinetti +# email: gaetanfacc@gmail.com +# +######################################################## + + +# First define some basic constants (the only one we will need) +const G_NEWTON::Float64 = 4.517103049894965e-48 # in Mpc^3 / Msun / s^(2) +const C_LIGHT::Float64 = 9.715611890180198e-15 # Mpc / s +const MPC_TO_KM::Float64 = 3.085677581491367e19 +const KM_TO_MPC::Float64 = 1.0/MPC_TO_KM + +export BkgCosmology, FLRW, FlatFLRW, planck18_bkg, edsPlanck18_bkg, convert_cosmo +export hubble_H0, hubble_E, hubble_E2 +export Ω, Ω_vs_a, Ω_m, Ω_r, Ω_Λ, mean_ρ, ρ_critical +export z_eq_mr, z_eq_Λm, z_to_a, a_to_z, temperature_CMB_K, k_eq_mr +export growth_factor, growth_factor_Carroll +export lookback_time, lookback_redshift, universe_age + +@doc raw""" + BkgCosmology + +Abstract type: generic background cosmology +""" +abstract type BkgCosmology end + +@doc raw""" + FLRW{T<:Real} <: BkgCosmology + +Defines a flat FLRW background cosmology +""" +struct FLRW{T<:Real} <: BkgCosmology + + # Hubble parameter + h::T + + # Abundances of the different components + Ω_χ0::T + Ω_b0::T + Ω_m0::T + Ω_r0::T + Ω_γ0::T + Ω_ν0::T + Ω_Λ0::T + Ω_k0::T + + # Derived quantities + z_eq_mr::T + z_eq_Λm::T + k_eq_mr::T + + # Quantity with units + ρ_c0::T + T0_CMB_K::T + +end + +Base.length(::BkgCosmology) = 1 +Base.iterate(iter::BkgCosmology) = (iter, nothing) +Base.iterate(::BkgCosmology, state::Nothing) = nothing + +""" Convert bkg_cosmo object attributes to another type """ +convert_cosmo(::Type{T}, bkg_cosmo::FLRW{<:Real} = planck18_bkg) where {T<:Real} = FLRW{T}([convert(T, getfield(bkg_cosmo, field)) for field in fieldnames(typeof(bkg_cosmo))]...) + +# First definition of the abundances +hubble_E2(z::Real, Ω_m0::Real, Ω_r0::Real, Ω_Λ0::Real, Ω_k0::Real) = Ω_m0 * (1+z)^3 + Ω_r0 * (1+z)^4 + Ω_Λ0 + Ω_k0 * (1+z)^2 + +Ω_m(z::Real, Ω_m0::Real, Ω_r0::Real, Ω_Λ0::Real, Ω_k0::Real) = Ω_m0 * (1+z)^3 / hubble_E2(z, Ω_m0, Ω_r0, Ω_Λ0, Ω_k0) +Ω_r(z::Real, Ω_m0::Real, Ω_r0::Real, Ω_Λ0::Real, Ω_k0::Real) = Ω_r0 * (1+z)^4 / hubble_E2(z, Ω_m0, Ω_r0, Ω_Λ0, Ω_k0) +Ω_Λ(z::Real, Ω_m0::Real, Ω_r0::Real, Ω_Λ0::Real, Ω_k0::Real) = Ω_Λ0 / hubble_E2(z, Ω_m0, Ω_r0, Ω_Λ0, Ω_k0) + +@doc raw""" + FlatFLRW(h, Ω_χ0, Ω_b0; T0_CMB_K = 2.72548, Neff = 3.04) + +Creates a FlatFLRW instance + +# Arguments +- `h::Real` : hubble parameter (dimensionless) +- `Ω_χ0::Real`: cold dark matter abundance today (dimensionless) +- `Ω_b0::Real`: baryon abundance today (dimensionless) +- `T0_CMB_K::Real`: temperature of the CMB today (in Kelvin) +- `Neff::Real`: effective number of neutrinos (dimensionless) +""" +function FLRW(h::Real, Ω_χ0::Real, Ω_b0::Real, Ω_k0::Real=0; T0_CMB_K::Real = 2.72548, Neff::Real = 3.04, EdS::Bool = false) + + # Derived abundances + Ω_γ0 = EdS ? 0.0 : 4.48131e-7 * T0_CMB_K^4 / h^2 + Ω_ν0 = EdS ? 0.0 : Neff * Ω_γ0 * (7 / 8) * (4 / 11)^(4 / 3) + Ω_r0 = Ω_γ0 + Ω_ν0 + Ω_m0 = Ω_χ0 + Ω_b0 + Ω_Λ0 = 1.0 - Ω_m0 - Ω_r0 + + ρ_c0 = 3/(8*π*G_NEWTON) * (100 * h * KM_TO_MPC)^2 + + z_eq_mr = 0.0 + z_eq_Λm = 0.0 + + try + z_eq_mr = EdS ? NaN : exp(Roots.find_zero( y -> Ω_r(exp(y), Ω_m0, Ω_r0, Ω_Λ0, Ω_k0) - Ω_m(exp(y), Ω_m0, Ω_r0, Ω_Λ0, Ω_k0), (-10, 10), Roots.Bisection())) + z_eq_Λm = exp(Roots.find_zero( y -> Ω_Λ(exp(y), Ω_m0, Ω_r0, Ω_Λ0, Ω_k0) - Ω_m(exp(y), Ω_m0, Ω_r0, Ω_Λ0, Ω_k0), (-10, 10), Roots.Bisection())) + catch e + println("Impossible to definez z_eq_mr and/or z_eq_Λm for this cosmology") + println("Error: ", e) + end + + k_eq_mr = EdS ? NaN : Ω_r0 / Ω_m0 * (100. * h * sqrt(Ω_m0 * (1+z_eq_mr)^3 + Ω_r0 * (1+z_eq_mr)^4 + Ω_Λ0 + Ω_k0 * (1+z_eq_mr)^2) * KM_TO_MPC / C_LIGHT) + + return FLRW(promote(h, Ω_χ0, Ω_b0, Ω_m0, Ω_r0, Ω_γ0, Ω_ν0, Ω_Λ0, Ω_k0, z_eq_mr, z_eq_Λm, k_eq_mr, ρ_c0, T0_CMB_K)...) +end + + +######################################################### +# Predefinition of some cosmologies +""" + planck18_bkg::FlatFLRW = FlatFLRW(0.6736, 0.26447, 0.04930) + +Flat FLRW cosmology with [*Planck18*](https://arxiv.org/abs/1807.06209) parameters +""" +const planck18_bkg::FLRW = FLRW(0.6736, 0.26447, 0.04930) + +""" + edsplanck18_bkg::FlatFLRW = FlatFLRW(0.6736, 0.3, 0) + +Flat FLRW cosmology with [*Planck18*](https://arxiv.org/abs/1807.06209) parameters and no baryons +""" +const edsPlanck18_bkg::FLRW = FLRW(0.6736, 0.3, 0, EdS = true) +######################################################### + + +######################################################### +# Definition of the hubble parameter and some other related quantities + +""" CMB temperature (in K) of the Universe at redshift `z` (by default z=0) for the cosmology `bkg_cosmo` """ +temperature_CMB_K(z::Real = 0.0, bkg_cosmo::BkgCosmology = planck18_bkg) = bkg_cosmo.T0_CMB_K * (1+z) + +""" Hubble constant H0 (in km/s/Mpc) for the cosmology `bkg_cosmo` """ +hubble_H0(bkg_cosmo::BkgCosmology = planck18_bkg) = 100 * bkg_cosmo.h + +""" Hubble constant H0 (in 1/s) for the cosmology `bkg_cosmo` """ +hubble_H0_s(bkg_cosmo::BkgCosmology = planck18_bkg) = hubble_H0(bkg_cosmo) * KM_TO_MPC + +""" Hubble evolution (no dimension) of the Universe squared at redshift `z` (by default z=0) for the cosmology `bkg_cosmo` """ +hubble_E2(z::Real = 0.0, bkg_cosmo::BkgCosmology = planck18_bkg) = hubble_E2(z, bkg_cosmo.Ω_m0, bkg_cosmo.Ω_r0, bkg_cosmo.Ω_Λ0, bkg_cosmo.Ω_k0) + +""" Hubble evolution (no dimension) of the Universe at redshift `z` (by default z=0) for the cosmology `bkg_cosmo` """ +hubble_E(z::Real = 0.0, bkg_cosmo::BkgCosmology = planck18_bkg) = sqrt(hubble_E2(z, bkg_cosmo)) + +""" Hubble rate H(z) (in km/s/Mpc) of the Universe at redshift `z` (by default z=0) for the cosmology `bkg_cosmo` """ +hubble_H(z::Real = 0.0, bkg_cosmo::BkgCosmology = planck18_bkg) = hubble_E(z, bkg_cosmo) .* hubble_H0(bkg_cosmo) + +""" k at equivalence between matter and radiation """ +k_eq_mr(bkg_cosmo::BkgCosmology) = bkg_cosmo.k_eq_mr +######################################################### + + +######################################################### +# Definition of the densities and abundances +# All densities are in units of Msun / Mpc^3 + +export Species, Radiation, Photons, Neutrinos, Matter, ColdDarkMatter, Baryons, DarkEnergy, Curvature + +abstract type Species end +abstract type Radiation <: Species end +abstract type Photons <: Radiation end +abstract type Neutrinos <: Radiation end +abstract type Matter <: Species end +abstract type ColdDarkMatter <: Matter end +abstract type Baryons <: Matter end +abstract type DarkEnergy <: Species end +abstract type Curvature <: Species end + +Ω(z::Real, n::Int, Ω0::Real, bkg_cosmo::BkgCosmology) = Ω0 * (1+z)^n / hubble_E2(z, bkg_cosmo) + +Ω0(::Type{Matter}, bkg_cosmo::BkgCosmology) = bkg_cosmo.Ω_m0 +Ω0(::Type{Radiation}, bkg_cosmo::BkgCosmology) = bkg_cosmo.Ω_r0 +Ω0(::Type{DarkEnergy}, bkg_cosmo::BkgCosmology) = bkg_cosmo.Ω_Λ0 +Ω0(::Type{Photons}, bkg_cosmo::BkgCosmology) = bkg_cosmo.Ω_γ0 +Ω0(::Type{Neutrinos}, bkg_cosmo::BkgCosmology) = bkg_cosmo.Ω_ν0 +Ω0(::Type{ColdDarkMatter}, bkg_cosmo::BkgCosmology) = bkg_cosmo.Ω_χ0 +Ω0(::Type{Baryons}, bkg_cosmo::BkgCosmology) = bkg_cosmo.Ω_b0 +Ω0(::Type{Curvature}, bkg_cosmo::BkgCosmology) = bkg_cosmo.Ω_k0 + +index_ρ(::Type{T}) where {T<:Matter} = 3 +index_ρ(::Type{T}) where {T<:Radiation} = 4 +index_ρ(::Type{DarkEnergy})::Int = 0 + +@doc raw""" + Ω(T, z = 0, bkg_cosmo = planck18_bkg) where {T<:Species} + +Abundance of species `T` at redshift `z` and for the background cosmology `bkg_cosmo` +""" +Ω(::Type{T}, z::Real = 0, bkg_cosmo::BkgCosmology = planck18_bkg) where {T<:Species} = Ω(z, index_ρ(T), Ω0(T, bkg_cosmo), bkg_cosmo) +Ω_vs_a(::Type{T}, a::Real = 1, bkg_cosmo::BkgCosmology = planck18_bkg) where {T<:Species} = Ω(T, a_to_z(a), bkg_cosmo) + + +""" Critical density (in Msun/Mpc^3) of the Universe at redshift `z` (by default z=0) for the cosmology `bkg_cosmo` """ +ρ_critical(z::Real = 0, bkg_cosmo::BkgCosmology = planck18_bkg) = bkg_cosmo.ρ_c0 * hubble_E2(z, bkg_cosmo) + +mean_ρ(::Type{Radiation}, z::Real, bkg_cosmo::BkgCosmology) = bkg_cosmo.Ω_r0 * bkg_cosmo.ρ_c0 * (1+z)^4 +mean_ρ(::Type{Photons}, z::Real, bkg_cosmo::BkgCosmology) = bkg_cosmo.Ω_γ0 * bkg_cosmo.ρ_c0 * (1+z)^4 +mean_ρ(::Type{Neutrinos}, z::Real, bkg_cosmo::BkgCosmology) = bkg_cosmo.Ω_ν0 * bkg_cosmo.ρ_c0 * (1+z)^4 +mean_ρ(::Type{ColdDarkMatter}, z::Real, bkg_cosmo::BkgCosmology) = bkg_cosmo.Ω_χ0 * bkg_cosmo.ρ_c0 * (1+z)^3 +mean_ρ(::Type{Matter}, z::Real, bkg_cosmo::BkgCosmology) = bkg_cosmo.Ω_m0 * bkg_cosmo.ρ_c0 * (1+z)^3 +mean_ρ(::Type{Baryons}, z::Real, bkg_cosmo::BkgCosmology) = bkg_cosmo.Ω_b0 * bkg_cosmo.ρ_c0 * (1+z)^3 +mean_ρ(::Type{Curvature}, z::Real, bkg_cosmo::BkgCosmology) = bkg_cosmo.Ω_k0 * bkg_cosmo.ρ_c0 * (1+z)^2 +mean_ρ(::Type{DarkEnergy}, z::Real, bkg_cosmo::BkgCosmology) = bkg_cosmo.Ω_Λ0 * bkg_cosmo.ρ_c0 + +@doc raw""" + mean_ρ(T, z = 0, bkg_cosmo = planck18_bkg) where {T<:Species} + +Energy density of species `T` at redshift `z` and for the background cosmology `bkg_cosmo` +""" +mean_ρ(::Type{T}, z::Real = 0, bkg_cosmo::BkgCosmology = planck18_bkg) where {T<:Species} = mean_ρ(T, z, bkg_cosmo) + +######################################################### + + +######################################################### +# Definition of specific time quantities + +""" Convert redshift to scale parameter """ +z_to_a(z::Real) = 1 / (1 + z) +""" Convert scale parameter to redshift """ +a_to_z(a::Real) = 1 / a - 1 + +z_eq_mr(bkg_cosmo::BkgCosmology = planck18_bkg) = bkg_cosmo.z_eq_mr +z_eq_Λm(bkg_cosmo::BkgCosmology = planck18_bkg) = bkg_cosmo.z_eq_Λm +a_eq_mr(bkg_cosmo::BkgCosmology = planck18_bkg) = z_to_a(bkg_cosmo.z_eq_mr) +a_eq_Λm(bkg_cosmo::BkgCosmology = planck18_bkg) = z_to_a(bkg_cosmo.z_eq_Λm) + +δt_s(a0::Real, a1::Real, bkg_cosmo::BkgCosmology = planck18_bkg; kws...) = QuadGK.quadgk(a -> 1 / hubble_E(a_to_z(a), bkg_cosmo) / a, a0, a1, rtol=1e-3; kws...)[1] / (hubble_H0(bkg_cosmo) * km_to_Mpc) + +""" age of the universe (s)""" +universe_age(z::Float64 = 0.0, bkg_cosmo::BkgCosmology = planck18_bkg; kws...) = δt_s(0, z_to_a(z), bkg_cosmo; kws...) + +""" lookback time of the Universe (s) """ +lookback_time(z::Real, bkg_cosmo::BkgCosmology = planck18_bkg; kws...) = δt_s(z_to_a(z), 1, bkg_cosmo; kws...) + +""" lookback redshift of the Universe for t in (s) """ +lookback_redshift(t::Real, bkg_cosmo::BkgCosmology = planck18_bkg; kws...) = exp(Roots.find_zero(lnz -> lookback_time(exp(lnz), bkg_cosmo; kws...)-t, (-5, 10), Roots.Bisection())) +######################################################### + + +######################################################### +# Functions related to perturbations +""" + growth_factor(z, bkg_cosmo) + + Exact growth factor in a matter-Λ Universe computed from an integral + Carroll et al. 1992 (Mo and White p. 172) + Corresponds to D1(a=1) with the definition of Dodelson 2003 + +# Arguments +- z: redshift +- bkg_cosmo: background cosmology (default Planck18) +""" +function growth_factor(z::Real, bkg_cosmo::BkgCosmology = planck18_bkg; rtol=1e-6, kws...)::Real + norm = 2.5 * bkg_cosmo.Ω_m0 * sqrt(bkg_cosmo.Ω_m0 * (1+z)^3 + bkg_cosmo.Ω_Λ0) + return norm * QuadGK.quadgk(a -> (bkg_cosmo.Ω_m0 * a^(-1) + bkg_cosmo.Ω_Λ0 * a^2)^(-3/2), 0, z_to_a(z), rtol=rtol; kws...)[1] +end + +""" + growth_factor_Carroll(z, bkg_cosmo) + + Approximate growth factor in a matter-Λ Universe (faster than growth_factor) + Carroll et al. 1992 (Mo and White p. 172) + Corresponds to D1(a=1) with the definition of Dodelson 2003 + +# Arguments +- z: redshift +- bkg_cosmo: background cosmology (default Planck18) +""" +function growth_factor_Carroll(z::Real, bkg_cosmo=BkgCosmology = planck18_bkg)::Real + # Abundances in a Universe without radiation + _Ω_m = bkg_cosmo.Ω_m0 * (1+z)^3 / (bkg_cosmo.Ω_m0 * (1+z)^3 + bkg_cosmo.Ω_Λ0) + _Ω_Λ = bkg_cosmo.Ω_Λ0 / (bkg_cosmo.Ω_m0 * (1+z)^3 + bkg_cosmo.Ω_Λ0) + return 2.5*_Ω_m/(_Ω_m^(4.0/7.0) - _Ω_Λ + (1.0 + 0.5*_Ω_m) * (1.0 + 1.0/70.0*_Ω_Λ))/(1+z) +end +######################################################### + + diff --git a/src/BackgroundCosmo.jl b/src/BackgroundCosmo.jl deleted file mode 100644 index 5cfd15c..0000000 --- a/src/BackgroundCosmo.jl +++ /dev/null @@ -1,271 +0,0 @@ -######################################################## -# Subfile of the module Cosmojuly.jl -# -# Contains functions related to the transfer functions -# -# author: Gaetan Facchinetti -# email: gaetanfacc@gmail.com -# -######################################################## - -export BkgCosmology, FLRW, FlatFLRW, planck18_bkg, edsPlanck18_bkg, convert_cosmo -export hubble_H0, hubble_E, hubble_E2 -export Ω, Ω_vs_a, Ω_m, Ω_r, Ω_Λ, mean_ρ, ρ_critical -export z_eq_mr, z_eq_Λm, z_to_a, a_to_z, temperature_CMB_K, k_eq_mr_Mpc -export growth_factor, growth_factor_Carroll -export lookback_time, lookback_redshift, universe_age - -@doc raw""" - BkgCosmology{T<:Real} - -Abstract type: generic background cosmology -""" -abstract type BkgCosmology{T<:Real} end - -@doc raw""" - FLRW{T<:Real} <: BkgCosmology{T} - -Abstract type: Generic background cosmology of type FLRW -""" -abstract type FLRW{T<:Real} <: BkgCosmology{T} end - -@doc raw""" - FlatFLRW{T<:Real} <: FLRW{T} - -Defines a flat FLRW cosmology -""" -struct FlatFLRW{T<:Real} <: FLRW{T} - - # Hubble parameter - h::T - - # Abundances of the different components - Ω_χ0::T - Ω_b0::T - Ω_m0::T - Ω_r0::T - Ω_γ0::T - Ω_ν0::T - Ω_Λ0::T - - # Derived quantities - z_eq_mr::T - z_eq_Λm::T - k_eq_mr_Mpc::T - - # Quantity with units - ρ_c0::T - T0_CMB_K::T - -end - -Base.length(::BkgCosmology) = 1 -Base.iterate(iter::BkgCosmology) = (iter, nothing) -Base.iterate(::BkgCosmology, state::Nothing) = nothing - - -""" Convert cosmo object attributes to another type """ -convert_cosmo(::Type{T}, cosmo::FlatFLRW = planck18_bkg) where {T<:Real} = FlatFLRW{T}([convert(T, getfield(cosmo, field)) for field in fieldnames(typeof(cosmo))]...) - - -# First definition of the abundances -hubble_E2(z::Real, Ω_m0::Real, Ω_r0::Real, Ω_Λ0::Real) = Ω_m0 * (1+z)^3 + Ω_r0 * (1+z)^4 + Ω_Λ0 - -Ω_m(z::Real, Ω_m0::Real, Ω_r0::Real, Ω_Λ0::Real) = Ω_m0 * (1+z)^3 / hubble_E2(z, Ω_m0, Ω_r0, Ω_Λ0) -Ω_r(z::Real, Ω_m0::Real, Ω_r0::Real, Ω_Λ0::Real) = Ω_r0 * (1+z)^4 / hubble_E2(z, Ω_m0, Ω_r0, Ω_Λ0) -Ω_Λ(z::Real, Ω_m0::Real, Ω_r0::Real, Ω_Λ0::Real) = Ω_Λ0 / hubble_E2(z, Ω_m0, Ω_r0, Ω_Λ0) - -@doc raw""" - FlatFLRW(h, Ω_χ0, Ω_b0; T0_CMB_K = 2.72548, Neff = 3.04) - -Creates a FlatFLRW instance - -# Arguments -- `h::Real` : hubble parameter (dimensionless) -- `Ω_χ0::Real`: cold dark matter abundance today (dimensionless) -- `Ω_b0::Real`: baryon abundance today (dimensionless) -- `T0_CMB_K::Real`: temperature of the CMB today (in Kelvin) -- `Neff::Real`: effective number of neutrinos (dimensionless) -""" -function FlatFLRW(h::Real, Ω_χ0::Real, Ω_b0::Real; T0_CMB_K::Real = 2.72548, Neff::Real = 3.04, EdS::Bool = false) - - # Derived abundances - Ω_γ0 = EdS ? 0.0 : 4.48131e-7 * T0_CMB_K^4 / h^2 - Ω_ν0 = EdS ? 0.0 : Neff * Ω_γ0 * (7 / 8) * (4 / 11)^(4 / 3) - Ω_r0 = EdS ? 0.0 : Ω_γ0 + Ω_ν0 - Ω_m0 = Ω_χ0 + Ω_b0 - Ω_Λ0 = 1 - Ω_m0 - Ω_r0 - - ρ_c0 = 3/(8*π*G_NEWTON) * (100 * h * Unitful.km / Unitful.s / UnitfulAstro.Mpc )^2 / (UnitfulAstro.Msun / UnitfulAstro.Mpc^3) - - z_eq_mr = 0.0 - z_eq_Λm = 0.0 - - try - z_eq_mr = EdS ? NaN : exp(Roots.find_zero( y -> Ω_r(exp(y), Ω_m0, Ω_r0, Ω_Λ0) - Ω_m(exp(y), Ω_m0, Ω_r0, Ω_Λ0), (-10, 10), Roots.Bisection())) - z_eq_Λm = exp(Roots.find_zero( y -> Ω_Λ(exp(y), Ω_m0, Ω_r0, Ω_Λ0) - Ω_m(exp(y), Ω_m0, Ω_r0, Ω_Λ0), (-10, 10), Roots.Bisection())) - catch e - println("Impossible to definez z_eq_mr and/or z_eq_Λm for this cosmology") - println("Error: ", e) - end - - k_eq_mr_Mpc = EdS ? NaN : Ω_r0 / Ω_m0 * (100. * h * sqrt(Ω_m0 * (1+z_eq_mr)^3 + Ω_r0 * (1+z_eq_mr)^4 + Ω_Λ0) * Unitful.km / Unitful.s / C_LIGHT) |> Unitful.NoUnits - - return FlatFLRW(promote(h, Ω_χ0, Ω_b0, Ω_m0, Ω_r0, Ω_γ0, Ω_ν0, Ω_Λ0, z_eq_mr, z_eq_Λm, k_eq_mr_Mpc, ρ_c0, T0_CMB_K)...) -end - - -######################################################### -# Predefinition of some cosmologies -""" - planck18_bkg::FlatFLRW = FlatFLRW(0.6736, 0.26447, 0.04930) - -Flat FLRW cosmology with [*Planck18*](https://arxiv.org/abs/1807.06209) parameters -""" -const planck18_bkg::FlatFLRW = FlatFLRW(0.6736, 0.26447, 0.04930) - -""" - edsplanck18_bkg::FlatFLRW = FlatFLRW(0.6736, 0.3, 0) - -Flat FLRW cosmology with [*Planck18*](https://arxiv.org/abs/1807.06209) parameters and no baryons -""" -const edsPlanck18_bkg::FlatFLRW = FlatFLRW(0.6736, 0.3, 0, EdS = true) -######################################################### - -######################################################### -# Definition of the densities - -""" CMB temperature (in K) of the Universe at redshift `z` (by default z=0) for the cosmology `cosmo` """ -temperature_CMB_K(z::Real = 0.0, cosmo::BkgCosmology = planck18_bkg) = cosmo.T0_CMB_K * (1+z) - -""" Hubble constant H0 (in km/s/Mpc) for the cosmology `cosmo` """ -hubble_H0(cosmo::BkgCosmology = planck18_bkg) = 100 * cosmo.h - -""" Hubble evolution (no dimension) of the Universe squared at redshift `z` (by default z=0) for the cosmology `cosmo` """ -hubble_E2(z::Real = 0.0, cosmo::BkgCosmology = planck18_bkg) = hubble_E2(z, cosmo.Ω_m0, cosmo.Ω_r0, cosmo.Ω_Λ0) - -""" Hubble evolution (no dimension) of the Universe at redshift `z` (by default z=0) for the cosmology `cosmo` """ -hubble_E(z::Real = 0.0, cosmo::BkgCosmology = planck18_bkg) = sqrt(hubble_E2(z, cosmo)) - -""" Hubble rate H(z) (in km/s/Mpc) of the Universe at redshift `z` (by default z=0) for the cosmology `cosmo` """ -hubble_H(z::Real = 0.0, cosmo::BkgCosmology = planck18_bkg) = hubble_E(z, cosmo) .* hubble_H0(cosmo) -######################################################### - -######################################################### -# Definition of the densities and abundances -# All densities are in units of Msun / Mpc^3 - -export Species, Radiation, Photons, Neutrinos, Matter, ColdDarkMatter, Baryons, DarkEnergy, Curvature - -abstract type Species end -abstract type Radiation <: Species end -abstract type Photons <: Radiation end -abstract type Neutrinos <: Radiation end -abstract type Matter <: Species end -abstract type ColdDarkMatter <: Matter end -abstract type Baryons <: Matter end -abstract type DarkEnergy <: Species end -abstract type Curvature <: Species end - -Ω(z::Real, n::Int, Ω0::Real, cosmo::BkgCosmology) = Ω0 * (1+z)^n / hubble_E2(z, cosmo) - -Ω0(::Type{Matter}, cosmo::BkgCosmology) = cosmo.Ω_m0 -Ω0(::Type{Radiation}, cosmo::BkgCosmology) = cosmo.Ω_r0 -Ω0(::Type{DarkEnergy}, cosmo::BkgCosmology) = cosmo.Ω_Λ0 -Ω0(::Type{Photons}, cosmo::BkgCosmology) = cosmo.Ω_γ0 -Ω0(::Type{Neutrinos}, cosmo::BkgCosmology) = cosmo.Ω_ν0 -Ω0(::Type{ColdDarkMatter}, cosmo::BkgCosmology) = cosmo.Ω_χ0 -Ω0(::Type{Baryons}, cosmo::BkgCosmology) = cosmo.Ω_b0 - -index_ρ(::Type{T}) where {T<:Matter} = 3 -index_ρ(::Type{T}) where {T<:Radiation} = 4 -index_ρ(::Type{DarkEnergy})::Int = 0 - -@doc raw""" - Ω(T, z = 0, cosmo = planck18_bkg) where {T<:Species} - -Abundance of species `T` at redshift `z` and for the background cosmology `cosmo` -""" -Ω(::Type{T}, z::Real = 0.0, cosmo::BkgCosmology = planck18_bkg) where {T<:Species} = Ω(z, index_ρ(T), Ω0(T, cosmo), cosmo) -Ω_vs_a(::Type{T}, a::Real = 1.0, cosmo::BkgCosmology = planck18_bkg) where {T<:Species} = Ω(T, a_to_z(a), cosmo) - - -""" Critical density (in Msun/Mpc^3) of the Universe at redshift `z` (by default z=0) for the cosmology `cosmo` """ -ρ_critical(z::Real = 0.0, cosmo::BkgCosmology = planck18_bkg) = cosmo.ρ_c0 * hubble_E2(z, cosmo) - -mean_ρ(::Type{ColdDarkMatter}, z::Real, cosmo::BkgCosmology) = cosmo.Ω_χ0 * cosmo.ρ_c0 * (1+z)^3 -mean_ρ(::Type{Radiation}, z::Real, cosmo::BkgCosmology) = cosmo.Ω_r0 * cosmo.ρ_c0 * (1+z)^4 -mean_ρ(::Type{Photons}, z::Real, cosmo::BkgCosmology) = cosmo.Ω_γ0 * cosmo.ρ_c0 * (1+z)^4 -mean_ρ(::Type{Neutrinos}, z::Real, cosmo::BkgCosmology) = cosmo.Ω_ν0 * cosmo.ρ_c0 * (1+z)^4 -mean_ρ(::Type{Matter}, z::Real, cosmo::BkgCosmology) = cosmo.Ω_m0 * cosmo.ρ_c0 * (1+z)^3 -mean_ρ(::Type{Baryons}, z::Real, cosmo::BkgCosmology) = cosmo.Ω_b0 * cosmo.ρ_c0 * (1+z)^3 -mean_ρ(::Type{DarkEnergy}, z::Real, cosmo::BkgCosmology) = cosmo.Ω_Λ0 * cosmo.ρ_c0 - -@doc raw""" - mean_ρ(T, z = 0, cosmo = planck18_bkg) where {T<:Species} - -Energy density of species `T` at redshift `z` and for the background cosmology `cosmo` -""" -mean_ρ(::Type{T}, z::Real = 0.0, cosmo::BkgCosmology = planck18_bkg) where {T<:Species} = mean_ρ(T, z, cosmo) - -######################################################### - -k_eq_mr_Mpc(cosmo::BkgCosmology) = cosmo.k_eq_mr_Mpc - -######################################################### -# Definition of specific time quantities - -z_to_a(z::Real)::Real = 1 /(1 + z) -a_to_z(a::Real)::Real = 1 / a - 1 -z_eq_mr(cosmo::BkgCosmology = planck18_bkg) = cosmo.z_eq_mr -z_eq_Λm(cosmo::BkgCosmology = planck18_bkg) = cosmo.z_eq_Λm -a_eq_mr(cosmo::BkgCosmology = planck18_bkg) = z_to_a(cosmo.z_eq_mr) -a_eq_Λm(cosmo::BkgCosmology = planck18_bkg) = z_to_a(cosmo.z_eq_Λm) -δt_s(a0::Real, a1::Real, cosmo::BkgCosmology = planck18_bkg; kws...) = QuadGK.quadgk(a -> 1.0 / hubble_E(a_to_z(a), cosmo) / a, a0, a1, rtol=1e-3; kws...)[1] / (hubble_H0(cosmo) * Unitful.km / Mpc) |> Uniful.NoUnits -""" age of the universe (s)""" -universe_age(z::Real=0, cosmo::BkgCosmology = planck18_bkg; kws...) = δt_s(0, z_to_a(z), cosmo; kws...) -""" lookback time of the Universe (s) """ -lookback_time(z::Real, cosmo::BkgCosmology = planck18_bkg; kws...) = δt_s(z_to_a(z), 1, cosmo; kws...) -""" lookback redshift of the Universe for t in (s) """ -lookback_redshift(t::Real, cosmo::BkgCosmology = planck18_bkg; kws...) = exp(Roots.find_zero(lnz -> lookback_time(exp(lnz), cosmo; kws...)-t, (-5, 10), Roots.Bisection())) -######################################################### - - -######################################################### -# Functions related to perturbations -""" - growth_factor(z, cosmo) - - Exact growth factor in a matter-Λ Universe computed from an integral - Carroll et al. 1992 (Mo and White p. 172) - Corresponds to D1(a=1) with the definition of Dodelson 2003 - -# Arguments -- z: redshift -- cosmo: background cosmology (default Planck18) -""" -function growth_factor(z::Real, cosmo::BkgCosmology = planck18_bkg; rtol=1e-6, kws...)::Real - norm = 2.5 * cosmo.Ω_m0 * sqrt(cosmo.Ω_m0 * (1+z)^3 + cosmo.Ω_Λ0) - return norm * QuadGK.quadgk(a -> (cosmo.Ω_m0 * a^(-1) + cosmo.Ω_Λ0 * a^2)^(-3/2), 0, z_to_a(z), rtol=rtol; kws...)[1] -end - -""" - growth_factor_Carroll(z, cosmo) - - Approximate growth factor in a matter-Λ Universe (faster than growth_factor) - Carroll et al. 1992 (Mo and White p. 172) - Corresponds to D1(a=1) with the definition of Dodelson 2003 - -# Arguments -- z: redshift -- cosmo: background cosmology (default Planck18) -""" -function growth_factor_Carroll(z::Real, cosmo=BkgCosmology = planck18_bkg)::Real - # Abundances in a Universe with no radiation - _Ω_m = cosmo.Ω_m0 * (1+z)^3 / (cosmo.Ω_m0 * (1+z)^3 + cosmo.Ω_Λ0) - _Ω_Λ = cosmo.Ω_Λ0 / (cosmo.Ω_m0 * (1+z)^3 + cosmo.Ω_Λ0) - return 2.5*_Ω_m/(_Ω_m^(4.0/7.0) - _Ω_Λ + (1.0 + 0.5*_Ω_m) * (1.0 + 1.0/70.0*_Ω_Λ))/(1+z) -end -######################################################### - - diff --git a/src/CosmoTools.jl b/src/CosmoTools.jl index c924fb4..61f73f5 100644 --- a/src/CosmoTools.jl +++ b/src/CosmoTools.jl @@ -1,10 +1,6 @@ module CosmoTools -import QuadGK, Roots, Unitful, UnitfulAstro -import PhysicalConstants.CODATA2018: G as G_NEWTON -import PhysicalConstants.CODATA2018: c_0 as C_LIGHT - -export Ω +import QuadGK, Roots include("MassFunction.jl") diff --git a/src/MassFunction.jl b/src/MassFunction.jl index 5432e62..830b4e8 100644 --- a/src/MassFunction.jl +++ b/src/MassFunction.jl @@ -15,6 +15,7 @@ abstract type SethTormen <: MassFunctionType end f_mass_function(ν::Real, ::Type{PressSchechter}) = sqrt(2.0 / π) * exp(-ν^2 / 2.0) function f_mass_function(ν::Real, ::Type{SethTormen}) + _a = 0.707; _νp = sqrt(_a) * ν; _q = 0.3; diff --git a/src/PowerSpectrum.jl b/src/PowerSpectrum.jl index e7bb53c..3d72087 100644 --- a/src/PowerSpectrum.jl +++ b/src/PowerSpectrum.jl @@ -9,14 +9,21 @@ export σ2_mps, dσ2_mps_dR, σ_mps, dσ_mps_dR, σ2_mps_M, dσ2_mps_dM, σ_mps_ # COSMOLOGY STRUCTURE ##################### + +@doc raw""" + Cosmology + +Defines a generic cosmology +""" struct Cosmology name::String - bkg::BkgCosmology{<:Real} + bkg::BkgCosmology power_spectrum::Function transfer_function_model::TransferFunctionModel end -function Cosmology(name::String, bkg::BkgCosmology{<:Real}, power_spectrum::Function, ::Type{T} = EH98) where {T<:TransferFunctionModel} + +function Cosmology(name::String, bkg::BkgCosmology, power_spectrum::Function, ::Type{T} = EH98) where {T<:TransferFunctionModel} return Cosmology(name, bkg, power_spectrum, T(bkg)) end @@ -47,7 +54,7 @@ function matter_power_spectrum( dimensionless = false, with_baryons::Bool = true) - _c_over_H0_Mpc = C_LIGHT / (hubble_H0(cosmology.bkg) * Unitful.km / Unitful.s / UnitfulAstro.Mpc) / UnitfulAstro.Mpc |> Unitful.NoUnits + _c_over_H0_Mpc = C_LIGHT / hubble_H0_s(cosmology.bkg) _D1_z = growth_function(z, cosmology.bkg) _tf = transfer_function(k, cosmology.transfer_function_model, with_baryons = with_baryons) _prefactor = !dimensionless ? 1 : k^3/(2*π^2) @@ -56,9 +63,9 @@ function matter_power_spectrum( end -############################################### -# QUANTITIES THAT DEPEND ON THE WINDOW FUNCTION -############################################### +################################################### +# QUANTITIES DERIVED FROM THE MATTER POWER SPECTRUM +################################################### abstract type Window end abstract type TopHatWindow <: Window end diff --git a/src/TransferFunction.jl b/src/TransferFunction.jl index e1a817c..f5ad7ff 100644 --- a/src/TransferFunction.jl +++ b/src/TransferFunction.jl @@ -8,7 +8,7 @@ # ######################################################## -include("BackgroundCosmo.jl") +include("Background.jl") export transfer_function, TrivialTF, EH98, TransferFunctionModel, EH98_planck18 @@ -36,7 +36,7 @@ end g_func(y::Real) = (-6.0*sqrt(1.0+y) + (2.0+3.0*y)*log((sqrt(1.0+y)+1)/(sqrt(1.0+y)-1.0))) * y -function EH98(Ω_m0_h2::Real, Ω_b0_h2::Real, Ω_χ0_h2::Real, z_eq_mr::Real, k_eq_mr_Mpc::Real, ::Type{T}; T0_CMB_K::Real = 2.72548) where {T<:Real} +function EH98(Ω_m0_h2::Real, Ω_b0_h2::Real, Ω_χ0_h2::Real, z_eq_mr::Real, k_eq_mr::Real, ::Type{T}; T0_CMB_K::Real = 2.72548) where {T<:Real} Θ27::T = T0_CMB_K / 2.7 @@ -48,7 +48,7 @@ function EH98(Ω_m0_h2::Real, Ω_b0_h2::Real, Ω_χ0_h2::Real, z_eq_mr::Real, k_ # sound horizon R_drag = 31.5 * Ω_b0_h2 * Θ27^(-4) * 1e+3 / z_drag R_eq = 31.5 * Ω_b0_h2 * Θ27^(-4) * 1e+3 / z_eq_mr - sound_horizon_Mpc::T = 2. / (3. * k_eq_mr_Mpc) * sqrt(6. / R_eq) * log((sqrt(1. + R_drag) + sqrt(R_drag + R_eq))/(1+sqrt(R_eq))) + sound_horizon_Mpc::T = 2. / (3. * k_eq_mr) * sqrt(6. / R_eq) * log((sqrt(1. + R_drag) + sqrt(R_drag + R_eq))/(1+sqrt(R_eq))) # α_c a1 = (46.9 * Ω_m0_h2)^0.670 * (1.0 + (32.1 * Ω_m0_h2)^(-0.532) ) @@ -62,22 +62,23 @@ function EH98(Ω_m0_h2::Real, Ω_b0_h2::Real, Ω_χ0_h2::Real, z_eq_mr::Real, k_ # k_silk, α_b, and β_b k_Silk::T = 1.6 * ( Ω_b0_h2^0.52 ) * (Ω_m0_h2^0.73 ) * ( 1.0 + (10.4 * Ω_m0_h2)^(-0.95) ) - α_b::T = 2.07 * k_eq_mr_Mpc * sound_horizon_Mpc * (1.0 + R_drag)^(-3.0/4.0) * g_func( (1.0 + z_eq_mr) / (1.0 + z_drag) ) + α_b::T = 2.07 * k_eq_mr * sound_horizon_Mpc * (1.0 + R_drag)^(-3.0/4.0) * g_func( (1.0 + z_eq_mr) / (1.0 + z_drag) ) β_b::T = 0.5 + Ω_b0_h2 / Ω_m0_h2 + (3.0 - 2.0 * Ω_b0_h2/Ω_m0_h2) * sqrt( 1.0 + (17.2 * Ω_m0_h2)^2 ) EH98(convert(T, Ω_m0_h2), convert(T, Ω_b0_h2), convert(T, Ω_χ0_h2), convert(T, Θ27), z_drag, sound_horizon_Mpc, α_c, α_b, β_c, β_b, k_Silk) end -function EH98(cosmo::BkgCosmology{<:Real}, ::Type{T} = Float64) where {T<:Real} - Ω_m0_h2::T = cosmo.Ω_m0 * cosmo.h^2 - Ω_b0_h2::T = cosmo.Ω_b0 * cosmo.h^2 - Ω_χ0_h2::T = cosmo.Ω_χ0 * cosmo.h^2 +function EH98(bkg_cosmo::FLRW{<:Real}, ::Type{T} = Float64) where {T<:Real} - z_eq = z_eq_mr(cosmo) - k_eq_Mpc = k_eq_mr_Mpc(cosmo) + Ω_m0_h2::T = bkg_cosmo.Ω_m0 * bkg_cosmo.h^2 + Ω_b0_h2::T = bkg_cosmo.Ω_b0 * bkg_cosmo.h^2 + Ω_χ0_h2::T = bkg_cosmo.Ω_χ0 * bkg_cosmo.h^2 - return EH98(Ω_m0_h2, Ω_b0_h2, Ω_χ0_h2, z_eq, k_eq_Mpc, T, T0_CMB_K = cosmo.T0_CMB_K) + z_eq = z_eq_mr(bkg_cosmo) + k_eq_Mpc = k_eq_mr(bkg_cosmo) + + return EH98(Ω_m0_h2, Ω_b0_h2, Ω_χ0_h2, z_eq, k_eq_Mpc, T, T0_CMB_K = bkg_cosmo.T0_CMB_K) end const EH98_planck18 = EH98(planck18_bkg) @@ -110,17 +111,19 @@ end +@doc raw""" + transfer_function(k, [bkg_cosmo::BkgCosmology], with_baryons = true) + +If bkg_cosmo is given computes the transfer function for + """ - transfer_function -""" -function transfer_function(k::Real, cosmo::BkgCosmology{<:Real}; with_baryons::Bool = true) - p = EH98(cosmo) +function transfer_function(k::Real, bkg_cosmo::FLRW{<:Real}; with_baryons::Bool = true) + p = EH98(bkg_cosmo) tc = transfer_cdm(k, p) tb = transfer_baryons(k, p) return with_baryons ? p.Ω_b0_h2 / p.Ω_m0_h2 * tb + p.Ω_χ0_h2/p.Ω_m0_h2 * tc : p.Ω_χ0_h2 / p.Ω_m0_h2 * tc end - function transfer_function(k::Real, p::EH98 = EH98_planck18; with_baryons::Bool = true) tc = transfer_cdm(k, p) tb = transfer_baryons(k, p) @@ -129,8 +132,8 @@ end ## Define here the trivial transfer function struct TrivialTF <: TransferFunctionModel end -TrivialTF(cosmo::BkgCosmology{<:Real}) = TrivialTF() -transfer_function(k::Real, p::TrivialTF; with_baryons::Bool = true) = 1.0 +TrivialTF(bkg_cosmo::BkgCosmology) = TrivialTF() +transfer_function(k::Real, p::TrivialTF; with_baryons::Bool = true) = 1 ## Other transfer functions can be implemented down here