From d60129ecd23e4a0a7e2aa138c412bafd0f283724 Mon Sep 17 00:00:00 2001 From: Sut Ieng Tam Date: Fri, 28 May 2021 18:32:12 +0800 Subject: [PATCH] Add cluster features --- skypy/halos/_colossus.py | 134 +++++++++++++++++++++++++++++ skypy/halos/tests/test_colossus.py | 52 ++++++++++- 2 files changed, 185 insertions(+), 1 deletion(-) diff --git a/skypy/halos/_colossus.py b/skypy/halos/_colossus.py index 2ecd048d7..1337e2fad 100644 --- a/skypy/halos/_colossus.py +++ b/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(= 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)