Halo circular velocity function from Maller & Bullock 2004

.github/CODEOWNERS

@@ -1,3 +1,5 @@
 skypy @skypyproject/skypy-infrastructure
 skypy/galaxy @skypyproject/galaxy-developers
+skypy/halos @skypyproject/halo-developers
+skypy/power_spectrum @skypyproject/power-spectrum-developers
 CODE_OF_CONDUCT.md @skypyproject/skypy-board

README.rst

@@ -8,6 +8,8 @@ This package contains methods for modelling the Universe, galaxies and the Milky
 Way. Also included are methods for generating observed data.
 
 * Galaxies_: morphology, luminosity and redshift distributions
+* Halo_ and subhalo mass distributions
+* `Power Spectra`_ using CAMB and Halofit
 * Pipelines_ to generate populations of astronomical objects
 
 The full list of features can be found in the `SkyPy Documentation`_.
@@ -16,6 +18,8 @@ If you use SkyPy for work or research presented in a publication please follow
 our `Citation Guidelines`_.
 
 .. _Galaxies: https://skypy.readthedocs.io/en/latest/galaxies.html
+.. _Halo: https://skypy.readthedocs.io/en/latest/halos/index.html
+.. _Power Spectra: https://skypy.readthedocs.io/en/latest/power_spectrum/index.html
 .. _Pipelines: https://skypy.readthedocs.io/en/latest/pipeline/index.html
 .. _SkyPy Documentation: https://skypy.readthedocs.io/en/latest/
 .. _Citation Guidelines: CITATION.rst

docs/halos/examples/halos.yml

@@ -0,0 +1,24 @@
+cosmology: !astropy.cosmology.default_cosmology.get 
+power_spectrum: !skypy.power_spectrum.eisenstein_hu
+  wavenumber: !numpy.logspace [-3, 1, 100]
+  A_s: 2.1982e-09
+  n_s: 0.969453
+  cosmology: $cosmology
+sheth-tormen: !skypy.halos.mass.sheth_tormen
+  m_min: 1.0E+10
+  m_max: 1.0E+13
+  resolution: 10000
+  wavenumber: !numpy.logspace [-3, 1, 100]
+  power_spectrum: $power_spectrum
+  growth_function: 1.0
+  cosmology: $cosmology
+  size: 10000
+press-schechter: !skypy.halos.mass.press_schechter
+  m_min: 1.0E+10
+  m_max: 1.0E+13
+  resolution: 10000
+  wavenumber: !numpy.logspace [-3, 1, 100]
+  power_spectrum: $power_spectrum
+  growth_function: 1.0
+  cosmology: $cosmology
+  size: 10000

docs/halos/index.rst

