"""
This module calculates the corrections to the standard lensed CMB power spectra results due to time delay and
emission angle, following `arXiv:1706.02673 <https://arxiv.org/abs/1706.02673>`_. This can be combined with the result
from the postborn module to estimate the leading corrections to the standard lensing B modes.
Corrections to T and E are negligible, and not calculated. The result for BB includes approximately contributions
from reionization, but this can optionally be turned off.
"""
from . import camb, model
import numpy as np
from .mathutils import threej
def cl_deflection_limber(results, PK, ls, nz, chi_source, emit_pow=2, lens_pow=0):
chi_source = np.float64(chi_source)
chis = np.linspace(0, chi_source, nz, dtype=np.float64)
zs = results.redshift_at_comoving_radial_distance(chis)
dchis = (chis[2:] - chis[:-2]) / 2
chis = chis[1:-1]
zs = zs[1:-1]
win = 1 / (chis * chi_source) ** emit_pow
if lens_pow:
win *= (-(1 - chis / chi_source) / chis ** 2) ** lens_pow
cl = np.zeros(ls.shape)
w = np.ones(chis.shape)
for i, l in enumerate(ls):
k = (l + 0.5) / chis
w[:] = 1
w[k < 1e-4] = 0
w[k >= PK.kmax] = 0
cl[i] = np.dot(dchis, w * PK.P(zs, k, grid=False) * win / k ** 4)
cl *= 4 * (ls * (ls + 1))
return cl
[docs]
def get_emission_angle_powers(camb_background, PK, chi_source, lmax=3000, acc=1, lsamp=None):
r"""
Get the power spectrum of :math:`\psi_d`, the potential for the emission angle, and its cross with standard lensing.
Uses the Limber approximation (and assumes flat universe).
:param camb_background: a CAMB results object, used for calling background functions
:param PK: a matter power spectrum interpolator (from camb.get_matter_power_interpolator)
:param chi_source: comoving radial distance of source in Mpc
:param lmax: maximum L
:param acc: accuracy parameter
:param lsamp: L sampling for the result
:return: a InterpolatedUnivariateSpline object containing :math:`L(L+1) C_L`
"""
from scipy.interpolate import InterpolatedUnivariateSpline
assert (isinstance(camb_background, camb.CAMBdata) and np.isclose(camb_background.Params.omk, 0))
nz = int(100 * acc)
ls = lsamp or np.hstack((np.arange(2, 60, 2), np.arange(60, min(lmax, 400), 10),
np.arange(min(lmax, 400), lmax, int(50. / acc)),
np.arange(lmax, lmax + 1))).astype(np.float64)
cl_psi_d = cl_deflection_limber(camb_background, PK, ls, nz, chi_source, emit_pow=2, lens_pow=0)
cl_psi_d_cross = cl_deflection_limber(camb_background, PK, ls, nz, chi_source, emit_pow=1, lens_pow=1)
return InterpolatedUnivariateSpline(ls, cl_psi_d), InterpolatedUnivariateSpline(ls, cl_psi_d_cross)
[docs]
def get_emission_delay_BB(params, kmax=100, lmax=3000, non_linear=True, CMB_unit='muK', raw_cl=False,
acc=1, lsamp=None, return_terms=False, include_reionization=True):
r"""
Get B modes from emission angle and time delay effects.
Uses full-sky result from appendix of `arXiv:1706.02673 <https://arxiv.org/abs/1706.02673>`_
:param params: :class:`.model.CAMBparams` instance with cosmological parameters etc.
:param kmax: maximum k (in :math:`{\rm Mpc}^{-1}` units)
:param lmax: maximum :math:`\ell`
:param non_linear: include non-linear corrections
:param CMB_unit: normalization for the result
:param raw_cl: if true return :math:`C_\ell`, else :math:`\ell(\ell+1)C_\ell/2\pi`
:param acc: accuracy setting, increase to test stability
:param lsamp: array of :math:`\ell` values to compute output at. If not set, set to sampling good for interpolation
:param return_terms: return the three sub-terms separately rather than the total
:param include_reionization: approximately include reionization terms by second scattering surface
:return: InterpolatedUnivariateSpline for :math:`C_\ell^{BB}`
"""
from scipy.interpolate import InterpolatedUnivariateSpline
assert (np.isclose(params.omk, 0))
camb_background = camb.get_background(params)
chi_source = camb_background.tau0 - camb_background.tau_maxvis
z_source = camb_background.redshift_at_comoving_radial_distance(chi_source)
PK = camb.get_matter_power_interpolator(params, nonlinear=non_linear,
hubble_units=False, k_hunit=False, kmax=kmax,
var1=model.Transfer_Weyl, var2=model.Transfer_Weyl, zmax=z_source)
assert (lmax > 250)
lsampvelcl = np.hstack(
(np.arange(2, 20, 2), np.arange(25, 200, 20 // acc), np.arange(220, lmax, 40 // acc)))
lmax_e = max(1500, lmax * 2)
pars = params.copy()
pars.set_for_lmax(lmax_e, lens_potential_accuracy=1)
cmb = get_source_cmb_cl(pars, CMB_unit=CMB_unit)
totautoB = np.zeros(lsampvelcl.shape)
totBEterm = np.zeros(lsampvelcl.shape)
totBxterm = np.zeros(lsampvelcl.shape)
for reion in [False, True]:
if reion:
if not include_reionization:
break
zreion = params.get_zre()
chi_source = camb_background.tau0 - camb_background.conformal_time(zreion)
lmax_e = 300
tag_E = 'E2'
tag_zeta = 'emit2'
lstep = 1
else:
tag_E = 'E1'
tag_zeta = 'emit1'
lstep = 5
cl_psi_d_sp, cl_psi_d_x_lens_sp = get_emission_angle_powers(camb_background, PK, chi_source, lmax_e, acc, lsamp)
lsarr = np.arange(2, lmax_e + 1, dtype=np.float64)
llp1 = lsarr * (lsarr + 1.)
cdd = cl_psi_d_sp(lsarr) * (lsarr + 2) * (lsarr - 1) * (2 * lsarr + 1)
cd = cl_psi_d_sp(lsarr) / llp1 * (2 * lsarr + 1)
cxdd = cl_psi_d_x_lens_sp(lsarr) * (2 * lsarr + 1)
cxd = cxdd / llp1
cEE = cmb['%sx%s' % (tag_E, tag_E)][2:lmax_e + 1] / llp1 * (2 * lsarr + 1) # raw CL
cEx = cmb['%sx%s' % (tag_E, tag_zeta)][2:lmax_e + 1] * (2 * lsarr + 1)
cExx = cEx / llp1
czeta = cmb['%sx%s' % (tag_zeta, tag_zeta)][2:lmax_e + 1] / llp1 * (2 * lsarr + 1)
for i, ll in enumerate(lsampvelcl):
if reion and ll > lmax_e:
break
for llp in range(2, lmax_e, lstep):
lp = np.float64(llp)
wig = threej(llp, ll, 2, -2)
minl = np.abs(llp - ll)
if minl < 2:
wig = wig[2 - minl:]
wigx = threej(llp, ll, 0, -2)
minl = max(2, np.abs(llp - ll))
maxl = min(lmax_e, np.abs(llp + ll))
off = 0
if (minl + llp + ll) % 2 == 0:
off = 1
wig2 = wig[off:maxl - minl + 1:2] ** 2
totautoB[i] += lstep * np.dot(wig2, czeta[minl + off - 2:maxl + 1 - 2:2]) * cdd[llp - 2]
totBEterm[i] += lstep * np.dot(wig2,
cd[minl + off - 2:maxl + 1 - 2:2] - cxdd[minl + off - 2:maxl + 1 - 2:2]
- (lp * (lp + 1) - ll * (ll + 1)) * cxd[minl + off - 2:maxl + 1 - 2:2]
) * cEE[llp - 2]
wigx2 = wigx[off:maxl - minl + 1:2] * wig[off:maxl - minl + 1:2]
totBxterm[i] += lstep * (np.dot(wigx2, cExx[minl + off - 2:maxl + 1 - 2:2])
* (2 * cd[llp - 2] - ((lp * (lp + 1) - ll * (ll + 1)) * cxd[llp - 2]))
- np.dot(wigx2, cEx[minl + off - 2:maxl + 1 - 2:2]) * cxd[llp - 2]
) * np.sqrt(lp * (lp + 1) * (lp + 2) * (lp - 1))
fac = 1 / 2. # (4 * np.pi) [CMB CL already have 1/2pi in ]
if not raw_cl:
fac *= lsampvelcl * (lsampvelcl + 1) / (2 * np.pi)
totautoB *= fac
totBEterm *= fac
totBxterm *= fac
if return_terms:
return InterpolatedUnivariateSpline(lsampvelcl, totautoB), \
InterpolatedUnivariateSpline(lsampvelcl, totBEterm), \
InterpolatedUnivariateSpline(lsampvelcl, totBxterm)
else:
return InterpolatedUnivariateSpline(lsampvelcl, totBxterm + totBEterm + totautoB)
[docs]
def get_source_cmb_cl(params, CMB_unit='muK'):
r"""
Get the angular power spectrum of emission angle and time delay sources :math:`\psi_t`, :math:`\psi_\zeta`,
as well as the perpendicular velocity and E polarization.
All are returned with 1 and 2 versions, for recombination and reionization respectively.
Note that this function destroys any custom sources currently configured.
:param params: :class:`.model.CAMBparams` instance with cosmological parameters etc.
:param CMB_unit: scale results from dimensionless, use 'muK' for :math:`\mu K^2` units
:return: dictionary of power spectra, with :math:`L(L+1)/2\pi` factors.
"""
import sympy
from sympy import diff
from . import symbolic as cs
assert (np.isclose(params.omk, 0))
angdist = cs.tau0 - cs.t
emission_sources = {
'vperp': -(cs.sigma + cs.v_b) * cs.visibility / cs.k / angdist,
'emit': 15 * diff(cs.polter * cs.visibility, cs.t) / 8 / cs.k ** 2 / angdist,
'delay': 15 * diff(cs.polter * cs.visibility, cs.t) / 8 / cs.k ** 2 / angdist ** 2 * (cs.tau0 - cs.tau_maxvis),
'E': cs.scalar_E_source}
sources = {}
for key, source in list(emission_sources.items()):
sources[key + '1'] = sympy.Piecewise((source, 1 / cs.a - 1 > 30), (0, True)) # recombination
sources[key + '2'] = sympy.Piecewise((source, 1 / cs.a - 1 <= 30), (0, True)) # reionization
params.set_custom_scalar_sources(sources, source_ell_scales={'E1': 2, 'E2': 2})
return camb.get_results(params).get_cmb_unlensed_scalar_array_dict(CMB_unit=CMB_unit)