Analytic Einasto and Hernquist normalization (#956)

* Add analytic _norm to Einasto and Hernquist profiles

CMakeLists.txt

@@ -22,7 +22,6 @@ set(CCL_SRC
     src/ccl_neutrinos.c
     src/ccl_emu17.c src/ccl_correlation.c
     src/ccl_halofit.c
-    src/ccl_haloprofile.c
     src/ccl_tracers.c
     src/ccl_mass_conversion.c
     src/ccl_fftlog.c)

include/ccl.h

@@ -41,7 +41,6 @@
 #include "ccl_eh.h"
 #include "ccl_halofit.h"
 #include "ccl_musigma.h"
-#include "ccl_haloprofile.h"
 #include "ccl_mass_conversion.h"
 
 CCL_BEGIN_DECLS

pyccl/ccl.i

@@ -51,7 +51,6 @@ from .errors import CCLError
 %include "ccl_covs.i"
 %include "ccl_neutrinos.i"
 %include "ccl_musigma.i"
-%include "ccl_haloprofile.i"
 %include "ccl_mass_conversion.i"
 %include "ccl_sigM.i"
 %include "ccl_f1d.i"

pyccl/halos/profiles.py

@@ -1,12 +1,10 @@
-from .. import ccllib as lib
-from ..core import check
 from ..background import h_over_h0, sigma_critical
 from ..power import sigmaM
 from ..pyutils import resample_array, _fftlog_transform
 from .concentration import Concentration
 from .massdef import MassDef
 import numpy as np
-from scipy.special import sici, erf
+from scipy.special import sici, erf, gamma, gammainc
 
 
 class HaloProfile(object):
@@ -873,6 +871,12 @@ def _get_alpha(self, cosmo, M, a, mdef):
         alpha = 0.155 + 0.0095 * nu * nu
         return alpha
 
+    def _norm(self, M, Rs, c, alpha):
+        # Einasto normalization from mass, radius, concentration and alpha
+        return M / (np.pi * Rs**3 * 2**(2-3/alpha) * alpha**(-1+3/alpha)
+                    * np.exp(2/alpha)
+                    * gamma(3/alpha) * gammainc(3/alpha, 2/alpha*c**alpha))
+
     def _real(self, cosmo, r, M, a, mass_def):
         r_use = np.atleast_1d(r)
         M_use = np.atleast_1d(M)
@@ -884,10 +888,7 @@ def _real(self, cosmo, r, M, a, mass_def):
 
         alpha = self._get_alpha(cosmo, M_use, a, mass_def)
 
-        status = 0
-        norm, status = lib.einasto_norm(R_s, R_M, alpha, M_use.size, status)
-        check(status, cosmo=cosmo)
-        norm = M_use / norm
+        norm = self._norm(M_use, R_s, c_M, alpha)
 
         x = r_use[None, :] / R_s[:, None]
         prof = norm[:, None] * np.exp(-2. * (x**alpha[:, None] - 1) /
@@ -945,6 +946,10 @@ def __init__(self, c_M_relation, truncated=True):
     def _get_cM(self, cosmo, M, a, mdef=None):
         return self.cM.get_concentration(cosmo, M, a, mdef_other=mdef)
 
+    def _norm(self, M, Rs, c):
+        # Hernquist normalization from mass, radius and concentration
+        return M / (2 * np.pi * Rs**3 * (c / (1 + c))**2)
+
     def _real(self, cosmo, r, M, a, mass_def):
         r_use = np.atleast_1d(r)
         M_use = np.atleast_1d(M)
@@ -954,10 +959,7 @@ def _real(self, cosmo, r, M, a, mass_def):
         c_M = self._get_cM(cosmo, M_use, a, mdef=mass_def)
         R_s = R_M / c_M
 
-        status = 0
-        norm, status = lib.hernquist_norm(R_s, R_M, M_use.size, status)
-        check(status, cosmo=cosmo)
-        norm = M_use / norm
+        norm = self._norm(M_use, R_s, c_M)
 
         x = r_use[None, :] / R_s[:, None]
         prof = norm[:, None] / (x * (1 + x)**3)

pyccl/tests/test_swig_interface.py

@@ -200,18 +200,3 @@ def test_swig_power():
         [0.5, 1.0],
         3,
         status)
-
-
-def test_swig_haloprofile():
-    status = 0
-    assert_raises(CCLError,
-                  ccllib.einasto_norm,
-                  [0.1, 1.0],
-                  [0.1, 1.0],
-                  [0.1, 1.0],
-                  4, status)
-    assert_raises(CCLError,
-                  ccllib.hernquist_norm,
-                  [0.1, 1.0],
-                  [0.1, 1.0],
-                  4, status)