@@ -0,0 +1,113 @@
+*********************************
+Dark Matter Halos (`skypy.halos`)
+*********************************
+
+You can reproduce figure 2 in Sheth and Tormen 1999
+and plot the collapse functions for different halo models. For this, you can use
+a python script, for example.
+
+.. plot::
+   :include-source: true
+   :nofigs:
+   :context: close-figs
+
+    import numpy as np
+    from astropy.cosmology import Planck15
+    from skypy.power_spectrum import eisenstein_hu
+    from skypy.halos.mass import _sigma_squared
+
+    # Power spectrum and amplitude of perturbations at redshift 0
+    growth_0 = 1.0
+    A_s, n_s = 2.1982e-09, 0.969453
+    k = np.logspace(-3, 1, num=1000, base=10.0)
+    mass = 10**np.arange(9.0, 15.0, 0.1)
+
+    pk0 = eisenstein_hu(k, A_s, n_s, Planck15, kwmap=0.02, wiggle=True)
+    sigma = np.sqrt(_sigma_squared(mass, k, pk0, growth_0, Planck15))
+
+    # Collapse functions
+    from skypy.halos import mass
+
+    params_ellipsoidal = (0.3, 0.7, 0.3, 1.686)
+
+    ST = mass.sheth_tormen_collapse_function(sigma)
+    PS = mass.press_schechter_collapse_function(sigma)
+    EM = mass.ellipsoidal_collapse_function(sigma, params=params_ellipsoidal)
+
+
+.. plot::
+   :include-source: false
+   :context: close-figs
+
+    import matplotlib.pyplot as plt
+
+    delta_c = 1.69
+    nu = np.square(delta_c / sigma)
+
+    # plot different collapse functions
+    plt.loglog(nu, ST, label='Sheth-Tormen')
+    plt.loglog(nu, PS, label='Press-Schechter')
+    plt.loglog(nu, EM, label='Ellipsoidal')
+
+    # axes labels and title
+    plt.xlabel(r'$\nu \equiv (\delta_c / \sigma)^2$')
+    plt.ylabel(r'$f_c(\nu)$')
+    plt.title('Collapse function')
+
+    # show plot labels
+    plt.legend()
+    plt.show()
+
+
+You can also sample halos using their mass function. For this, you can use a config
+file and run the pipeline, for example.
+
+.. literalinclude:: examples/halos.yml
+   :language: yaml
+
+.. plot::
+   :include-source: false
+   :context: close-figs
+
+    from skypy.pipeline import Pipeline
+    pipeline = Pipeline.read('examples/halos.yml')
+    pipeline.execute()
+
+    # Draw from different halo mass samplers
+    halo_massST = pipeline['sheth-tormen']
+    halo_massPS = pipeline['press-schechter']
+
+    plt.hist(np.log10(halo_massST), histtype='step', label='Sheth-Tormen')
+    plt.hist(np.log10(halo_massPS), histtype='step', label='Press-Schechter')
+
+    # axis label and title
+    plt.xlabel(r'$log(mass)$')
+    plt.title('Halo sampler')
+
+    # show plot labels
+    plt.legend()
+    plt.show()
+
+
+Halos (`skypy.halos`)
+=====================
+
+.. automodule:: skypy.halos
+
+
+Abundance Matching (`skypy.halos.abundance_matching`)
+=====================================================
+
+.. automodule:: skypy.halos.abundance_matching
+
+
+Mass (`skypy.halos.mass`)
+=========================
+
+.. automodule:: skypy.halos.mass
+
+
+Quenching (`skypy.halos.quenching`)
+===================================
+
+.. automodule:: skypy.halos.quenching

docs/index.rst

@@ -32,6 +32,8 @@ Packages
    :maxdepth: 1
 
    galaxies
+   halos/index
+   power_spectrum/index
    utils/index
 
 Pipeline

docs/power_spectrum/examples/power_spectrum.yml

@@ -0,0 +1,14 @@
+cosmology: !astropy.cosmology.default_cosmology.get
+wavenumber: !numpy.logspace [-3, 1, 100]
+eisenstein_hu_wiggle: !skypy.power_spectrum.eisenstein_hu
+  wavenumber: $wavenumber
+  A_s: 2.1982e-09
+  n_s: 0.969453
+  cosmology: $cosmology
+  kwmap: 0.02
+  wiggle: True
+halofit: !skypy.power_spectrum.halofit_smith
+  wavenumber: $wavenumber
+  redshift: 0.0
+  linear_power_spectrum: $eisenstein_hu_wiggle
+  cosmology: $cosmology

docs/power_spectrum/index.rst

@@ -0,0 +1,38 @@
+***************************************
+Power Spectrum (`skypy.power_spectrum`)
+***************************************
+This module contains methods to model the matter power spectrum.
+
+SkyPy provides wrappers to a number of external codes for calculating the
+matter power spectrum, including `~skypy.power_spectrum.camb` and
+`~skypy.power_spectrum.classy`. Here we demonstrate calculating the linear
+matter power spectrum using `~skypy.power_spectrum.eisenstein_hu` and the
+non-linear corrections using `~skypy.power_spectrum.halofit_smith`:
+
+.. literalinclude:: examples/power_spectrum.yml
+   :language: yaml
+
+.. plot::
+   :include-source: false
+   :context: close-figs
+
+    import matplotlib.pyplot as plt
+    from skypy.pipeline import Pipeline
+
+    pipeline = Pipeline.read('examples/power_spectrum.yml')
+    pipeline.execute()
+
+    # Eisenstein and Hu power spectrum and Halofit matter power spectra
+    k = pipeline['wavenumber']
+    power_EH_w = pipeline['eisenstein_hu_wiggle']
+    hf_Smith = pipeline['halofit']
+
+    plt.loglog(k, power_EH_w, label='Eisenstein & Hu')
+    plt.loglog(k, hf_Smith, '--', label='Halofit')
+    plt.xlabel(r'Wavenumber $(1/Mpc)$')
+    plt.ylabel(r'Power spectrum $(Mpc^3)$')
+    plt.legend(frameon=False, loc='lower left');
+
+
+
+.. automodule:: skypy.power_spectrum

