Homogenise linear collapse threshold (#1101)

* implemented, documented, and tested

* missing deltas

include/ccl_massfunc.h

@@ -38,15 +38,6 @@ double ccl_sigmaM(ccl_cosmology *cosmo, double log_halomass, double a, int *stat
  */
 double ccl_dlnsigM_dlogM(ccl_cosmology *cosmo, double log_halomass, double a, int *status);
 
-/**
- * Fitting function for the spherical-model critical linear density for collapse
- * Fitting formula from Nakamura & Suto (1997; arXiv:astro-ph/9710107)
- * @param cosmo Cosmological parameters
- * @param a, scale factor, normalized to a=1 today
- * @param status Status flag. 0 if there are no errors, nonzero otherwise.
- */
-double dc_NakamuraSuto(ccl_cosmology *cosmo, double a, int *status);
-
 CCL_END_DECLS
 
 #endif

pyccl/halos/concentration/__init__.py

@@ -1,4 +1,4 @@
-from ..halo_model_base import Concentration
+from ..halo_model_base import Concentration, get_delta_c
 from .bhattacharya13 import *
 from .constant import *
 from .diemer15 import *

pyccl/halos/concentration/bhattacharya13.py

@@ -1,7 +1,6 @@
 __all__ = ("ConcentrationBhattacharya13",)
 
-from ... import lib
-from . import Concentration
+from . import Concentration, get_delta_c
 
 
 class ConcentrationBhattacharya13(Concentration):
@@ -32,8 +31,7 @@ def _setup(self):
 
     def _concentration(self, cosmo, M, a):
         gz = cosmo.growth_factor(a)
-        status = 0
-        delta_c, status = lib.dc_NakamuraSuto(cosmo.cosmo, a, status)
+        delta_c = get_delta_c(cosmo, a, kind='NakamuraSuto97')
         sig = cosmo.sigmaM(M, a)
         nu = delta_c / sig
         return self.A * gz**self.B * nu**self.C

pyccl/halos/concentration/diemer15.py

@@ -2,7 +2,7 @@
 
 import numpy as np
 
-from . import Concentration
+from . import Concentration, get_delta_c
 
 
 class ConcentrationDiemer15(Concentration):
@@ -42,7 +42,7 @@ def _concentration(self, cosmo, M, a):
         n = pk(k_R, a, derivative=True)
 
         sig = cosmo.sigmaM(M, a)
-        delta_c = 1.68647
+        delta_c = get_delta_c(cosmo, a, kind='EdS')
         nu = delta_c / sig
 
         floor = self.phi_0 + n * self.phi_1

pyccl/halos/concentration/ishiyama21.py

@@ -5,7 +5,7 @@
 
 from ... import lib
 from ... import check
-from . import Concentration
+from . import Concentration, get_delta_c
 
 
 class ConcentrationIshiyama21(Concentration):
@@ -85,7 +85,7 @@ def _G_inv(self, arg, n_eff):
         return np.asarray(roots)
 
     def _concentration(self, cosmo, M, a):
-        nu = 1.686 / cosmo.sigmaM(M, a)
+        nu = get_delta_c(cosmo, a, 'EdS_approx') / cosmo.sigmaM(M, a)
         n_eff = -2 * self._dlsigmaR(cosmo, M, a) - 3
         alpha_eff = cosmo.growth_rate(a)
 

pyccl/halos/halo_model_base.py

@@ -1,4 +1,5 @@
-__all__ = ("HMIngredients", "Concentration", "MassFunc", "HaloBias", )
+__all__ = ("HMIngredients", "Concentration", "MassFunc", "HaloBias",
+           "get_delta_c", )
 
 from abc import abstractmethod
 
@@ -8,6 +9,43 @@
 from .. import physical_constants as const
 
 
+def get_delta_c(cosmo, a, kind='EdS'):
+    """Returns the linear collapse threshold.
+
+    Args:
+        cosmo (:class:`~pyccl.cosmology.Cosmology`): A Cosmology object.
+        a (:obj:`float` or `array`): scale factor.
+        kind (:obj:`str`): prescription to use. Should be one of
+
+            * 'EdS': the SC prediction in Einstein de-Sitter, :math:`\\delta_c=(3/20)(12\\pi)^{2/3}`.
+            * 'EdS_approx': a common approximation to the EdS result :math:`\\delta_c=1.686`.
+            * 'NakamuraSuto97': the prescription from `Nakamura & Suto 1997 <https://arxiv.org/abs/astro-ph/9612074>`_.
+            * 'Mead16': the prescription from `Mead et al. 2016 <https://arxiv.org/abs/1602.02154>`_.
+
+    Returns:
+        (:obj:`float` or `array`): linear collapse threshold.
+    """  # noqa
+    # This is the linear collapse threshold in Einstein de-Sitter:
+    # delta_c = 3/20*(12*pi)^(2/3)
+    dc0 = 1.68647019984
+
+    if kind == 'EdS':
+        return dc0
+    elif kind == 'EdS_approx':
+        return 1.686
+    elif kind == 'NakamuraSuto97':
+        Om = cosmo.omega_x(a, 'matter')
+        return dc0*(1+0.012299*np.log10(Om))
+    elif kind == 'Mead16':
+        Om = cosmo.omega_x(a, 'matter')
+        s8 = cosmo.sigma8()*cosmo.growth_factor(a)
+        facs8 = (1.59+0.0314*np.log(s8))
+        facOm = (1+0.0123*np.log10(Om))
+        return facs8*facOm
+    else:
+        raise ValueError(f"Unknown threshold kind {kind}")
+
+
 class HMIngredients(CCLAutoRepr, CCLNamedClass):
     """Base class for halo model ingredients."""
     __repr_attrs__ = __eq_attrs__ = ("mass_def", "mass_def_strict",)

pyccl/halos/hbias/__init__.py

@@ -1,4 +1,4 @@
-from ..halo_model_base import HaloBias
+from ..halo_model_base import HaloBias, get_delta_c
 from .bhattacharya11 import *
 from .sheth01 import *
 from .sheth99 import *

pyccl/halos/hbias/bhattacharya11.py

@@ -1,6 +1,6 @@
 __all__ = ("HaloBiasBhattacharya11",)
 
-from . import HaloBias
+from . import HaloBias, get_delta_c
 
 
 class HaloBiasBhattacharya11(HaloBias):
@@ -29,7 +29,7 @@ def _setup(self):
         self.az = 0.01
         self.p = 0.807
         self.q = 1.795
-        self.dc = 1.68647
+        self.dc = get_delta_c(None, None, kind='EdS')
 
     def _get_bsigma(self, cosmo, sigM, a):
         nu = self.dc / sigM

pyccl/halos/hbias/sheth01.py

@@ -1,6 +1,6 @@
 __all__ = ("HaloBiasSheth01",)
 
-from . import HaloBias
+from . import HaloBias, get_delta_c
 
 
 class HaloBiasSheth01(HaloBias):
@@ -29,7 +29,7 @@ def _setup(self):
         self.sqrta = 0.84083292038
         self.b = 0.5
         self.c = 0.6
-        self.dc = 1.68647
+        self.dc = get_delta_c(None, None, kind='EdS')
 
     def _get_bsigma(self, cosmo, sigM, a):
         nu = self.dc/sigM

pyccl/halos/hbias/sheth99.py

@@ -1,7 +1,6 @@
 __all__ = ("HaloBiasSheth99",)
 
-from ... import check, lib
-from . import HaloBias
+from . import HaloBias, get_delta_c
 
 
 class HaloBiasSheth99(HaloBias):
@@ -40,11 +39,9 @@ def _setup(self):
 
     def _get_bsigma(self, cosmo, sigM, a):
         if self.use_delta_c_fit:
-            status = 0
-            delta_c, status = lib.dc_NakamuraSuto(cosmo.cosmo, a, status)
-            check(status, cosmo=cosmo)
+            delta_c = get_delta_c(cosmo, a, kind='NakamuraSuto97')
         else:
-            delta_c = 1.68647
+            delta_c = get_delta_c(cosmo, a, kind='EdS')
 
         nu = delta_c / sigM
         anu2 = self.a * nu**2

pyccl/halos/hbias/tinker10.py

@@ -2,7 +2,7 @@
 
 import numpy as np
 
-from . import HaloBias
+from . import HaloBias, get_delta_c
 
 
 class HaloBiasTinker10(HaloBias):
@@ -33,7 +33,7 @@ def _setup(self):
         self.B = 0.183
         self.b = 1.5
         self.c = 2.4
-        self.dc = 1.68647
+        self.dc = get_delta_c(None, None, kind='EdS')
 
     def _get_bsigma(self, cosmo, sigM, a):
         nu = self.dc / sigM

pyccl/halos/hmfunc/__init__.py

@@ -1,4 +1,4 @@
-from ..halo_model_base import MassFunc
+from ..halo_model_base import MassFunc, get_delta_c
 from .angulo12 import *
 from .bocquet16 import *
 from .despali16 import *

pyccl/halos/hmfunc/despali16.py

@@ -2,8 +2,7 @@
 
 import numpy as np
 
-from ... import check, lib
-from . import MassFunc
+from . import MassFunc, get_delta_c
 
 
 class MassFuncDespali16(MassFunc):
@@ -46,9 +45,7 @@ def _setup(self):
         self.poly_A, self.poly_a, self.poly_p = map(np.poly1d, coeffs)
 
     def _get_fsigma(self, cosmo, sigM, a, lnM):
-        status = 0
-        delta_c, status = lib.dc_NakamuraSuto(cosmo.cosmo, a, status)
-        check(status, cosmo=cosmo)
+        delta_c = get_delta_c(cosmo, a, 'NakamuraSuto97')
 
         Dv = self.mass_def.get_Delta_vir(cosmo, a)
         x = np.log10(self.mass_def.get_Delta(cosmo, a) / Dv)

pyccl/halos/hmfunc/press74.py

@@ -2,7 +2,7 @@
 
 import numpy as np
 
-from . import MassFunc
+from . import MassFunc, get_delta_c
 
 
 class MassFuncPress74(MassFunc):
@@ -28,6 +28,6 @@ def _check_mass_def_strict(self, mass_def):
         return mass_def.Delta != "fof"
 
     def _get_fsigma(self, cosmo, sigM, a, lnM):
-        delta_c = 1.68647
+        delta_c = get_delta_c(cosmo, a, kind='EdS')
         nu = delta_c/sigM
         return self._norm * nu * np.exp(-0.5 * nu**2)

pyccl/halos/hmfunc/sheth99.py

@@ -2,8 +2,7 @@
 
 import numpy as np
 
-from ... import check, lib
-from . import MassFunc
+from . import MassFunc, get_delta_c
 
 
 class MassFuncSheth99(MassFunc):
@@ -43,11 +42,9 @@ def _setup(self):
 
     def _get_fsigma(self, cosmo, sigM, a, lnM):
         if self.use_delta_c_fit:
-            status = 0
-            delta_c, status = lib.dc_NakamuraSuto(cosmo.cosmo, a, status)
-            check(status, cosmo=cosmo)
+            delta_c = get_delta_c(cosmo, a, 'NakamuraSuto97')
         else:
-            delta_c = 1.68647
+            delta_c = get_delta_c(cosmo, a, 'EdS')
 
         nu = delta_c / sigM
         return nu * self.A * (1. + (self.a * nu**2)**(-self.p)) * (

pyccl/halos/hmfunc/tinker10.py

@@ -3,7 +3,7 @@
 import numpy as np
 from scipy.interpolate import interp1d
 
-from . import MassFunc
+from . import MassFunc, get_delta_c
 
 
 class MassFuncTinker10(MassFunc):
@@ -69,7 +69,7 @@ def _setup(self):
 
     def _get_fsigma(self, cosmo, sigM, a, lnM):
         ld = np.log10(self.mass_def._get_Delta_m(cosmo, a))
-        nu = 1.686 / sigM
+        nu = get_delta_c(cosmo, a, 'EdS_approx') / sigM
         # redshift evolution only up to z=3
         a = np.clip(a, 0.25, 1)
         pa = self.pa0(ld) * a**(-0.27)

pyccl/halos/profiles/einasto.py

@@ -3,7 +3,7 @@
 import numpy as np
 from scipy.special import gamma, gammainc
 
-from .. import MassDef, mass_translator
+from .. import MassDef, mass_translator, get_delta_c
 from . import HaloProfileMatter
 
 
@@ -71,7 +71,7 @@ def _get_alpha(self, cosmo, M, a):
         if self.alpha == 'cosmo':
             Mvir = self._to_virial_mass(cosmo, M, a)
             sM = cosmo.sigmaM(Mvir, a)
-            nu = 1.686 / sM
+            nu = get_delta_c(cosmo, a, kind='EdS_approx') / sM
             return 0.155 + 0.0095 * nu * nu
         return np.full_like(M, self.alpha)
 

pyccl/tests/test_massfunction.py

@@ -30,3 +30,18 @@ def test_sigmaM_smoke(m):
     s = ccl.sigmaM(COSMO, m, a)
     assert np.all(np.isfinite(s))
     assert np.shape(s) == np.shape(m)
+
+
+def test_deltac():
+    cosmo = ccl.Cosmology(Omega_c=0.25, Omega_b=0.05,
+                          h=0.7, n_s=0.96, sigma8=0.8,
+                          transfer_function='bbks')
+    # Test EdS value
+    dca = 3*(12*np.pi)**(2/3)/20
+    dcb = ccl.halos.get_delta_c(cosmo, 1.0, kind='EdS')
+    assert np.fabs(dcb/dca-1) < 1E-5
+
+    # Test Mead et al. value
+    dca = (1.59+0.0314*np.log(0.8))*(1+0.0123*np.log10(0.3))
+    dcb = ccl.halos.get_delta_c(cosmo, 1.0, kind='Mead16')
+    assert np.fabs(dcb/dca-1) < 1E-5

src/ccl_massfunc.c

@@ -10,23 +10,6 @@
 
 #include "ccl.h"
 
-
-/*----- ROUTINE: dc_NakamuraSuto -----
-INPUT: cosmology, scale factor
-TASK: Computes the peak threshold: delta_c(z) assuming LCDM.
-Cosmology dependence of the critical linear density according to the spherical-collapse model.
-Fitting function from Nakamura & Suto (1997; arXiv:astro-ph/9710107).
-*/
-double dc_NakamuraSuto(ccl_cosmology *cosmo, double a, int *status){
-
-  double Om_mz = ccl_omega_x(cosmo, a, ccl_species_m_label, status);
-  double dc0 = (3./20.)*pow(12.*M_PI,2./3.);
-  double dc = dc0*(1.+0.012299*log10(Om_mz));
-
-  return dc;
-
-}
-
 static double sigmaM_m2r(ccl_cosmology *cosmo, double halomass, int *status)
 {
   double rho_m, smooth_radius;