Added most functions for hubble's constant and distance calculation.

.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:

.github/workflows/CI.yml

@@ -16,7 +16,7 @@ jobs:
         version:
           - '1.6'
           - '1.9'
-          - 'nightly'
+          # - 'nightly'
         os:
           - ubuntu-latest
         arch:

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"

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
+

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,
+        χ
 
 
 

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
+