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 250 evaluations per sample.
 Range (min … max):  298.964 ns …  35.423 μs  ┊ GC (min … max): 0.00% … 98.60%
 Time  (median):     313.348 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   337.884 ns ± 460.263 ns  ┊ GC (mean ± σ):  5.67% ±  5.66%

   ▃▅▆▆█▇▄▂▁▁          ▁▃▃▂    ▁                                ▂
  ███████████▇▅▆▄▃▄▅▅▇███████▇██▇▆▆▅▅▃▃▅▅▃▅▆▄▃▄▃▃▁▄▃▄▄▁▁▄▃▃▁▁▁▃ █
  299 ns        Histogram: log(frequency) by time        454 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 per sample.
 Range (min … max):  3.396 ns … 21.811 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     3.407 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   3.452 ns ±  0.583 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

         █                                                    
  ▂▁▁▁▁▁▁█▁▁▁▁▁▁▇▅▁▁▁▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▂▁▁▁▁▁▁▂▂▁▁▁▁▁▁▃▁▁▁▁▁▁▂ ▂
  3.4 ns         Histogram: frequency by time        3.48 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.7333137424245774e-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.8640630854892 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 934 evaluations per sample.
 Range (min … max):  106.742 ns … 199.851 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     106.829 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   107.892 ns ±   3.257 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

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

 Memory estimate: 0 bytes, allocs estimate: 0.

@benchmark E($c_masslessnu,$0.8)

BenchmarkTools.Trial: 10000 samples with 997 evaluations per sample.
 Range (min … max):  18.339 ns … 41.391 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     19.315 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   19.448 ns ±  1.395 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

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

 Memory estimate: 0 bytes, allocs estimate: 0.

@benchmark E($c_nonu,$0.8)

BenchmarkTools.Trial: 10000 samples with 998 evaluations per sample.
 Range (min … max):  17.919 ns … 107.407 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     18.120 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   18.357 ns ±   1.680 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▃█▃                                                          ▁
  ███▁▄▄█▅▄▃▁▃▁▁▄▅▅▄▁▁▁▁▁▁▁▁▃▃▁▁▁▁▁▃▁▁▁▁▃▃▁▁▁▁▁▁▁▁▁▁▁▁▃▁▄▆▇▇█▇ █
  17.9 ns       Histogram: log(frequency) by time      26.9 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.