Merge branch 'main' into galaxyq

Author: itrharrison

docs/galaxies.rst

@@ -95,6 +95,18 @@ The following models are found in the `skypy.galaxies.stellar_mass` package.
   schechter_smf_phi_satellite_quenched
 
 
+Velocity dispersion
+-------------------
+
+The following models are found in the `skypy.galaxies.velocity_dispersion` package.
+
+.. currentmodule:: skypy.galaxies.velocity_dispersion
+.. autosummary::
+  :nosignatures:
+
+  schechter_vdf
+
+
 Reference/API
 =============
 
@@ -107,3 +119,4 @@ Reference/API
 .. automodapi:: skypy.galaxies.spectrum
    :include-all-objects:
 .. automodapi:: skypy.galaxies.stellar_mass
+.. automodapi:: skypy.galaxies.velocity_dispersion

skypy/galaxies/__init__.py

@@ -13,6 +13,7 @@
 from . import redshift  # noqa F401,F403
 from . import spectrum  # noqa F401,F403
 from . import stellar_mass  # noqa F401,F403
+from . import velocity_dispersion  # noqa F401,F403
 
 from ._schechter import schechter_lf  # noqa
 from ._schechter import schechter_smf  # noqa

skypy/galaxies/tests/test_velocity_dispersion.py

@@ -0,0 +1,38 @@
+import numpy as np
+from scipy.stats import kstest
+from scipy.special import gamma
+from scipy.integrate import cumtrapz
+
+
+def test_schechter_vdf():
+    from skypy.galaxies.velocity_dispersion import schechter_vdf
+
+    # test schechter velocity dispersion function
+    vd_min = 1
+    vd_max = 300
+    size = 5000
+    alpha = 2.32
+    beta = 2.67
+    vd_star = 161
+    resolution = 3000
+    phi_star = 8e-3
+    samples = schechter_vdf(alpha=alpha, beta=beta, vd_star=vd_star,
+                            vd_min=vd_min, vd_max=vd_max, size=size, resolution=resolution)
+
+    # test output is within limits and size
+    assert len(samples) == size
+    assert vd_min <= samples.all() <= vd_max
+
+    # test that sampling corresponds to sampling from the right pdf
+    def calc_pdf(vd):
+        return phi_star*(vd/vd_star)**alpha*np.exp(-(vd/vd_star)**beta) * \
+            (beta/gamma(alpha/beta))*(1/vd)
+
+    def calc_cdf(m):
+        pdf = calc_pdf(m)
+        cdf = cumtrapz(pdf, m, initial=0)
+        cdf /= cdf[-1]
+        return cdf
+
+    p_value = kstest(samples, calc_cdf)[1]
+    assert p_value > 0.01

skypy/galaxies/velocity_dispersion.py

@@ -0,0 +1,68 @@
+r"""Models of galaxy velocity dispersion.
+
+"""
+
+import numpy as np
+from skypy.utils.random import schechter
+
+__all__ = [
+    'schechter_vdf',
+]
+
+
+def schechter_vdf(alpha, beta, vd_star, vd_min, vd_max, size=None, resolution=1000):
+    r"""Sample velocity dispersion of elliptical galaxies in the local universe
+    following a Schecter function.
+
+    Parameters
+    ----------
+    alpha: int
+        The alpha parameter in the modified Schechter equation.
+    beta: int
+        The beta parameter in the modified Schechter equation.
+    vd_star: int
+        The characteristic velocity dispersion.
+    vd_min, vd_max: int
+        Lower and upper bounds of random variable x. Samples are drawn uniformly from bounds.
+    resolution: int
+        Resolution of the inverse transform sampling spline. Default is 100.
+    size: int
+        Number of samples returned. Default is 1.
+
+    Returns
+    -------
+    velocity_dispersion: array_like
+        Velocity dispersion drawn from Schechter function.
+
+    Notes
+    -----
+    The probability distribution function :math:`p(\sigma)` for velocity dispersion :math:`\sigma`
+    can be described by a Schechter function (see eq. (4) in [1]_)
+
+    .. math::
+
+        \phi = \phi_* \left(\frac{\sigma}{\sigma_*}\right)^\alpha
+            \exp\left[-\left( \frac{\sigma}{\sigma_*} \right)^\beta\right]
+            \frac{\beta}{\Gamma(\alpha/\beta)} \frac{1}{\sigma} \mathrm{d}\sigma \;.
+
+    where :math:`\Gamma` is the gamma function, :math:`\sigma_*` is the
+    characteristic velocity dispersion, :math:`\phi_*` is
+    number density of all spiral galaxies and
+    :math:`\alpha` and :math:`\beta` are free parameters.
+
+    References
+    ----------
+    .. [1] Choi, Park and Vogeley, (2007), astro-ph/0611607, doi:10.1086/511060
+
+    """
+
+    if np.ndim(alpha) > 0:
+        raise NotImplementedError('only scalar alpha is supported')
+
+    alpha_prime = alpha/beta - 1
+    x_min, x_max = (vd_min/vd_star)**beta, (vd_max/vd_star)**beta
+
+    samples = schechter(alpha_prime, x_min, x_max, resolution=resolution, size=size)
+    samples = samples**(1/beta) * vd_star
+
+    return samples