Getting Started with CosmicPy

First, start an ipython sesion in pylab mode:

$ ipython --pylab

This will automatically load numpy and matplolib and enable the interactive plotting environment.

Once ipython is running, you can import the cosmicpy package:

In [1]: from cosmicpy import *

Create and modify a cosmology

The next step is to create a cosmology object. Let us start with a default cosmology:

In [2]: cosmo = cosmology()

To get a summary of the cosmological parameters defined in this object, you can use the print() fucntion:

In [3]: print cosmo
FLRW Cosmology with the following parameters: 
    h:        0.7 
    Omega_b:  0.045 
    Omega_m:  0.25 
    Omega_de: 0.75 
    w0:       -0.95 
    wa:       0.0 
    n:        1.0 
    tau:      0.09 
    sigma8:   0.8

All the cosmological parameters are accessible as attributes of the cosmology object. For instance, to check the value of the matter density today \(\Omega_m\) run:

In [4]: cosmo.Omega_m
Out[4]: 0.25

The values for the cosmological parameters can be modified which will automatically trigger an update of all the derived quantities.

To illustrate this, let’s see the value of the sound horizon at recombination \(s_h\) with the default cosmology, update \(\Omega_m\) and check that the value of \(s_h\) has changed:

In [5]: cosmo.sh_r
Out[5]: 105.27381988350771

In [6]: cosmo.Omega_m = 0.3

In [7]: cosmo.sh_r
Out[7]: 100.59308203228834

The cosmological parameters can also be set at the creation of the cosmology object by specifying their values to the constructor. For instance, to create a cosmology with \(\Omega_m = 0.2\) and \(\Omega_{de} = 0.8\):

In [8]: cosmo2 = cosmology(Omega_m=0.2, Omega_de=0.8)

In [9]: print cosmo2
FLRW Cosmology with the following parameters: 
    h:        0.7 
    Omega_b:  0.045 
    Omega_m:  0.2 
    Omega_de: 0.8 
    w0:       -0.95 
    wa:       0.0 
    n:        1.0 
    tau:      0.09 
    sigma8:   0.8

Compute distances

Now that we have seen how to create a cosmology object with the desired cosmological parameters, let’s use it to compute some distances.

cosmology offers a number of different methods to compute various cosmological distances:

cosmology.a2chi(a) Radial comoving distance in [Mpc/h] for a given scale factor.
cosmology.f_k(a) Transverse comoving distance in [Mpc/h] for a given scale factor.
cosmology.d_A(a) Angular diameter distance in [Mpc/h] for a given scale factor.

These methods all compute a distance as a function of the scale factor \(a\).

It is often usefull to work with redshifts instead of scale factors. The conversion between the two quantities is performed by the 2 utility functions:

cosmicpy.utils.z2a(z) converts from redshift to scale factor
cosmicpy.utils.a2z(a) converts from scale factor to redshift

Here is an example to compute the radial comoving distance for an array of redshifts between 0 and 10:

In [10]: z = arange(0,10,0.01)

In [11]: figure();

# The z2a function is used to convert z to scale factors
In [12]: plot(z,cosmo.a2chi(z2a(z))); grid();

In [13]: title('Radial comoving distance as a function of redshift');

Compute the matter power spectrum

cosmology implements several methods to compute the linear and non linear matter power spectrum:

cosmology.pk_lin(k[, a]) Computes the linear matter power spectrum.[, a, nl_type]) Computes the full non linear matter power spectrum.

By default, these two functions will compute the power spectrum at \(a=1\) and using the Eisenstein & Hu transfer function with oscillations. The non-linear correction implemented in is computed based on [SmithPeacockJenkins+03].

Here is an example to compute these two power spectra:

# Creating a vector of wavenumbers
In [14]: k = logspace(-3,0,100)

In [15]: loglog(k,, label='Smith et al. $P(k)$');

In [16]: loglog(k,cosmo.pk_lin(k), label='Linear $P(k)$');

In [17]: legend(); xlabel('k [h/Mpc]');

Both power spectra can be computed with or without baryonic oscillations by setting the type keyword to either:

  • ‘eisenhu_osc’ : Include BAOs (default behavior)
  • ‘eisenhu’ : Without BAOs

This example shows the matter power spectrum, with non-linear corrections, with and without BAO wiggles:

In [18]: loglog(k,, type='eisenhu'), label='eisenhu');

In [19]: loglog(k,, type='eisenhu_osc'), label='eisenhu_osc');

In [20]: legend(); xlabel('k [h/Mpc]'); title('Matter power spectrum');

Finally, the power spectrum can be computed at different scale factors:

# Create an array of scale factors corresponding to z=0, 1 ,2 and 5
In [21]: a = z2a(array([0, 1, 2, 5]))

In [22]: loglog(k,,a));

In [23]: xlabel('k [h/Mpc]'); title('Matter power spectrum at z=0, 1, 2, 5');


[SmithPeacockJenkins+03]R. E. Smith, J. A. Peacock, A. Jenkins, S. D. M. White, C. S. Frenk, F. R. Pearce, P. A. Thomas, G. Efstathiou, and H. M. P. Couchman. Stable clustering, the halo model and non-linear cosmological power spectra. MNRAS, 341:1311–1332, June 2003. arXiv:astro-ph/0207664, doi:10.1046/j.1365-8711.2003.06503.x.