ISW (#861)

* ISW implemented

* benchmarked

* fixed docs

* whitespace

* test name

benchmarks/test_isw.py

@@ -0,0 +1,49 @@
+import numpy as np
+import pyccl as ccl
+from scipy.integrate import simps
+
+
+def test_iswcl():
+    # Cosmology
+    Ob = 0.05
+    Oc = 0.25
+    h = 0.7
+    COSMO = ccl.Cosmology(
+        Omega_b=Ob,
+        Omega_c=Oc,
+        h=h,
+        n_s=0.96,
+        sigma8=0.8,
+        transfer_function='bbks')
+
+    # CCL calculation
+    ls = np.arange(2, 100)
+    zs = np.linspace(0, 0.6, 256)
+    nz = np.exp(-0.5*((zs-0.3)/0.05)**2)
+    bz = np.ones_like(zs)
+    tr_n = ccl.NumberCountsTracer(COSMO, has_rsd=False,
+                                  dndz=(zs, nz), bias=(zs, bz))
+    tr_i = ccl.ISWTracer(COSMO)
+    cl = ccl.angular_cl(COSMO, tr_n, tr_i, ls)
+
+    # Benchmark from Eq. 6 in 1710.03238
+    pz = nz / simps(nz, x=zs)
+    H0 = h / ccl.physical_constants.CLIGHT_HMPC
+    # Prefactor
+    prefac = 3*COSMO['T_CMB']*(Oc+Ob)*H0**3/(ls+0.5)**2
+    # H(z)/H0
+    ez = ccl.h_over_h0(COSMO, 1./(1+zs))
+    # Linear growth and derivative
+    dz = ccl.growth_factor(COSMO, 1./(1+zs))
+    gz = np.gradient(dz*(1+zs), zs[1]-zs[0])/dz
+    # Comoving distance
+    chi = ccl.comoving_radial_distance(COSMO, 1/(1+zs))
+    # P(k)
+    pks = np.array([ccl.nonlin_matter_power(COSMO, (ls+0.5)/(c+1E-6), 1./(1+z))
+                    for c, z in zip(chi, zs)]).T
+    # Limber integral
+    cl_int = pks[:, :]*(pz*ez*gz)[None, :]
+    clbb = simps(cl_int, x=zs)
+    clbb *= prefac
+
+    assert np.all(np.fabs(cl/clbb-1) < 1E-3)

pyccl/__init__.py

@@ -45,7 +45,7 @@
 
 # Cl's and tracers
 from .tracers import Tracer, NumberCountsTracer, WeakLensingTracer, CMBLensingTracer, \
-    tSZTracer, get_density_kernel, get_kappa_kernel, get_lensing_kernel
+    tSZTracer, 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/tracers.py

@@ -1,7 +1,7 @@
 from . import ccllib as lib
 from .core import check
 from .background import comoving_radial_distance, growth_rate, \
-    growth_factor, scale_factor_of_chi
+    growth_factor, scale_factor_of_chi, h_over_h0
 from .pyutils import _check_array_params, NoneArr, _vectorize_fn6
 import numpy as np
 
@@ -701,6 +701,45 @@ def __init__(self, cosmo, z_max=6., n_chi=1024):
         self.add_tracer(cosmo, kernel=(chi_arr, w_arr))
 
 
+class ISWTracer(Tracer):
+    """Specific :class:`Tracer` associated with the integrated Sachs-Wolfe
+    effect (ISW). Useful when cross-correlating any low-redshift probe with
+    the primary CMB anisotropies. The ISW contribution to the temperature
+    fluctuations is:
+
+    .. math::
+        \\Delta T_{\\rm CMB} =
+        2T_{\\rm CMB} \\int_0^{\\chi_{LSS}}d\\chi a\\,\\dot{\\phi}
+
+    Any angular power spectra computed with this tracer, should use
+    a three-dimensional power spectrum involving the matter power spectrum.
+    The current implementation of this tracers assumes a standard Poisson
+    equation relating :math:`\\phi` and :math:`\\delta`, and linear structure
+    growth. Although this should be valid in :math:`\\Lambda` CDM and on
+    the large scales the ISW is sensitive to, these approximations must be
+    borne in mind.
+
+    Args:
+        cosmo (:class:`~pyccl.core.Cosmology`): Cosmology object.
+        zmax (float): maximum redshift up to which we define the
+            kernel.
+        n_chi (float): number of intervals in the radial comoving
+            distance on which we sample the kernel.
+    """
+    def __init__(self, cosmo, z_max=6., n_chi=1024):
+        self.chi_max = comoving_radial_distance(cosmo, 1./(1+z_max))
+        chi = np.linspace(0, self.chi_max, n_chi)
+        a_arr = scale_factor_of_chi(cosmo, chi)
+        H0 = cosmo['h'] / lib.cvar.constants.CLIGHT_HMPC
+        OM = cosmo['Omega_c']+cosmo['Omega_b']
+        Ez = h_over_h0(cosmo, a_arr)
+        fz = growth_rate(cosmo, a_arr)
+        w_arr = 3*cosmo['T_CMB']*H0**3*OM*Ez*chi**2*(1-fz)
+
+        self._trc = []
+        self.add_tracer(cosmo, kernel=(chi, w_arr), der_bessel=-1)
+
+
 def _check_returned_tracer(return_val):
     """Wrapper to catch exceptions when tracers are spawned from C.
     """