Using Cosmology.jl

Constructing Cosmology Instances

To create an instance of a cosmology that you can use to do computations, you can use the convenience constructor cosmology. By default this constructor will use cosmological parameters from the Planck 2018 results; no arguments are necessary unless you want to override the defaults.

using Cosmology
c = cosmology()

Cosmology.FlatLCDM{Float64, 3}(0.6766, 0.6888463055445425, 0.30966, 0.04897, 2.7255, 3.046, (0.0, 0.0, 0.06))

Alternatively, you can use one of the pre-constructed instances, like Planck18 or WMAP9. These are not exported and so must be imported as

import Cosmology: Planck18
Planck18

Cosmology.FlatLCDM{Float64, 3}(0.6766, 0.6888463055445425, 0.30966, 0.04897, 2.7255, 3.046, (0.0, 0.0, 0.06))

import Cosmology: WMAP9
WMAP9

Cosmology.FlatLCDM{Float64, 3}(0.6932, 0.7134367388797183, 0.2865, 0.04628, 2.725, 3.04, (0.0, 0.0, 0.0))

Tip

cosmology is type-unstable because it returns different concrete subtypes of AbstractCosmology depending on the parameters you give it. This results in a significant performance hit; for example, c=cosmology() takes ~300 ns, while FlatLCDM(c.h,c.Ω_Λ,c.Ω_m,c.Ω_b,c.Tcmb0,c.Neff,c.m_nu) takes ~1 ns. If you want to create many cosmology instances very quickly, it is recommended that you use the base constructors. However, these basic constructors do not include much in the way of argument checking; for example, see the section below on neutrino handling.

import BenchmarkTools: @benchmark
import Cosmology: FlatLCDM
@benchmark cosmology(h=0.6766,OmegaK=0.0,OmegaM=0.30966,OmegaB=0.04897,OmegaG=nothing,Tcmb0=2.7255,w0=-1,wa=0,N_eff=3.046,m_ν=(0.0,0.0,0.06))

BenchmarkTools.Trial: 10000 samples with 226 evaluations.
 Range (min … max):  328.765 ns …  18.780 μs  ┊ GC (min … max): 0.00% … 96.54%
 Time  (median):     355.757 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   390.646 ns ± 672.642 ns  ┊ GC (mean ± σ):  6.78% ±  3.86%

  ▁▄▃  ▆█▆▅▇▇▅▄▂▁ ▁                                ▁▁           ▂
  ████▇██████████████▆█▇▇▇▇█▇█▇▇▇▇▆▆▆▄▅▅▅▄▅▅▄▄▆▇███████▇▇█▇▆▆▇▇ █
  329 ns        Histogram: log(frequency) by time        511 ns <

 Memory estimate: 336 bytes, allocs estimate: 16.

@benchmark FlatLCDM($0.6766, $0.6888463055445425, $0.30966, $0.04897, $2.7255, $3.046, $(0.0,0.0,0.06))

BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
 Range (min … max):  4.200 ns … 23.100 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     4.400 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   4.367 ns ±  0.314 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

                              ▅                            █  
  ▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▂
  4.2 ns         Histogram: frequency by time         4.4 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

Calling Methods

With a cosmological type constructed, we can use it to call methods. For example, to calculate the critical density at a redshift of 1,

ρ_c(c,1.0)

2.7333137424233683e-29 g cm^-3

We can transform this to different units as well,

import Unitful as u
import UnitfulAstro as ua
ρ_c(ua.Msun/ua.kpc^3,c,1.0)

403.86406308531053 M⊙ kpc^-3

Methods with return values that have attached units should also have this conversion interface. See the section on methods for a full list of defined methods.

Neutrinos

