ENH: Add concentration, radius, Delta_Sigma functions to halos module

skypy/halos/_colossus.py

@@ -203,3 +203,137 @@ def colossus_mf(redshift, model, mdef, m_min, m_max, sky_area, cosmology,
                                         cosmology, sigma8, ns, size, resolution)
 
     return z, m
+
+
+def concentration(mass, mdef, redshift, model, cosmology, sigma8, ns):
+    r'''Halo concentration calculator.
+
+    This function calculates halo concentration(s) using the model of c-M relation
+    available in colossus.
+
+    Parameters
+    ----------
+    mass : float or array_like
+        Spherical overdensity halo mass in units of solar mass/h, corresponding
+        to the mass definition, mdef.
+    mdef : str
+        Halo mass definition for spherical overdensities used by colossus.
+    redshift : float
+        Halo redshift
+    model : string
+        The model of the c-M relation which is available in colossus.
+    cosmology : astropy.cosmology.Cosmology
+        Astropy cosmology object
+    sigma8 : float
+        Cosmology parameter, amplitude of the (linear) power spectrum on the
+        scale of 8 Mpc/h.
+    ns : float
+        Cosmology parameter, spectral index of scalar perturbation power spectrum.
+
+    Returns
+    -------
+    concentration : float or array_like
+        Halo concentration(s); has the same dimensions as mass.
+
+    '''
+    from colossus.cosmology.cosmology import fromAstropy
+    from colossus.halo import concentration
+
+    fromAstropy(cosmology, sigma8=sigma8, ns=ns)
+
+    c = concentration.concentration(mass, mdef, redshift, model)
+
+    return c
+
+
+def radius(mass, concentration, redshift, mdef, Delta, cosmology, sigma8, ns):
+    r'''Calculate the scale radius and the spherical overdensity radius of halo by assuming
+    the NFW model.
+
+    This function calculates the scale radius and any spherical overdensity radius
+    for NFW dark matter halo.
+
+    Parameters
+    ----------
+    mass : float
+        A spherical overdensity halo mass in units of solar mass/h, corresponding
+        to the mass definition, mdef.
+    concentration : float
+        The concentration corresponding to the given halo mass and mass definition.
+    redshift : float
+        The halo redshift value.
+    mdef : str
+        Halo mass definition for spherical overdensities used by colossus.
+    Delta : str
+        The mass definition for which the spherical overdensity radius is computed.
+    cosmology : astropy.cosmology.Cosmology
+        Astropy cosmology object
+    sigma8 : float
+        Cosmology parameter, amplitude of the (linear) power spectrum on the
+        scale of 8 Mpc/h.
+    ns : float
+        Cosmology parameter, spectral index of scalar perturbation power spectrum.
+
+    Returns
+    -------
+    rs : float
+        The scale radius in physical kpc/h.
+    RDelta : float
+        Spherical overdensity radius of a given mass definition Delta.
+
+    '''
+    from colossus.cosmology.cosmology import fromAstropy
+    from colossus.halo import profile_nfw
+
+    fromAstropy(cosmology, sigma8=sigma8, ns=ns)
+
+    prof = profile_nfw.NFWProfile(M=mass, c=concentration, z=redshift, mdef=mdef)
+    rs = prof.par['rs']
+    RDelta = prof.RDelta(redshift, Delta)
+
+    return rs, RDelta
+
+
+def Delta_Sigma(mass, concentration, redshift, mdef, radius, cosmology, sigma8, ns):
+    r'''The excess surface density at given radius by assuming the NFW model.
+
+    This function uses Colossus routines to compute the excess surface density profile,
+    which is defined as Delta_Sigma(R) = Sigma(<R) − Sigma(R).
+
+    Parameters
+    ----------
+    mass : float
+        A spherical overdensity halo mass in units of solar mass/h, corresponding
+        to the mass definition, mdef.
+    concentration : float
+        The concentration corresponding to the given halo mass and mass definition.
+    redshift : float
+        Halo redshift
+    mdef : str
+        Halo mass definition for spherical overdensities used by colossus.
+    radius : float or array_like
+        Radius in physical kpc/h.
+    cosmology : astropy.cosmology.Cosmology
+        Astropy cosmology object
+    sigma8 : float
+        Cosmology parameter, amplitude of the (linear) power spectrum on the
+        scale of 8 Mpc/h.
+    ns : float
+        Cosmology parameter, spectral index of scalar perturbation power spectrum.
+
+    Returns
+    -------
+    DeltaSigma: float or array_like
+        The excess surface density at the given radius, in units of h physical Msun/kpc^2;
+        has the same dimensions as radius.
+
+    '''
+    from colossus.cosmology.cosmology import fromAstropy
+    from colossus.halo import profile_nfw
+
+    fromAstropy(cosmology, sigma8=sigma8, ns=ns)
+
+    prof = profile_nfw.NFWProfile(M=mass, c=concentration, z=redshift, mdef=mdef)
+    deltaSigma = prof.deltaSigma(radius)
+
+    return deltaSigma