setup.cfg

@@ -34,9 +34,11 @@ test =
     speclite>=0.11
 all =
     speclite>=0.11
+    colossus
 docs =
     sphinx-astropy
     matplotlib
+    colossus
     speclite>=0.11
 
 [options.package_data]

skypy/__init__.py

@@ -15,4 +15,6 @@
 __all__ = []
 
 from . import galaxies  # noqa
+from . import halos  # noqa
 from . import pipeline  # noqa
+from . import power_spectrum  # noqa

skypy/halo/properties.py

@@ -0,0 +1,71 @@
+"""Halo properties module.
+
+This module provides methods to add simple properties to halos
+
+Models
+======
+
+.. autosummary::
+   :nosignatures:
+   :toctree: ../api/
+
+   halo_circular_velocity
+
+"""
+
+import numpy as np
+from astropy.units import Unit
+
+__all__ = [
+    'halo_circular_velocity',
+ ]
+
+def halo_circular_velocity(halo_virial_mass, Delta_v, redshift, cosmology):
+    """Halo circular velocity.
+    This function computes the halo circular velocity, setting it 
+    equal to the virial velocity using equation (3) and footnote 2 from [1]_.
+
+    Parameters
+    ----------
+    halo_virial_mass : (nm,) array_like
+    	Array for the virial mass, in units of solar mass.
+    Delta_v : (nm,) array_like
+        The mean overdensity of the halo.
+    redshift : (nm,) array_like
+        The redshift of each halo.
+    cosmology : astropy.cosmology.Cosmology
+        Cosmology object providing methods for the evolution history of
+        omega_matter and omega_lambda with redshift.
+
+    Returns
+    --------
+    circular_velocity: (nm,) array_like
+        Halo circular velocity for a given mass array, cosmology and redshift, in
+        units of km s-1.
+
+    Examples
+    ---------
+    >>> import numpy as np
+    >>> from skypy.halo import properties
+
+    This example will compute the halo circular velocity, for a Planck15 cosmology at redshift 0.
+
+    >>> from astropy.cosmology import Planck15
+    >>> cosmology = Planck15
+    >>> halo_virial_mass = 10**np.arange(9.0, 12.0, 2) * Unit.Msun
+    >>> Delta_v = np.arange(201.0, 201.1, 0.1)
+    >>> redshift = np.arange(0.3, 1, 0.5)
+    >>> properties.halo_circular_velocity(halo_virial_mass, Delta_v, redshift, cosmology)
+    <Quantity [ 6.11303684, 36.72661831] km2 / (Mpc2 s2)>
+
+    References
+    ----------
+    .. [1] Barnes and Haehnelt 2010 equation 3 https://arxiv.org/pdf/1403.1873.pdf
+    """
+
+    virial_velocity = 96.6 * Unit('km s-1') * \
+        np.power(Delta_v * cosmology.Om0 * cosmology.h**2 / 24.4, 1.0/6.0) * \
+        np.sqrt((1 + redshift) / 3.3) * \
+        np.power(halo_virial_mass / (1.0e11 * Unit('Msun')), 1.0/3.0)
+
+    return virial_velocity

skypy/halos/__init__.py

@@ -0,0 +1,22 @@
+"""Halos module.
+
+This module contains methods that model the properties of dark matter halo
+populations.
+
+Models
+======
+.. autosummary::
+   :nosignatures:
+   :toctree: ../api/
+
+   colossus_mf
+"""
+
+__all__ = [
+    'colossus_mf',
+]
+
+from . import abundance_matching  # noqa F401,F403
+from . import mass  # noqa F401,F403
+from . import quenching  # noqa F401,F403
+from ._colossus import colossus_mf  # noqa F401,F403