From 73bc4e28851abc6951e3257ffafbbfdabc6e1cc1 Mon Sep 17 00:00:00 2001 From: Kazi Abu Rousan Date: Thu, 18 Jan 2024 19:18:29 +0530 Subject: [PATCH] Added most functions for hubble's constant and distance calculation. --- .appveyor.yml | 2 +- .github/workflows/CI.yml | 2 +- Project.toml | 3 +- src/basic_cos.jl | 156 +++++++++++++++++++++++++++++++++++++++ src/cosmic.jl | 17 +---- src/test_f.jl | 1 - src/thermo.jl | 0 test/runtests.jl | 12 ++- 8 files changed, 173 insertions(+), 20 deletions(-) create mode 100644 src/basic_cos.jl delete mode 100644 src/test_f.jl create mode 100644 src/thermo.jl diff --git a/.appveyor.yml b/.appveyor.yml index 249ceab..d2e12d7 100644 --- a/.appveyor.yml +++ b/.appveyor.yml @@ -3,7 +3,7 @@ environment: matrix: - julia_version: 1.0 - julia_version: 1.9 - - julia_version: nightly + # - julia_version: nightly platform: - x64 cache: diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index c073b54..0837e9a 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -16,7 +16,7 @@ jobs: version: - '1.6' - '1.9' - - 'nightly' + # - 'nightly' os: - ubuntu-latest arch: diff --git a/Project.toml b/Project.toml index 8067728..6d639af 100644 --- a/Project.toml +++ b/Project.toml @@ -5,11 +5,10 @@ version = "1.0.0-DEV" [deps] DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" -Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" [compat] julia = "1" diff --git a/src/basic_cos.jl b/src/basic_cos.jl new file mode 100644 index 0000000..3d9bf7b --- /dev/null +++ b/src/basic_cos.jl @@ -0,0 +1,156 @@ +using QuadGK + + + + +abstract type AbstractCosmology end +abstract type AbstractClosedCosmology <: AbstractCosmology end +abstract type AbstractFlatCosmology <: AbstractCosmology end +abstract type AbstractOpenCosmology <: AbstractCosmology end + + +# For flat ΛCDM +struct FlatLCDM{T <: Real} <: AbstractFlatCosmology + h::T + Ω_Λ::T + Ω_m::T + Ω_r::T +end +FlatLCDM(h::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real) = FlatLCDM(promote(float(h), float(Ω_Λ), float(Ω_m), float(Ω_r))...) +H_a2_by_H0(c::FlatLCDM, a) = sqrt(c.Ω_r + c.Ω_m * a + c.Ω_Λ * a^4) +# mostly we will use this + +# For just an extension later will work on this +struct ClosedLCDM{T <: Real} <: AbstractClosedCosmology + h::T + Ω_k::T + Ω_Λ::T + Ω_m::T + Ω_r::T +end +ClosedLCDM(h::Real, Ω_k::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real) = ClosedLCDM(promote(float(h), float(Ω_k), float(Ω_Λ), float(Ω_m),float(Ω_r))...) + +struct OpenLCDM{T <: Real} <: AbstractOpenCosmology + h::T + Ω_k::T + Ω_Λ::T + Ω_m::T + Ω_r::T +end +OpenLCDM(h::Real, Ω_k::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real) = OpenLCDM(promote(float(h), float(Ω_k), float(Ω_Λ), float(Ω_m),float(Ω_r))...) + +function H_a2_by_H0(c::Union{ClosedLCDM,OpenLCDM}, a) + a2 = a * a + return sqrt(c.Ω_r + c.Ω_m * a + (c.Ω_k + c.Ω_Λ * a2) * a2) +end +# end here. + +#For understanding w0 and wa see the paper. +#https://doi.org/10.1016/j.aop.2023.169244 +for c in ("Flat", "Open", "Closed") + name = Symbol("$(c)WCDM") + @eval begin + struct $(name){T <: Real} <: $(Symbol("Abstract$(c)Cosmology")) + h::T + Ω_k::T + Ω_Λ::T + Ω_m::T + Ω_r::T + w0::T + wa::T + end + function $(name)(h::Real, Ω_k::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real, + w0::Real, wa::Real) + $(name)(promote(float(h), float(Ω_k), float(Ω_Λ), float(Ω_m), + float(Ω_r), float(w0), float(wa))...) + end + end +end +#For a more appropriate modle +function WCDM(h::Real, Ω_k::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real, w0::Real, wa::Real) + if Ω_k < 0 + ClosedWCDM(h, Ω_k, Ω_Λ, Ω_m, Ω_r, w0, wa) + elseif Ω_k > 0 + OpenWCDM(h, Ω_k, Ω_Λ, Ω_m, Ω_r, w0, wa) + else + FlatWCDM(h, Ω_k, Ω_Λ, Ω_m, Ω_r, w0, wa) + end +end + +function H_a2_by_H0(c::Union{FlatWCDM,ClosedWCDM,OpenWCDM}, a) + ade = exp((1 - 3 * (c.w0 + c.wa)) * log(a) + 3 * c.wa * (a - 1)) + return sqrt(c.Ω_r + (c.Ω_m + c.Ω_k * a) * a + c.Ω_Λ * ade) +end + +# Final Cosmology Function +# All predefined values are taken from Baumann's book table-2.1 +function cosmology(;h = 0.6774,Neff = 3.04, + Ωk = 0, + Ωm = 0.3089, + Ωr = nothing, + Tcmb = 2.7255, + w0 = -1,wa = 0) + if Ωr === nothing + Ωγ = 4.48131e-7 * Tcmb^4 / h^2 + Ων = Neff * Ωγ * (7 / 8) * (4 / 11)^(4 / 3) + Ωr = Ωγ + Ων + end + + ΩΛ = 1 - Ωk - Ωm - Ωr + + if !(w0 == -1 && wa == 0) + return WCDM(h, Ωk, ΩΛ, Ωm, Ωr, w0, wa) + end + + if Ωk < 0 + return ClosedLCDM(h, Ωk, ΩΛ, Ωm, Ωr) + elseif Ωk > 0 + return OpenLCDM(h, Ωk, ΩΛ, Ωm, Ωr) + else + return FlatLCDM(h, ΩΛ, Ωm, Ωr) + end +end +#-------------------------------------------------------------------------------- +#Conversion +Mpc_to_km(mp_km) = mp_km * 3.262e6 * 9.461e12 +sec_to_year(sec_ye) = 3.171e-8*sec_ye +#-------------------------------------------------------------------------------- +# Scale Factor +a(z) = 1/(1+z)#Sacle factor as a function of z +z(a) = 1/a - 1# Inverse function +#-------------------------------------------------------------------------------- +# Hubble Rate +H_by_H0(c::AbstractCosmology, z) = (a = a(z); H_a2_by_H0(c, a) / a^2) +H(c::AbstractCosmology, z) = 100 * c.h * H_by_H0(c, z)# in km /s / Mpc +#-------------------------------------------------------------------------------- +# Distance +χ0(c::AbstractCosmology) = 2997.92458 / c.h# in Mpc +hubble_distance(c::AbstractCosmology, z) = χ0(c) / H_by_H0(c, z)# 1/H -> hubble radius +χ(c::AbstractCosmology, z::Real, ::Nothing; kws...) = QuadGK.quadgk(a->1 / H_a2_by_H0(c, a), scale_factor(z), 1; kws...)[1] +χ(c::AbstractCosmology, z₁::Real, z₂::Real; kws...) = QuadGK.quadgk(a->1 / H_a2_by_H0(c, a), scale_factor(z₂), scale_factor(z₁); kws...)[1] +# Times + + +# Constants all are in MeV and C = 1, ħ = 1, kᵦ = 1 +H0 = (67.74/(3.26*9.5))*1e-18 +Ωm = 0.32 +ΩΛ = 0.68 +Ωγ = 5.35e-5 +Ων = 3.64e-5 +Ωr = Ωγ + Ων +zeq = 3395 +# All constants are from Baumann's book + + + +function inte_for_age(a) + deno = √(Ωr/a^2 + Ωm/a + ΩΛ*a^2) + return 1/deno +end + + +age(a) = QuadGK.quadgk(inte_for_age, 0, a)[1]/H0 # In sec +my_f(x,y) = 2x+7y + +my_g(x) = x^2 + diff --git a/src/cosmic.jl b/src/cosmic.jl index e286745..5267176 100644 --- a/src/cosmic.jl +++ b/src/cosmic.jl @@ -1,19 +1,10 @@ module cosmic -include("test_f.jl") - -export my_f - - - - - - - - - - +include("basic_cos.jl") +export a, z, cosmology, Mpc_to_km, sec_to_year, H, my_f, my_g, + hubble_distance, + χ diff --git a/src/test_f.jl b/src/test_f.jl deleted file mode 100644 index a62f716..0000000 --- a/src/test_f.jl +++ /dev/null @@ -1 +0,0 @@ -my_f(x,y) = 2x+3y \ No newline at end of file diff --git a/src/thermo.jl b/src/thermo.jl new file mode 100644 index 0000000..e69de29 diff --git a/test/runtests.jl b/test/runtests.jl index 40f699d..f1282dd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,6 +2,14 @@ using cosmic using Test @testset "cosmic.jl" begin - @test my_f(2,1) == 7 - @test my_f(2,3) == 13 + @test my_f(2,1) == 11 + @test my_g(2) == 4 + @test z(2.9e-4) ≈ 3447.27586 end + +Tcmb = 2.7260 +h = 0.740 +OmegaG = 4.48131e-7 * Tcmb^4 / h^2 +OmegaN = Neff * OmegaG * (7 / 8) * (4 / 11)^(4 / 3) +OmegaR = OmegaG + OmegaN +