skypy/halos/tests/test_colossus.py

@@ -11,7 +11,7 @@ def test_colossus_mass_sampler():
     from astropy.cosmology import WMAP9
     from colossus.cosmology.cosmology import fromAstropy
     from colossus.lss import mass_function
-    from skypy.halos.mass import colossus_mass_sampler
+    from skypy.halos._colossus import colossus_mass_sampler
     m_min, m_max, size = 1e+12, 1e+16, 100
     sample = colossus_mass_sampler(redshift=0.1, model='despali16',
                                    mdef='500c', m_min=m_min, m_max=m_max,
@@ -105,3 +105,53 @@ def test_colossus_mf():
     assert np.all(z_halo <= np.max(z))
     assert np.all(m_halo >= m_min)
     assert np.all(m_halo <= m_max)
+
+
+@pytest.mark.skipif(not HAS_COLOSSUS, reason='test requires colossus')
+@pytest.mark.flaky
+def test_colossus_concentration():
+    from skypy.halos._colossus import concentration
+    from astropy.cosmology import Planck15
+
+    mass = np.logspace(np.log10(1e+14), np.log10(1e+15), 10)
+    model, mdef = 'diemer19', '200c'
+    redshift = 0.5
+    cosmo = Planck15
+    sigma8, ns = 0.82, 0.96
+    c = concentration(mass, mdef, redshift, model, cosmo, sigma8, ns)
+
+    assert len(c) == len(mass)
+
+
+@pytest.mark.skipif(not HAS_COLOSSUS, reason='test requires colossus')
+@pytest.mark.flaky
+def test_colossus_radius():
+    from skypy.halos._colossus import radius
+    from astropy.cosmology import Planck15
+
+    mass, c, redshift = 1e+15, 3, 0.5
+    Delta, mdef = '500c', '200c'
+    cosmo = Planck15
+    sigma8, ns = 0.82, 0.96
+
+    rs, RDelta = radius(mass, c, redshift, mdef, Delta, cosmo, sigma8, ns)
+
+    assert rs > 0
+    assert RDelta > 0
+
+
+@pytest.mark.skipif(not HAS_COLOSSUS, reason='test requires colossus')
+@pytest.mark.flaky
+def test_colossus_Delta_Sigma():
+    from skypy.halos._colossus import Delta_Sigma
+    from astropy.cosmology import Planck15
+
+    mass, c, redshift = 1e+15, 3, 0.5
+    mdef = '200c'
+    cosmo = Planck15
+    sigma8, ns = 0.82, 0.96
+    radius = np.logspace(np.log10(300), np.log10(3000), 10)
+
+    DeltaSigma = Delta_Sigma(mass, c, redshift, mdef, radius, cosmo, sigma8, ns)
+
+    assert len(DeltaSigma) == len(radius)