Update few errors and add time functions

.travis.yml

@@ -4,7 +4,7 @@ notifications:
   email: false
 julia:
   - 1.6
-  - 1.9
+  - 1.10
   - nightly
 os:
   - linux

Project.toml

@@ -8,13 +8,15 @@ DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 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"
 
 [extras]
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 
 [targets]
 test = ["Test"]

src/basic_cos.jl

@@ -84,12 +84,8 @@ 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)
+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)
@@ -117,28 +113,31 @@ 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
+redshift(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_by_H0(c::AbstractCosmology, z) = (scal = a(z); H_a2_by_H0(c, scal) / scal^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]
+χ(c::AbstractCosmology, z::Real, ::Nothing; kws...) = QuadGK.quadgk(a->1 / H_a2_by_H0(c, a), a(z), 1; kws...)[1]
+χ(c::AbstractCosmology, z₁::Real, z₂::Real; kws...) = QuadGK.quadgk(a->1 / H_a2_by_H0(c, a), a(z₂), a(z₁); kws...)[1]
 # Times
-
-
+η0(c::AbstractCosmology) = 3.262*9.46*1e18 / (1e2*c.h)# sec
+hubble_time(c::AbstractCosmology, z) = η0(c) / H_by_H0(c, z)
+T(c::AbstractCosmology, a0, a1; kws...) = QuadGK.quadgk(x->x / H_a2_by_H0(c, x), a0, a1; kws...)[1]
+age(c::AbstractCosmology, z; kws...) = η0(c) * T(c, 0, a(z); kws...)
+lookback_time(c::AbstractCosmology, z; kws...) = η0(c) * T(c, a(z), 1; kws...)
 # 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
+# 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
 
 

src/cosmic.jl

@@ -2,18 +2,7 @@ module cosmic
 
 include("basic_cos.jl")
 
-export a, z, cosmology, Mpc_to_km, sec_to_year, H, my_f, my_g,
-        hubble_distance,
-        χ
-
-
-
-
-
-
-
-
-
-
+export a, redshift, cosmology, Mpc_to_km, sec_to_year, H, my_f, my_g,
+        hubble_distance, χ, hubble_time, T, age, lookback_time
 
 end#  module end

test/runtests.jl

@@ -4,12 +4,15 @@ using Test
 @testset "cosmic.jl" begin
     @test my_f(2,1) == 11
     @test my_g(2) == 4
-    @test z(2.9e-4) ≈ 3447.27586
+    @test redshift(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
+using Plots, LaTeXStrings
 
+c = cosmology()
+z_vals = range(0, 2.5, length=1000)
+H_1z(z) = H(c,z)/(1+z)
+H_vals = H_1z.(z_vals)
+plot(z_vals,H_vals,lw=2.5, label="")
+xlabel!(L"Redshift ($z$)")
+ylabel!(L"$\frac{H(z)}{1+z}$")