This package supports massive neutrinos by assuming that the total mass in neutrinos is split between the species. For greatest efficiency you should pass m_ν to cosmology or one of the basic constructors as an NTuple{N,T}; that is, a tuple where all the elements are of the same concrete type. Internally they are stored in eV but with the units stripped, so a good argument would be m_ν=(0.0,0.0,0.06) for 3 neutrino species with one massive species with mass 0.06 eV. You can provide neutrino masses to cosmology and the basic constructors with energy or mass units as well;

cosmology(m_ν=(0.0,0.0,0.06) .* u.eV)

Cosmology.FlatLCDM{Float64, 3}(0.6766, 0.6888463055445425, 0.30966, 0.04897, 2.7255, 3.046, (0.0, 0.0, 0.06))

FlatLCDM(0.6766, 0.6888463055445425, 0.30966, 0.04897, 2.7255, 3.046, (0.0,0.0,0.06) .* u.eV)

Cosmology.FlatLCDM{Float64, 3}(0.6766, 0.6888463055445425, 0.30966, 0.04897, 2.7255, 3.046, (0.0, 0.0, 0.06))

There is currently some weirdness with how massless neutrinos are dealt with; for now, I recommend that you make the length of the m_ν iterable you provide equal to the result of Cosmology.n_nu(N_eff) for the effective number of neutrino species you choose. For example, if N_eff=3.046 then Cosmology.n_nu(N_eff)=3 and your m_ν iterable should have length 3. This is done for you in the pre-constructed instances. A warning will be issued by cosmology if length(m_ν) != Cosmology.n_nu(N_eff) and !iszero(N_eff), but none of the basic constructors contain such checks. The implementation of massive neutrinos is open to change.

Tip

Inclusion of massive neutrinos is expensive. For example, for the default massive neutrino parameters c=cosmology(), the evaluation of E(c, 0.8) takes 114.613 ns, while E( cosmology(m_ν=(0.0,),N_eff=3.046), 0.8) takes 6.986 ns and E( cosmology(m_ν=(0.0,),N_eff=0), 0.8) takes 6.095 ns. This makes a significant difference in methods that involve integrations (e.g., comoving_radial_dist). If speed is a concern, consider if you can neglect neutrinos for your calculation.

c_massivenu = cosmology()
c_masslessnu = cosmology(m_ν=(0.0,),N_eff=3.046)
c_nonu = cosmology(m_ν=(0.0,),N_eff=0.0)

┌ Warning: Length of `m_ν` should typically be equal to `Int(floor(N_eff))` but is not. Did you make a mistake?
└ @ Cosmology ~/work/Cosmology.jl/Cosmology.jl/src/Cosmology.jl:339

@benchmark E($c_massivenu,$0.8)

BenchmarkTools.Trial: 10000 samples with 923 evaluations.
 Range (min … max):  116.361 ns … 144.313 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     116.579 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   116.793 ns ±   1.038 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▅██▁                                                         ▂
  █████▁▁▁▃▃▁▁▁▁▁▁▅▆▁▇▇▇▆▆▁▅▄▄▄▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▃▅▇▇▁███▇▇▆ █
  116 ns        Histogram: log(frequency) by time        122 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

@benchmark E($c_masslessnu,$0.8)

BenchmarkTools.Trial: 10000 samples with 996 evaluations.
 Range (min … max):  22.691 ns … 40.463 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     22.791 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   22.861 ns ±  0.510 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

                               █                               
  ▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▂
  22.7 ns         Histogram: frequency by time        22.9 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

@benchmark E($c_nonu,$0.8)

BenchmarkTools.Trial: 10000 samples with 997 evaluations.
 Range (min … max):  22.066 ns … 39.920 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     22.267 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   22.278 ns ±  0.463 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

                                         █                     
  ▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▇▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂ ▂
  22.1 ns         Histogram: frequency by time        22.4 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

Integrated Packages

Cosmology.jl is envisioned as the base for a collection of packages for things like cosmological power spectra, growth factors, halo mass functions, etc. As they become available, some examples will be shown here and on the integrated packages page.