CIB profiles (#905)

* added Shang et al. 2012 CIB profile

* tracer added

* benchmarked

* docs

* bringing back ISW

* more tests

* more tests and bugfix

* Nick's comments

benchmarks/data/codes/cib_cl_benchmark.sh

@@ -0,0 +1,10 @@
+#!/bin/bash
+
+# This generates the CIB power spectrum benchmarks
+git clone git@github.com:borisbolliet/class_sz.git
+cd class_sz
+make clean ; make -j8 class
+cd ..
+
+./class_sz/class class_sz_cib.ini
+rm ../cib_class_sz_cib.txt ../cib_class_sz_pk.dat ../cib_class_sz_redshift_dependent_functions.txt ../cib_class_sz_tk.dat

benchmarks/data/codes/class_sz_cib.ini

@@ -0,0 +1,54 @@
+path_to_class = /home/damonge/Science/Codes/CCL/benchmarks/data/codes/class_sz
+h = 0.7
+tau_reio = 0.07
+n_s = 0.9645
+k_pivot = 0.05
+Omega_b = 0.05
+Omega_cdm = 0.25
+A_s = 2.02e-9
+N_ncdm = 0
+N_ur = 3.046
+deg_ncdm = 3
+m_ncdm = 0.0
+T_ncdm = 0.71611
+output = cib_cib_1h, cib_cib_2h
+M_min = 1e8
+M_max = 1e16
+z_min = 0.07
+z_max = 6.
+freq_min = 10.
+freq_max = 5e4
+root = ../cib_class_sz_
+write sz results to files = yes
+z_max_pk = 6.
+mass function = T10
+concentration parameter = D08
+delta for cib = 200m
+hm_consistency = 0
+damping_1h_term = 0
+M_min_HOD = 1E10
+M1_prime_HOD = 3.3E12
+Redshift evolution of dust temperature =  0.36
+Dust temperature today in Kelvins = 24.4
+Emissivity index of sed = 1.75
+Power law index of SED at high frequency = 1.7
+Redshift evolution of L − M normalisation = 3.6
+Most efficient halo mass in Msun = 3.981072E+12
+Normalisation of L − M relation in [Jy MPc2/Msun/Hz] = 6.4e-8
+Size of of halo masses sourcing CIB emission = 0.5
+cib_frequency_list_num = 1
+cib_frequency_list_in_GHz = 217
+# cib_frequency_list_num = 5
+# cib_frequency_list_in_GHz = '217,353,545,857,3000'
+Frequency_id nu for cib in GHz (to save in file) = 0
+Frequency_id nu^prime for cib in GHz (to save in file) = 0
+dlogell = 0.4
+ell_max = 5000
+ell_min = 2
+mass_epsrel = 1e-4
+mass_epsabs = 1e-40
+redshift_epsrel = 1e-4
+redshift_epsabs = 1e-40
+sz_verbose = 2
+pressure_profile_epsabs = 1e-8
+pressure_profile_epsrel = 1e-3

benchmarks/test_cib.py

@@ -0,0 +1,48 @@
+import numpy as np
+import pyccl as ccl
+
+
+def test_cibcl():
+    # Read benchmarks
+    bm = np.loadtxt("benchmarks/data/cib_class_sz_szpowerspectrum.txt",
+                    unpack=True)
+
+    # Initialize cosmology and halo model ingredients
+    cosmo = ccl.Cosmology(Omega_b=0.05,
+                          Omega_c=0.25,
+                          h=0.7,
+                          n_s=0.9645,
+                          A_s=2.02E-9,
+                          Neff=3.046)
+    mdef = ccl.halos.MassDef200m()
+    cM = ccl.halos.ConcentrationDuffy08(mdef)
+    nM = ccl.halos.MassFuncTinker10(cosmo, mdef, norm_all_z=True)
+    bM = ccl.halos.HaloBiasTinker10(cosmo, mdef)
+    hmc = ccl.halos.HMCalculator(cosmo, nM, bM, mdef)
+    pr = ccl.halos.HaloProfileCIBShang12(cM, 217, Mmin=1E10)
+    pr.update_parameters(nu_GHz=217,
+                         alpha=0.36,
+                         T0=24.4,
+                         beta=1.75,
+                         gamma=1.7,
+                         s_z=3.6,
+                         log10meff=12.6,
+                         sigLM=np.sqrt(0.5),
+                         Mmin=1E10,
+                         L0=6.4E-8)
+
+    pr2pt = ccl.halos.Profile2ptCIB()
+
+    # CIB tracer
+    tr = ccl.CIBTracer(cosmo, z_min=0.07)
+
+    # 3D power spectrum
+    pk = ccl.halos.halomod_Pk2D(cosmo, hmc, pr, prof_2pt=pr2pt)
+
+    # Angular power spectrum
+    ls = bm[0]
+    cl = ccl.angular_cl(cosmo, tr, tr, ls, p_of_k_a=pk)
+    dl = cl*ls*(ls+1)/(2*np.pi)
+
+    # Compare
+    assert np.all(np.fabs(dl/(bm[31]+bm[32])-1) < 1E-2)

pyccl/__init__.py

@@ -45,7 +45,7 @@
 
 # Cl's and tracers
 from .tracers import Tracer, NumberCountsTracer, WeakLensingTracer, CMBLensingTracer, \
-    tSZTracer, ISWTracer, get_density_kernel, get_kappa_kernel, get_lensing_kernel
+    tSZTracer, CIBTracer, ISWTracer, get_density_kernel, get_kappa_kernel, get_lensing_kernel
 from .cls import angular_cl
 from .covariances import angular_cl_cov_cNG, angular_cl_cov_SSC, sigma2_B_disc, sigma2_B_from_mask
 

pyccl/cls.py

@@ -34,7 +34,7 @@ def angular_cl(cosmo, cltracer1, cltracer2, ell, p_of_k_a=None,
         limber_integration_method (string) : integration method to be used
             for the Limber integrals. Possibilities: 'qag_quad' (GSL's `qag`
             method backed up by `quad` when it fails) and 'spline' (the
-            integrand is splined and then integrated analytically).
+            integrand is splined and then integrated numerically).
 
     Returns:
         float or array_like: Angular (cross-)power spectrum values, \

pyccl/covariances.py

@@ -132,12 +132,12 @@ def angular_cl_cov_cNG(cosmo, cltracer1, cltracer2, ell, tkka, fsky=1.,
 
 
 def sigma2_B_disc(cosmo, a=None, fsky=1., p_of_k_a=None):
-    """ Returns the variance of the projected linear density field
+    """Returns the variance of the projected linear density field
     over a circular disc covering a sky fraction `fsky` as a function
     of scale factor. This is given by
 
     .. math::
-        \\sigma^2_B(z) = \\int_0^\\infty \frac{k\\,dk}{2\\pi}
+        \\sigma^2_B(z) = \\int_0^\\infty \\frac{k\\,dk}{2\\pi}
             P_L(k,z)\\,\\left[\\frac{2J_1(k R(z))}{k R(z)}\\right]^2,
 
     where :math:`R(z)` is the corresponding radial aperture as a

pyccl/halos/__init__.py

@@ -52,3 +52,8 @@
     halomod_trispectrum_1h,
     halomod_Tk3D_1h,
     halomod_Tk3D_SSC)
+
+# CIB profiles
+from .profiles_cib import (  # noqa
+    HaloProfileCIBShang12,
+    Profile2ptCIB)