cosmo_compare.py 3.89 KB
Newer Older
1
"""cosmo_compare -- Comparison values to test cosmo.el
2

Francesco Montanari's avatar
Francesco Montanari committed
3 4 5 6 7 8 9 10 11 12 13 14 15 16
Copyright (C) 2017 Francesco Montanari

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
17
"""
Francesco Montanari's avatar
Francesco Montanari committed
18

19
from astropy.cosmology import FlatLambdaCDM
20
from astropy.cosmology import LambdaCDM
Francesco Montanari's avatar
Francesco Montanari committed
21
import numpy as np
22

23

Francesco Montanari's avatar
Francesco Montanari committed
24 25 26 27
def sexp_print(reds, vals, name=' '):
    """Print redshift list and respective values list as sexp."""
    vals_str = '(' + ' '.join(map(str, vals)) + ')'
    print('{}:\n  {}'.format(name, vals_str))
28

Francesco Montanari's avatar
Francesco Montanari committed
29 30
    reds_str = '(' + ' '.join(map(str, reds)) + ')'
    print('redshifts:\n  {}'.format(reds_str))
31

Francesco Montanari's avatar
Francesco Montanari committed
32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
    print('')


def print_curves(reds, cosmo):
    """Print cosmological functions given redshifts and a cosmology."""
    sexp_print(reds, cosmo.efunc(reds), name='efunc')
    sexp_print(reds, cosmo.inv_efunc(reds), name='inv-efunc')
    sexp_print(reds, cosmo.H(reds).value,
               name='hubble [km/s/Mpc]')
    sexp_print(['-'], [cosmo.hubble_distance.value],
               name='hubble-distance [Mpc]')
    sexp_print(['-'], [cosmo.hubble_time.value],
               name='hubble-time [Gyr]')
    sexp_print(reds, cosmo.comoving_distance(reds).value,
               name='los-comoving-distance [Mpc]')
    sexp_print(reds, cosmo.comoving_transverse_distance(reds).value,
               name='transverse-comoving-distance [Mpc]')
    sexp_print(reds, cosmo.luminosity_distance(reds).value,
               name='luminosity-distance [Mpc]')
    sexp_print(reds, cosmo.angular_diameter_distance(reds).value,
               name='angular-diameter-distance [Mpc]')
53 54
    sexp_print(reds, cosmo.comoving_volume(reds).value,
               name='comoving-volume [Mpc^3]')
55 56 57 58
    sexp_print(reds, cosmo.lookback_time(reds).value,
               name='lookback-time [Gyr]')
    sexp_print(reds, cosmo.age(reds).value,
               name='age [Gyr]')
Francesco Montanari's avatar
Francesco Montanari committed
59 60 61 62 63 64 65 66 67


def print_open_cosmo():
    """Print cosmological functions for an open cosmology."""
    H0 = 70.
    Om0 = 0.31
    Ode0 = 0.7
    cosmo = LambdaCDM(H0=H0, Om0=Om0, Ode0=Ode0)

68 69 70
    print('===============')
    print('Open cosmology:')
    print('===============')
Francesco Montanari's avatar
Francesco Montanari committed
71 72 73
    print('H0={} [km/s/Mpc], Om0={}, Ode0={}, Orel0={}, Ok0={}\n'
          .format(H0, Om0, Ode0, cosmo.Ogamma(0.)+cosmo.Onu(0.),
                  cosmo.Ok(0.)))
74 75

    redshifts = (0., 0.1, 10., 1000.)
Francesco Montanari's avatar
Francesco Montanari committed
76 77 78
    print_curves(redshifts, cosmo)


79 80 81 82 83 84 85
def print_close_cosmo():
    """Print cosmological functions for a close cosmology."""
    H0 = 70.
    Om0 = 0.27
    Ode0 = 0.7
    cosmo = LambdaCDM(H0=H0, Om0=Om0, Ode0=Ode0)

86 87 88
    print('================')
    print('Close cosmology:')
    print('================')
89 90 91 92 93 94 95 96 97 98 99 100 101 102
    print('H0={} [km/s/Mpc], Om0={}, Ode0={}, Orel0={}, Ok0={}\n'
          .format(H0, Om0, Ode0, cosmo.Ogamma(0.)+cosmo.Onu(0.),
                  cosmo.Ok(0.)))

    redshifts = (0., 0.1, 10., 1000.)
    print_curves(redshifts, cosmo)


def print_flat_cosmo():
    """Print cosmological functions for a flat cosmology."""
    H0 = 70.
    Om0 = 0.3
    cosmo = FlatLambdaCDM(H0=H0, Om0=Om0, Tcmb0=0.)

103 104 105
    print('===============')
    print('Flat cosmology:')
    print('===============')
106 107 108 109 110 111 112 113
    print('H0={} [km/s/Mpc], Om0={}, Ode0={}, Orel0={}, Ok0={}\n'
          .format(H0, Om0, cosmo.Ode(0.), cosmo.Ogamma(0.)+cosmo.Onu(0.),
                  cosmo.Ok(0.)))

    redshifts = (0., 0.1, 10., 1000.)
    print_curves(redshifts, cosmo)


Francesco Montanari's avatar
Francesco Montanari committed
114 115
if __name__ == '__main__':
    print_open_cosmo()
116 117
    print_close_cosmo()
    print_flat_cosmo()