In [1]:
# Notebook demonstrating the use of quantities defined incamb.symbolic, and examples of usage
# for defining custom sources and plotting different quantities.

# Set of scalar equations implemented in CAMB, and calculation of the line-of-sight sources
# indices are:
# g - photons, r- massless neutrinos, c - CDM, b - baryons, de - dark energy, nu - massive neutrinos

# kappa = 8 pi G, Pi is anisotropic stress, q =(rho+p)v the heat flux
# (and for components rho_i q_i = (rho_i+p_i)v_i)
# z is the perturbation to the expansion, h perturbation to the scale factor, sigma the shear
# phi the Weyl potential, and eta the 3-curvature. Equations are in a general gauge,
# and are implemented in CAMB in the CDM frame (synchronous gauge, but using variables above).
# There are functions to convert into Newtonian and synchronous gauge metric variables

%matplotlib inline
import sys, platform, os
from matplotlib import pyplot as plt
import numpy as np
from IPython.display import display
import six
import warnings
warnings.filterwarnings('ignore', category=DeprecationWarning, module='.*/IPython/.*')
print('Using CAMB installed at %s'%(os.path.realpath(os.path.join(os.getcwd(),'..'))))
sys.path.insert(0,os.path.realpath(os.path.join(os.getcwd(),'..')))
import camb
from camb.symbolic import *
sympy.init_printing()
print('CAMB: %s, Sympy: %s'%(camb.__version__,sympy.__version__))

Using CAMB installed at c:\work\dist\git\camb\pycamb
CAMB: 0.1.6.1, Sympy: 1.0

In [2]:
display('background_eqs',background_eqs)
display('contraints',constraints)
display('var_subs',var_subs)
display('q_sub',q_sub)
display('pert_eqs',pert_eqs)
display('total_eqs',total_eqs)
#can use tot_eqs as combination of total_eqs + pert_eqs + background_eqs

'background_eqs'
$$\left [ \frac{d}{d t} a{\left (t \right )} = H{\left (t \right )} a{\left (t \right )}, \quad \frac{d}{d t} H{\left (t \right )} = - \frac{\kappa}{6} \left(3 P{\left (t \right )} + \rho{\left (t \right )}\right) a^{2}{\left (t \right )}, \quad \frac{d}{d t} \operatorname{exptau}{\left (t \right )} = \operatorname{visibility}{\left (t \right )}\right ]$$
'contraints'
$$\left [ Kf_{1} k^{3} \phi{\left (t \right )} + \frac{k \kappa}{2} \left(Kf_{1} \Pi{\left (t \right )} + \delta{\left (t \right )}\right) a^{2}{\left (t \right )} + \frac{3 \kappa}{2} H{\left (t \right )} a^{2}{\left (t \right )} q{\left (t \right )}, \quad k^{2} \eta{\left (t \right )} + 2 k H{\left (t \right )} z{\left (t \right )} - \kappa a^{2}{\left (t \right )} \delta{\left (t \right )}, \quad \frac{2 k^{2}}{3 a^{2}{\left (t \right )}} \left(- Kf_{1} \sigma{\left (t \right )} + z{\left (t \right )}\right) + \kappa q{\left (t \right )}, \quad - \frac{k}{3} z{\left (t \right )} + A{\left (t \right )} H{\left (t \right )} + \dot{h}{\left (t \right )}\right ]$$
'var_subs'
$$\left \{ \dot{h}{\left (t \right )} : \frac{1}{6 H{\left (t \right )}} \left(- k^{2} \eta{\left (t \right )} + \kappa a^{2}{\left (t \right )} \delta{\left (t \right )} - 6 A{\left (t \right )} H^{2}{\left (t \right )}\right), \quad \phi{\left (t \right )} : - \frac{\kappa a^{2}{\left (t \right )}}{2 Kf_{1} k^{3}} \left(Kf_{1} k \Pi{\left (t \right )} + k \delta{\left (t \right )} + 3 H{\left (t \right )} q{\left (t \right )}\right), \quad \sigma{\left (t \right )} : \frac{1}{2 Kf_{1} k^{2} H{\left (t \right )}} \left(k \left(- k^{2} \eta{\left (t \right )} + \kappa a^{2}{\left (t \right )} \delta{\left (t \right )}\right) + 3 \kappa H{\left (t \right )} a^{2}{\left (t \right )} q{\left (t \right )}\right), \quad z{\left (t \right )} : \frac{1}{2 k H{\left (t \right )}} \left(- k^{2} \eta{\left (t \right )} + \kappa a^{2}{\left (t \right )} \delta{\left (t \right )}\right)\right \}$$
'q_sub'
$$q{\left (t \right )} = \frac{2 k^{2}}{3 \kappa a^{2}{\left (t \right )}} \left(Kf_{1} \sigma{\left (t \right )} - \frac{1}{2 k H{\left (t \right )}} \left(- k^{2} \eta{\left (t \right )} + \kappa a^{2}{\left (t \right )} \delta{\left (t \right )}\right)\right)$$
'pert_eqs'
$$\left [ \frac{d}{d t} z{\left (t \right )} = Kf_{1} k A{\left (t \right )} - H{\left (t \right )} z{\left (t \right )} + \frac{3 \kappa}{2 k} \left(P{\left (t \right )} + \rho{\left (t \right )}\right) A{\left (t \right )} a^{2}{\left (t \right )} - \frac{\kappa}{2 k} \left(\delta{\left (t \right )} + 3 \delta_{P}{\left (t \right )}\right) a^{2}{\left (t \right )}, \quad \frac{d}{d t} \sigma{\left (t \right )} = k \left(A{\left (t \right )} + \phi{\left (t \right )}\right) - H{\left (t \right )} \sigma{\left (t \right )} - \frac{\Pi{\left (t \right )}}{2 k} \kappa a^{2}{\left (t \right )}, \quad \frac{d}{d t} \eta{\left (t \right )} = - \frac{1}{k} \left(2 K z{\left (t \right )} + 2 Kf_{1} k A{\left (t \right )} H{\left (t \right )} + \kappa a^{2}{\left (t \right )} q{\left (t \right )}\right), \quad \frac{d}{d t} \phi{\left (t \right )} = - H{\left (t \right )} \phi{\left (t \right )} + \frac{\kappa}{2 k^{2}} \left(k \left(P{\left (t \right )} + \rho{\left (t \right )}\right) \sigma{\left (t \right )} + k q{\left (t \right )} - H{\left (t \right )} \Pi{\left (t \right )} - \frac{d}{d t} \Pi{\left (t \right )}\right) a^{2}{\left (t \right )}\right ]$$
'total_eqs'
$$\left [ \frac{d}{d t} q{\left (t \right )} = - \frac{2 Kf_{1}}{3} k \Pi{\left (t \right )} - k \left(P{\left (t \right )} + \rho{\left (t \right )}\right) A{\left (t \right )} + k \delta_{P}{\left (t \right )} - 4 H{\left (t \right )} q{\left (t \right )}, \quad \frac{d}{d t} \delta{\left (t \right )} = - k q{\left (t \right )} - 3 \left(P{\left (t \right )} + \rho{\left (t \right )}\right) \dot{h}{\left (t \right )} - 3 \left(\delta{\left (t \right )} + \delta_{P}{\left (t \right )}\right) H{\left (t \right )}\right ]$$

Newtonian gauge variables $\Psi_N$ and $\Phi_N$ (not used in CAMB but may be useful) defined for metrix sign choices so flat metric is $$ds^2 = a(\eta)^2\left( (1+2\Psi_N)d\eta^2 - (1-2\Phi_N)\delta_{ij}dx^idx^j\right)$$ (default, as defined by Ma and Bertschinger, number count and 21cm papers, etc.) Definitions are also valid for non-flat.

The alternative definition $$ds^2 = a(\eta)^2\left( (1+2\Psi_N)d\eta^2 - (1+2\Phi_N)\delta_{ij}dx^idx^j\right)$$ is used by Hu et.al., lensing review, etc, corresponding to a sign change in $\Phi_N$.

In [3]:
#Relatations for going to/from the Newtonian gauge
display('Newt_vars',Newt_vars)
display('Newtonian_subs',Newtonian_subs)

'Newt_vars'
$$\left [ \Psi_{N}{\left (t \right )} = - A{\left (t \right )} + \frac{1}{k} \left(H{\left (t \right )} \sigma{\left (t \right )} + \frac{d}{d t} \sigma{\left (t \right )}\right), \quad \Phi_{N}{\left (t \right )} = - \frac{1}{k} H{\left (t \right )} \sigma{\left (t \right )} - \frac{\eta{\left (t \right )}}{2 Kf_{1}}\right ]$$
'Newtonian_subs'
$$\left [ A{\left (t \right )} = - \Psi_{N}{\left (t \right )}, \quad \frac{d^{2}}{d t^{2}} \sigma{\left (t \right )} = 0, \quad \frac{d}{d t} \sigma{\left (t \right )} = 0, \quad \sigma{\left (t \right )} = 0, \quad \phi{\left (t \right )} = \frac{1}{2} \Phi_{N}{\left (t \right )} + \frac{1}{2} \Psi_{N}{\left (t \right )}, \quad \eta{\left (t \right )} = - 2 Kf_{1} \Phi_{N}{\left (t \right )}, \quad \dot{h}{\left (t \right )} = - \frac{d}{d t} \Phi_{N}{\left (t \right )}, \quad z{\left (t \right )} = \frac{1}{k} \left(- 3 H{\left (t \right )} \Psi_{N}{\left (t \right )} - 3 \frac{d}{d t} \Phi_{N}{\left (t \right )}\right)\right ]$$
In [4]:
#e.g. get the Newtonian gauge equation for diff(Phi,t) + H*Psi
#[eq 23b of Ma and Bertschinger,
# noting that their (rho+P)theta = k q since q = (rho+P)v]
newtonian_gauge(subs(tot_eqs,subs(Newtonian_var_subs,
(diff(Phi_N,t) + H*Psi_N)).simplify().doit()))

Out[4]:
$$\frac{\kappa}{2 k} a^{2}{\left (t \right )} q{\left (t \right )}$$
In [5]:
#CDM frame is used by CAMB, corresponding to zero acceleration or CDM velocity
display('cdm_subs', cdm_subs)


'cdm_subs'
$$\left [ \frac{d}{d t} A{\left (t \right )} = 0, \quad A{\left (t \right )} = 0, \quad \frac{d}{d t} \operatorname{v_{c}}{\left (t \right )} = 0, \quad \operatorname{v_{c}}{\left (t \right )} = 0\right ]$$

Define synchonous gauge variables in Ma and Bertschinger notation (generalized to non-flat)

In terms of Hu et al variables $h_L+ h_T/3 = \eta_s$ and $h_L = -h_s/6$

In [6]:
#General gauge-invariant form of synchronous gauge metric variables
display('synchronous_vars',synchronous_vars)

'synchronous_vars'
$$\left [ \eta_{s}{\left (t \right )} = \frac{1}{Kf_{1}} \left(\frac{Kf_{1}}{k} H{\left (t \right )} \operatorname{v_{c}}{\left (t \right )} - \frac{1}{2} \eta{\left (t \right )}\right), \quad \dot{h}_{s}{\left (t \right )} = 6 A{\left (t \right )} H{\left (t \right )} + 6 \dot{h}{\left (t \right )} + \frac{6}{k} \left(\frac{Kf_{1} k^{2}}{3} + \frac{\kappa}{2} \left(P{\left (t \right )} + \rho{\left (t \right )}\right) a^{2}{\left (t \right )}\right) \operatorname{v_{c}}{\left (t \right )}\right ]$$
In [7]:
#To convert from general to synchronous variables can use these
display('synchronous_subs',synchronous_subs)

'synchronous_subs'
$$\left [ \eta{\left (t \right )} = - 2 Kf_{1} \eta_{s}{\left (t \right )}, \quad \dot{h}{\left (t \right )} = \frac{1}{6} \dot{h}_{s}{\left (t \right )}, \quad \phi{\left (t \right )} = - \frac{\kappa a^{2}{\left (t \right )}}{2 Kf_{1} k^{3}} \left(Kf_{1} k \Pi{\left (t \right )} + k \delta{\left (t \right )} + 3 H{\left (t \right )} q{\left (t \right )}\right), \quad z{\left (t \right )} = \frac{\dot{h}_{s}{\left (t \right )}}{2 k}, \quad \sigma{\left (t \right )} = \frac{1}{2 k} \left(\dot{h}_{s}{\left (t \right )} + 6 \frac{d}{d t} \eta_{s}{\left (t \right )}\right), \quad \frac{d}{d t} A{\left (t \right )} = 0, \quad A{\left (t \right )} = 0, \quad \frac{d}{d t} \operatorname{v_{c}}{\left (t \right )} = 0, \quad \operatorname{v_{c}}{\left (t \right )} = 0\right ]$$
In [8]:
#Alternative pure-metric expression for phi
Eq(phi,subs(Eq(Pi,-(3*diff(eta_s,t,t)+diff(hdot_s,t)/2 + 2*(3*H*diff(eta_s,t) + H*hdot_s/2) - k**2*eta_s)/kappa/a**2),
synchronous_gauge(subs(solve(constraints[1:],[delta,q,A]), phi_sub).rhs)).simplify())

Out[8]:
$$\phi{\left (t \right )} = \frac{1}{4 k^{2}} \left(2 k^{2} \eta_{s}{\left (t \right )} + \frac{d}{d t} \dot{h}_{s}{\left (t \right )} + 6 \frac{d^{2}}{d t^{2}} \eta_{s}{\left (t \right )}\right)$$
In [9]:
#Check the four synchronous gauge equations
synchronous_gauge(subs(var_subs,subs(synchronous_vars, K_fac*k**2*eta_s-H*hdot_s/2)).simplify())

Out[9]:
$$- \frac{\kappa}{2} a^{2}{\left (t \right )} \delta{\left (t \right )}$$
In [10]:
synchronous_gauge(subs(pert_eqs,subs(synchronous_vars, -K/2*hdot_s/k**2 + K_fac*diff(eta_s,t)).doit()))

Out[10]:
$$\frac{\kappa}{2 k} a^{2}{\left (t \right )} q{\left (t \right )}$$
In [11]:
synchronous_gauge(subs(K_sub,subs(Friedmann,subs(var_subs,subs(total_eqs,subs(background_eqs,subs(pert_eqs,subs(var_subs,
subs(synchronous_vars, -(diff(hdot_s,t) + H*hdot_s)/6).doit()).doit())).simplify().expand())
).simplify().expand())).simplify())

Out[11]:
$$\frac{\kappa}{6} \left(\delta{\left (t \right )} + 3 \delta_{P}{\left (t \right )}\right) a^{2}{\left (t \right )}$$
In [12]:
#Seems to be a factor of 2 missing in last line of eq A8 of Hu et al.
synchronous_gauge(subs(q_sub,subs(K_sub,subs(Friedmann_Kfac_subs,subs(var_subs,
subs(tot_eqs+background_eqs, subs(var_subs,subs(pert_eqs,subs(synchronous_vars,
3*diff(eta_s,t,t)+diff(hdot_s,t)/2 + 2*(3*H*diff(eta_s,t) + H*hdot_s/2) - k**2*eta_s).doit()).doit()).doit()))
.simplify()))))

Out[12]:
$$- \kappa \Pi{\left (t \right )} a^{2}{\left (t \right )}$$
In [13]:
#Fluid components
display('density_eqs',density_eqs)
display('delta_eqs', delta_eqs)
display('vel_eqs',vel_eqs)
#can use component_eqs as combination of density_eqs + delta_eqs + vel_eqs

'density_eqs'
$$\left [ \frac{d}{d t} \rho_{b}{\left (t \right )} = - 3 \left(\operatorname{p_{b}}{\left (t \right )} + \rho_{b}{\left (t \right )}\right) H{\left (t \right )}, \quad \frac{d}{d t} \rho_{c}{\left (t \right )} = - 3 H{\left (t \right )} \rho_{c}{\left (t \right )}, \quad \frac{d}{d t} \rho_{g}{\left (t \right )} = - 4 H{\left (t \right )} \rho_{g}{\left (t \right )}, \quad \frac{d}{d t} \rho_{r}{\left (t \right )} = - 4 H{\left (t \right )} \rho_{r}{\left (t \right )}, \quad \frac{d}{d t} \rho_{\nu}{\left (t \right )} = - 3 \left(\operatorname{p_{\nu}}{\left (t \right )} + \rho_{\nu}{\left (t \right )}\right) H{\left (t \right )}, \quad \frac{d}{d t} \rho_{de}{\left (t \right )} = - 3 \left(\operatorname{w_{de}}{\left (t \right )} + 1\right) H{\left (t \right )} \rho_{de}{\left (t \right )}\right ]$$
'delta_eqs'
$$\left [ \frac{d}{d t} \Delta_{r}{\left (t \right )} = - k \operatorname{q_{r}}{\left (t \right )} - 4 \dot{h}{\left (t \right )}, \quad \frac{d}{d t} \Delta_{g}{\left (t \right )} = - k \operatorname{q_{g}}{\left (t \right )} - 4 \dot{h}{\left (t \right )}, \quad \frac{d}{d t} \Delta_{b}{\left (t \right )} = \left(k \operatorname{v_{b}}{\left (t \right )} + 3 \dot{h}{\left (t \right )}\right) \left(- \frac{\operatorname{p_{b}}{\left (t \right )}}{\rho_{b}{\left (t \right )}} - 1\right) + \left(- 3 \operatorname{c^{2}_{sb}}{\left (t \right )} + \frac{3 \operatorname{p_{b}}{\left (t \right )}}{\rho_{b}{\left (t \right )}}\right) \Delta_{b}{\left (t \right )} H{\left (t \right )}, \quad \frac{d}{d t} \Delta_{c}{\left (t \right )} = - k \operatorname{v_{c}}{\left (t \right )} - 3 \dot{h}{\left (t \right )}, \quad \frac{d}{d t} \Delta_{\nu}{\left (t \right )} = - k \operatorname{q_{\nu}}{\left (t \right )} + \left(- \frac{3 \operatorname{p_{\nu}}{\left (t \right )}}{\rho_{\nu}{\left (t \right )}} - 3\right) \dot{h}{\left (t \right )} + 3 \left(- \Delta_{P \nu}{\left (t \right )} + \frac{\Delta_{\nu}{\left (t \right )} \operatorname{p_{\nu}}{\left (t \right )}}{\rho_{\nu}{\left (t \right )}}\right) H{\left (t \right )}, \quad \frac{d}{d t} \Delta_{de}{\left (t \right )} = - k \left(\operatorname{w_{de}}{\left (t \right )} + 1\right) \operatorname{v_{de}}{\left (t \right )} - 3 \left(\Delta_{de}{\left (t \right )} + \frac{3}{k} \left(\operatorname{w_{de}}{\left (t \right )} + 1\right) H{\left (t \right )} \operatorname{v_{de}}{\left (t \right )}\right) \left(\hat{c}^{2}_{sde}{\left (t \right )} - \operatorname{w_{de}}{\left (t \right )}\right) H{\left (t \right )} + \left(- 3 \operatorname{w_{de}}{\left (t \right )} - 3\right) \dot{h}{\left (t \right )} - \frac{3}{k} H{\left (t \right )} \operatorname{v_{de}}{\left (t \right )} \frac{d}{d t} \operatorname{w_{de}}{\left (t \right )}\right ]$$
'vel_eqs'
$$\left [ \frac{d}{d t} \operatorname{q_{r}}{\left (t \right )} = - \frac{2 Kf_{1}}{3} k \pi_{r}{\left (t \right )} - \frac{4 k}{3} A{\left (t \right )} + \frac{k}{3} \Delta_{r}{\left (t \right )}, \quad \frac{d}{d t} \operatorname{q_{g}}{\left (t \right )} = - \frac{2 Kf_{1}}{3} k \pi_{g}{\left (t \right )} - \frac{4 k}{3} A{\left (t \right )} + \frac{k}{3} \Delta_{g}{\left (t \right )} + \left(- \operatorname{q_{g}}{\left (t \right )} + \frac{4}{3} \operatorname{v_{b}}{\left (t \right )}\right) \operatorname{opacity}{\left (t \right )}, \quad \frac{d}{d t} \operatorname{v_{b}}{\left (t \right )} = - k A{\left (t \right )} + \frac{1}{\operatorname{p_{b}}{\left (t \right )} + \rho_{b}{\left (t \right )}} \left(k \Delta_{b}{\left (t \right )} \operatorname{c^{2}_{sb}}{\left (t \right )} \rho_{b}{\left (t \right )} - \left(- \operatorname{q_{g}}{\left (t \right )} + \frac{4}{3} \operatorname{v_{b}}{\left (t \right )}\right) \operatorname{opacity}{\left (t \right )} \rho_{g}{\left (t \right )}\right) - H{\left (t \right )} \operatorname{v_{b}}{\left (t \right )} - \frac{\operatorname{v_{b}}{\left (t \right )} \frac{d}{d t} \operatorname{p_{b}}{\left (t \right )}}{\operatorname{p_{b}}{\left (t \right )} + \rho_{b}{\left (t \right )}}, \quad \frac{d}{d t} \operatorname{v_{c}}{\left (t \right )} = - k A{\left (t \right )} - H{\left (t \right )} \operatorname{v_{c}}{\left (t \right )}, \quad \frac{d}{d t} \operatorname{q_{\nu}}{\left (t \right )} = - \frac{k}{3} \left(2 Kf_{1} \pi_{\nu}{\left (t \right )} - 3 \Delta_{P \nu}{\left (t \right )}\right) - k \left(\frac{\operatorname{p_{\nu}}{\left (t \right )}}{\rho_{\nu}{\left (t \right )}} + 1\right) A{\left (t \right )} - \left(- \frac{3 \operatorname{p_{\nu}}{\left (t \right )}}{\rho_{\nu}{\left (t \right )}} + 1\right) H{\left (t \right )} \operatorname{q_{\nu}}{\left (t \right )}, \quad \frac{d}{d t} \operatorname{v_{de}}{\left (t \right )} = - k A{\left (t \right )} + \frac{k \Delta_{de}{\left (t \right )} \hat{c}^{2}_{sde}{\left (t \right )}}{\operatorname{w_{de}}{\left (t \right )} + 1} - \left(- 3 \hat{c}^{2}_{sde}{\left (t \right )} + 1\right) H{\left (t \right )} \operatorname{v_{de}}{\left (t \right )}\right ]$$
In [14]:
#e.g. can check we recover standard Newtonian gauge equations
#(note all equations above are valid in any frame)
newtonian_gauge(delta_eqs)

Out[14]:
$$\left [ \frac{d}{d t} \Delta_{r}{\left (t \right )} = - k \operatorname{q_{r}}{\left (t \right )} + 4 \frac{d}{d t} \Phi_{N}{\left (t \right )}, \quad \frac{d}{d t} \Delta_{g}{\left (t \right )} = - k \operatorname{q_{g}}{\left (t \right )} + 4 \frac{d}{d t} \Phi_{N}{\left (t \right )}, \quad \frac{d}{d t} \Delta_{b}{\left (t \right )} = - \frac{1}{\rho_{b}{\left (t \right )}} \left(\left(k \operatorname{v_{b}}{\left (t \right )} - 3 \frac{d}{d t} \Phi_{N}{\left (t \right )}\right) \left(\operatorname{p_{b}}{\left (t \right )} + \rho_{b}{\left (t \right )}\right) + 3 \left(\operatorname{c^{2}_{sb}}{\left (t \right )} \rho_{b}{\left (t \right )} - \operatorname{p_{b}}{\left (t \right )}\right) \Delta_{b}{\left (t \right )} H{\left (t \right )}\right), \quad \frac{d}{d t} \Delta_{c}{\left (t \right )} = - k \operatorname{v_{c}}{\left (t \right )} + 3 \frac{d}{d t} \Phi_{N}{\left (t \right )}, \quad \frac{d}{d t} \Delta_{\nu}{\left (t \right )} = \frac{1}{\rho_{\nu}{\left (t \right )}} \left(- k \operatorname{q_{\nu}}{\left (t \right )} \rho_{\nu}{\left (t \right )} - 3 \left(\Delta_{P \nu}{\left (t \right )} \rho_{\nu}{\left (t \right )} - \Delta_{\nu}{\left (t \right )} \operatorname{p_{\nu}}{\left (t \right )}\right) H{\left (t \right )} + 3 \left(\operatorname{p_{\nu}}{\left (t \right )} + \rho_{\nu}{\left (t \right )}\right) \frac{d}{d t} \Phi_{N}{\left (t \right )}\right), \quad \frac{d}{d t} \Delta_{de}{\left (t \right )} = \frac{1}{k} \left(k \left(- k \operatorname{v_{de}}{\left (t \right )} + 3 \frac{d}{d t} \Phi_{N}{\left (t \right )}\right) \left(\operatorname{w_{de}}{\left (t \right )} + 1\right) - 3 \left(k \Delta_{de}{\left (t \right )} + 3 \left(\operatorname{w_{de}}{\left (t \right )} + 1\right) H{\left (t \right )} \operatorname{v_{de}}{\left (t \right )}\right) \left(\hat{c}^{2}_{sde}{\left (t \right )} - \operatorname{w_{de}}{\left (t \right )}\right) H{\left (t \right )} - 3 H{\left (t \right )} \operatorname{v_{de}}{\left (t \right )} \frac{d}{d t} \operatorname{w_{de}}{\left (t \right )}\right)\right ]$$
In [15]:
#Can use make_frame_invariant function to get explicit gauge-invariant equations for Newtonian (or other gauge) quantities
delta_N, v_b_N, sigma_N, eta_N = make_frame_invariant([delta,v_b, sigma, eta], 'Newtonian')
delta_sync, v_b_sync, sigma_sync, eta_sync = make_frame_invariant([delta,v_b, sigma, eta], 'CDM')

display('Gaude-dependent:', [delta, v_b, sigma, eta])
display('Newtonian:', [delta_N, v_b_N, sigma_N, eta_N])
display('Synchronous (CDM frame):', [delta_sync, v_b_sync, sigma_sync, eta_sync])

'Gaude-dependent:'
$$\left [ \delta{\left (t \right )}, \quad \operatorname{v_{b}}{\left (t \right )}, \quad \sigma{\left (t \right )}, \quad \eta{\left (t \right )}\right ]$$
'Newtonian:'
$$\left [ \delta{\left (t \right )} - \frac{3}{k} \left(P{\left (t \right )} + \rho{\left (t \right )}\right) H{\left (t \right )} \sigma{\left (t \right )}, \quad \sigma{\left (t \right )} + \operatorname{v_{b}}{\left (t \right )}, \quad 0, \quad \frac{2 Kf_{1}}{k} H{\left (t \right )} \sigma{\left (t \right )} + \eta{\left (t \right )}\right ]$$
'Synchronous (CDM frame):'
$$\left [ \delta{\left (t \right )} + \frac{3}{k} \left(P{\left (t \right )} + \rho{\left (t \right )}\right) H{\left (t \right )} \operatorname{v_{c}}{\left (t \right )}, \quad \operatorname{v_{b}}{\left (t \right )} - \operatorname{v_{c}}{\left (t \right )}, \quad \sigma{\left (t \right )} + \operatorname{v_{c}}{\left (t \right )}, \quad - \frac{2 Kf_{1}}{k} H{\left (t \right )} \operatorname{v_{c}}{\left (t \right )} + \eta{\left (t \right )}\right ]$$
In [16]:
def show_gauges(x):
print('Newtonian gauge version:')
display(newtonian_gauge(x))
print('CDM frame version:')
display(cdm_gauge(x))
print('Synchronous gauge variable version:')
display(synchronous_gauge(x))

def check_equation(camb_eq, view = True, p_b_zero=False):
if view:
show_gauges(camb_eq)
res = simplify(subs(background_eqs +component_eqs,camb_eq.simplify()))
res = subs(var_subs,res).simplify()
res = subs(Friedmann_Kfac_subs, res).simplify()
if p_b_zero: res = subs(Eq(p_b,0),res).doit().simplify()
if res==0:
print('OK')
else:
print('Non-zero, equal to:')
display(res)

#For example CAMB implements the equation d Delta_g/dt= -k*(4/3*z+qg) in synchronous gauge. Check this works.
Delta_g_sync, z_sync, q_g_sync = make_frame_invariant([Delta_g, z, q_g], 'CDM')
check_equation(diff(Delta_g_sync,t) + k*(4*z_sync/3 + q_g_sync))

Newtonian gauge version:

$$\frac{1}{k} \left(\frac{k^{2}}{3} \left(4 Kf_{1} \operatorname{v_{c}}{\left (t \right )} + 3 \operatorname{q_{g}}{\left (t \right )} - 4 \operatorname{v_{c}}{\left (t \right )}\right) - 4 k \left(H{\left (t \right )} \Psi_{N}{\left (t \right )} + \frac{d}{d t} \Phi_{N}{\left (t \right )}\right) + k \frac{d}{d t} \Delta_{g}{\left (t \right )} + 2 \kappa \left(P{\left (t \right )} + \rho{\left (t \right )}\right) a^{2}{\left (t \right )} \operatorname{v_{c}}{\left (t \right )} + 4 H{\left (t \right )} \frac{d}{d t} \operatorname{v_{c}}{\left (t \right )} + 4 \operatorname{v_{c}}{\left (t \right )} \frac{d}{d t} H{\left (t \right )}\right)$$
CDM frame version:

$$k \operatorname{q_{g}}{\left (t \right )} + \frac{4 k}{3} z{\left (t \right )} + \frac{d}{d t} \Delta_{g}{\left (t \right )}$$
Synchronous gauge variable version:

$$k \operatorname{q_{g}}{\left (t \right )} + \frac{2}{3} \dot{h}_{s}{\left (t \right )} + \frac{d}{d t} \Delta_{g}{\left (t \right )}$$
OK

In [17]:
csq_b_sync, Delta_b_sync = make_frame_invariant([csq_b, Delta_b], 'CDM')
check_equation(diff(v_b_sync,t)+H*v_b_sync
-csq_b_sync*k*Delta_b_sync+rho_g/rho_b*opacity*(4*v_b_sync/3-q_g_sync),p_b_zero=True)

Newtonian gauge version:

$$\frac{1}{\rho_{b}{\left (t \right )}} \left(- k \Delta_{b}{\left (t \right )} \operatorname{c^{2}_{sb}}{\left (t \right )} \rho_{b}{\left (t \right )} - \frac{1}{3} \left(3 \operatorname{q_{g}}{\left (t \right )} - 4 \operatorname{v_{b}}{\left (t \right )}\right) \operatorname{opacity}{\left (t \right )} \rho_{g}{\left (t \right )} + \left(\left(\operatorname{v_{b}}{\left (t \right )} - \operatorname{v_{c}}{\left (t \right )}\right) H{\left (t \right )} + \frac{d}{d t} \operatorname{v_{b}}{\left (t \right )} - \frac{d}{d t} \operatorname{v_{c}}{\left (t \right )}\right) \rho_{b}{\left (t \right )} + \operatorname{v_{c}}{\left (t \right )} \frac{d}{d t} \operatorname{p_{b}}{\left (t \right )}\right)$$
CDM frame version:

$$\frac{1}{\rho_{b}{\left (t \right )}} \left(- \frac{1}{3} \left(3 \operatorname{q_{g}}{\left (t \right )} - 4 \operatorname{v_{b}}{\left (t \right )}\right) \operatorname{opacity}{\left (t \right )} \rho_{g}{\left (t \right )} + \left(- k \Delta_{b}{\left (t \right )} \operatorname{c^{2}_{sb}}{\left (t \right )} + H{\left (t \right )} \operatorname{v_{b}}{\left (t \right )} + \frac{d}{d t} \operatorname{v_{b}}{\left (t \right )}\right) \rho_{b}{\left (t \right )}\right)$$
Synchronous gauge variable version:

$$\frac{1}{\rho_{b}{\left (t \right )}} \left(- \frac{1}{3} \left(3 \operatorname{q_{g}}{\left (t \right )} - 4 \operatorname{v_{b}}{\left (t \right )}\right) \operatorname{opacity}{\left (t \right )} \rho_{g}{\left (t \right )} + \left(- k \Delta_{b}{\left (t \right )} \operatorname{c^{2}_{sb}}{\left (t \right )} + H{\left (t \right )} \operatorname{v_{b}}{\left (t \right )} + \frac{d}{d t} \operatorname{v_{b}}{\left (t \right )}\right) \rho_{b}{\left (t \right )}\right)$$
OK

In [18]:
#can see how a synchronous gauge variable used in CAMB can be obtained in Newtonian gauge, e.g.
show_gauges(eta_sync)

Newtonian gauge version:

$$- \frac{2 Kf_{1}}{k} \left(k \Phi_{N}{\left (t \right )} + H{\left (t \right )} \operatorname{v_{c}}{\left (t \right )}\right)$$
CDM frame version:

$$\eta{\left (t \right )}$$
Synchronous gauge variable version:

$$- 2 Kf_{1} \eta_{s}{\left (t \right )}$$
In [19]:
Delta_g_N = make_frame_invariant(Delta_g, 'Newtonian')
show_gauges(Delta_g_N)

Newtonian gauge version:

$$\Delta_{g}{\left (t \right )}$$
CDM frame version:

$$\Delta_{g}{\left (t \right )} - \frac{4}{k} H{\left (t \right )} \sigma{\left (t \right )}$$
Synchronous gauge variable version:

$$\frac{1}{k^{2}} \left(k^{2} \Delta_{g}{\left (t \right )} - 2 \left(\dot{h}_{s}{\left (t \right )} + 6 \frac{d}{d t} \eta_{s}{\left (t \right )}\right) H{\left (t \right )}\right)$$
In [20]:
#Relations between components and totals
display('tot_subs',tot_subs)
display('tot_pert_subs',tot_pert_subs)

'tot_subs'
$$\left [ \rho{\left (t \right )} = \rho_{b}{\left (t \right )} + \rho_{c}{\left (t \right )} + \rho_{de}{\left (t \right )} + \rho_{g}{\left (t \right )} + \rho_{\nu}{\left (t \right )} + \rho_{r}{\left (t \right )}, \quad P{\left (t \right )} = \operatorname{p_{b}}{\left (t \right )} + \operatorname{p_{\nu}}{\left (t \right )} + \rho_{de}{\left (t \right )} \operatorname{w_{de}}{\left (t \right )} + \frac{1}{3} \rho_{g}{\left (t \right )} + \frac{1}{3} \rho_{r}{\left (t \right )}\right ]$$
'tot_pert_subs'
$$\left [ \Pi{\left (t \right )} = \pi_{g}{\left (t \right )} \rho_{g}{\left (t \right )} + \pi_{\nu}{\left (t \right )} \rho_{\nu}{\left (t \right )} + \pi_{r}{\left (t \right )} \rho_{r}{\left (t \right )}, \quad \delta{\left (t \right )} = \Delta_{b}{\left (t \right )} \rho_{b}{\left (t \right )} + \Delta_{c}{\left (t \right )} \rho_{c}{\left (t \right )} + \Delta_{de}{\left (t \right )} \rho_{de}{\left (t \right )} + \Delta_{g}{\left (t \right )} \rho_{g}{\left (t \right )} + \Delta_{\nu}{\left (t \right )} \rho_{\nu}{\left (t \right )} + \Delta_{r}{\left (t \right )} \rho_{r}{\left (t \right )}, \quad q{\left (t \right )} = \left(\operatorname{p_{b}}{\left (t \right )} + \rho_{b}{\left (t \right )}\right) \operatorname{v_{b}}{\left (t \right )} + \left(\operatorname{w_{de}}{\left (t \right )} + 1\right) \rho_{de}{\left (t \right )} \operatorname{v_{de}}{\left (t \right )} + \operatorname{q_{g}}{\left (t \right )} \rho_{g}{\left (t \right )} + \operatorname{q_{\nu}}{\left (t \right )} \rho_{\nu}{\left (t \right )} + \operatorname{q_{r}}{\left (t \right )} \rho_{r}{\left (t \right )} + \rho_{c}{\left (t \right )} \operatorname{v_{c}}{\left (t \right )}, \quad \delta_{P}{\left (t \right )} = \left(\Delta_{de}{\left (t \right )} \hat{c}^{2}_{sde}{\left (t \right )} + \frac{3}{k} \left(\operatorname{w_{de}}{\left (t \right )} + 1\right) \left(\hat{c}^{2}_{sde}{\left (t \right )} - \operatorname{w_{de}}{\left (t \right )} + \frac{\frac{d}{d t} \operatorname{w_{de}}{\left (t \right )}}{3 \left(\operatorname{w_{de}}{\left (t \right )} + 1\right) H{\left (t \right )}}\right) H{\left (t \right )} \operatorname{v_{de}}{\left (t \right )}\right) \rho_{de}{\left (t \right )} + \Delta_{P \nu}{\left (t \right )} \rho_{\nu}{\left (t \right )} + \Delta_{b}{\left (t \right )} \operatorname{c^{2}_{sb}}{\left (t \right )} \rho_{b}{\left (t \right )} + \frac{1}{3} \Delta_{g}{\left (t \right )} \rho_{g}{\left (t \right )} + \frac{1}{3} \Delta_{r}{\left (t \right )} \rho_{r}{\left (t \right )}\right ]$$
In [21]:
#First few equations in Boltzmann hierarchies for L>=2 (J photon, G neutrino)
display('hierarchies',hierarchies)

'hierarchies'
$$\left [ \frac{d}{d t} \pi_{g}{\left (t \right )} = - \frac{k}{5} \left(3 \operatorname{J_{3}}{\left (t \right )} Kf_{2} - 2 \operatorname{q_{g}}{\left (t \right )}\right) + \frac{8 k}{15} \sigma{\left (t \right )} - \operatorname{opacity}{\left (t \right )} \pi_{g}{\left (t \right )} + \operatorname{opacity}{\left (t \right )} \operatorname{polter}{\left (t \right )}, \quad \frac{d}{d t} \pi_{r}{\left (t \right )} = - \frac{k}{5} \left(3 \operatorname{G_{3}}{\left (t \right )} Kf_{2} - 2 \operatorname{q_{r}}{\left (t \right )}\right) + \frac{8 k}{15} \sigma{\left (t \right )}, \quad \frac{d}{d t} \operatorname{E_{2}}{\left (t \right )} = - \frac{k Kf_{2}}{3} \operatorname{E_{3}}{\left (t \right )} - \operatorname{E_{2}}{\left (t \right )} \operatorname{opacity}{\left (t \right )} + \operatorname{opacity}{\left (t \right )} \operatorname{polter}{\left (t \right )}, \quad \frac{d}{d t} \operatorname{J_{3}}{\left (t \right )} = - \frac{k}{7} \left(4 \operatorname{J_{4}}{\left (t \right )} Kf_{3} - 3 \pi_{g}{\left (t \right )}\right) - \operatorname{J_{3}}{\left (t \right )} \operatorname{opacity}{\left (t \right )}, \quad \frac{d}{d t} \operatorname{G_{3}}{\left (t \right )} = - \frac{k}{7} \left(4 \operatorname{G_{4}}{\left (t \right )} Kf_{3} - 3 \pi_{r}{\left (t \right )}\right), \quad \frac{d}{d t} \operatorname{E_{3}}{\left (t \right )} = - \frac{k}{7} \left(- 3 \operatorname{E_{2}}{\left (t \right )} + 3 \operatorname{E_{4}}{\left (t \right )} Kf_{3}\right) - \operatorname{E_{3}}{\left (t \right )} \operatorname{opacity}{\left (t \right )}, \quad \frac{d}{d t} \operatorname{J_{4}}{\left (t \right )} = - \frac{k}{9} \left(- 4 \operatorname{J_{3}}{\left (t \right )} + 5 \operatorname{J_{5}}{\left (t \right )} Kf_{4}\right) - \operatorname{J_{4}}{\left (t \right )} \operatorname{opacity}{\left (t \right )}, \quad \frac{d}{d t} \operatorname{G_{4}}{\left (t \right )} = - \frac{k}{9} \left(- 4 \operatorname{G_{3}}{\left (t \right )} + 5 \operatorname{G_{5}}{\left (t \right )} Kf_{4}\right), \quad \frac{d}{d t} \operatorname{E_{4}}{\left (t \right )} = - \frac{k}{9} \left(- 4 \operatorname{E_{3}}{\left (t \right )} + \frac{21 Kf_{4}}{5} \operatorname{E_{5}}{\left (t \right )}\right) - \operatorname{E_{4}}{\left (t \right )} \operatorname{opacity}{\left (t \right )}\right ]$$
In [22]:
# polter is the quadrupole source. Need its derivatives for line-of-sight solution.
polterdot = subs(hierarchies,diff(polter_t,t))
polterdot

Out[22]:
$$- \frac{k}{50} \left(3 \operatorname{J_{3}}{\left (t \right )} Kf_{2} - 2 \operatorname{q_{g}}{\left (t \right )}\right) - \frac{k Kf_{2}}{5} \operatorname{E_{3}}{\left (t \right )} + \frac{4 k}{75} \sigma{\left (t \right )} - \frac{3}{5} \operatorname{E_{2}}{\left (t \right )} \operatorname{opacity}{\left (t \right )} - \frac{1}{10} \operatorname{opacity}{\left (t \right )} \pi_{g}{\left (t \right )} + \frac{7}{10} \operatorname{opacity}{\left (t \right )} \operatorname{polter}{\left (t \right )}$$
In [23]:
polterddot = subs(phi_sub,diff(polterdot,t).subs(diff(sigma,t),dsigma).simplify()).simplify()
polterddot

Out[23]:
$$\frac{1}{150 Kf_{1} k} \left(8 Kf_{1} k^{3} A{\left (t \right )} - 8 Kf_{1} k^{2} H{\left (t \right )} \sigma{\left (t \right )} + Kf_{1} k \left(- 3 k \left(3 \frac{d}{d t} \operatorname{J_{3}}{\left (t \right )} Kf_{2} - 2 \frac{d}{d t} \operatorname{q_{g}}{\left (t \right )}\right) - 30 k \frac{d}{d t} \operatorname{E_{3}}{\left (t \right )} Kf_{2} - 4 \kappa \Pi{\left (t \right )} a^{2}{\left (t \right )} - 90 \operatorname{E_{2}}{\left (t \right )} \frac{d}{d t} \operatorname{opacity}{\left (t \right )} - 90 \operatorname{opacity}{\left (t \right )} \frac{d}{d t} \operatorname{E_{2}}{\left (t \right )} - 15 \operatorname{opacity}{\left (t \right )} \frac{d}{d t} \pi_{g}{\left (t \right )} + 105 \operatorname{opacity}{\left (t \right )} \frac{d}{d t} \operatorname{polter}{\left (t \right )} - 15 \pi_{g}{\left (t \right )} \frac{d}{d t} \operatorname{opacity}{\left (t \right )} + 105 \operatorname{polter}{\left (t \right )} \frac{d}{d t} \operatorname{opacity}{\left (t \right )}\right) - 4 \kappa \left(Kf_{1} k \Pi{\left (t \right )} + k \delta{\left (t \right )} + 3 H{\left (t \right )} q{\left (t \right )}\right) a^{2}{\left (t \right )}\right)$$
In [24]:
monopole_source, ISW, doppler, quadrupole_source = get_scalar_temperature_sources()

In [25]:
#These are definitions used in CAMB to get the various sources for the temperature
print(camb_fortran(dphi, 'phidot'))
print(camb_fortran(diff(polter_t,t), 'polterdot'))
print(camb_fortran(polterddot, 'polterddot'))
print(camb_fortran(2*diff(phi,t)*exptau,'ISW'))
print(camb_fortran(monopole_source, 'monopole_source'))
print(camb_fortran(doppler, 'doppler'))

phidot = (1.0d0/2.0d0)*(-adotoa*dgpi - 2*adotoa*k**2*phi + dgq*k &
-diff_rhopi + k*sigma*(gpres + grho))/k**2
polterdot = (1.0d0/10.0d0)*pigdot + (3.0d0/5.0d0)*Edot(2)
2.0d0/75.0d0*dgrho/Kf(1) + dopacity*(-1.0d0/10.0d0*pig + &
(7.0d0/10.0d0)*polter - 3.0d0/5.0d0*E(2)) &
-3.0d0/50.0d0*k*octgdot*Kf(2) + (1.0d0/25.0d0)*k*qgdot - &
1.0d0/5.0d0*k*Edot(3)*Kf(2) + opacity*(-1.0d0/10.0d0*pigdot + &
(7.0d0/10.0d0)*polterdot - 3.0d0/5.0d0*Edot(2))
ISW = 2*exptau*phidot
monopole_source = (1.0d0/4.0d0)*visibility*(-4*etak + k*(clxg + &
8*phi)*Kf(1))/(k*Kf(1))
doppler = (dvisibility*(sigma + vb) + visibility*(sigmadot + vbdot))/k
6*dvisibility*polterdot + visibility*(k**2*polter + &
3*polterddot))/k**2

In [26]:
#internal_consistency_checks()

In [27]:
pars = camb.set_params(H0=67.5, ombh2=0.022, omch2=0.122, As=2e-9, ns=0.95, tau=0.055)
from matplotlib import rcParams
rcParams.update( {'axes.labelsize': 14,
'font.size': 14,
'legend.fontsize': 14,
'xtick.labelsize': 13,
'ytick.labelsize': 13})
cl_label= r'$\ell(\ell+1)C_\ell/2\pi\quad [\mu {\rm K}^2]$'

In [28]:
#Example of plotting the time evolution of Newtonian gauge variables and the monopole sources

data= camb.get_background(pars)
conformal_times = np.linspace(1, 800, 300)
ks = [0.01,0.05]
Delta_g_N = make_frame_invariant(Delta_g, 'Newtonian')
display('Temperature monopole source in general', monopole_source)
display('Temperature monopole source in Newtonian gauge', newtonian_gauge(monopole_source))

ev = data.get_time_evolution(ks, conformal_times, ['delta_photon', Delta_g_N, Psi_N, monopole_source, ISW])
_, axs= plt.subplots(1,2, figsize=(14,6))
for i, ax in enumerate(axs):
ax.plot(conformal_times,ev[i,:, 0])
ax.plot(conformal_times,ev[i,:, 1])
ax.plot(conformal_times,ev[i,:, 2])

ax.plot(conformal_times,ev[i,:, 3]*1000)
ax.plot(conformal_times,ev[i,:, 4]*1000)

ax.set_title('$k= %s$'%ks[i])
ax.set_xlabel(r'$\eta/\rm{Mpc}$');
ax.set_xlim(conformal_times[0], conformal_times[-1]);
plt.legend([r'$\Delta_\gamma (CDM)$', r'$\Delta_\gamma (Newtonian)$',r'$\Psi_N$',
r'Monopole (x1000)', r'ISW (x1000)']);

'Temperature monopole source in general'
$$\left(\frac{1}{4} \Delta_{g}{\left (t \right )} + 2 \phi{\left (t \right )} + \frac{\eta{\left (t \right )}}{2 Kf_{1}}\right) \operatorname{visibility}{\left (t \right )}$$
'Temperature monopole source in Newtonian gauge'
$$\frac{1}{4} \left(\Delta_{g}{\left (t \right )} + 4 \Psi_{N}{\left (t \right )}\right) \operatorname{visibility}{\left (t \right )}$$
In [29]:
# You can also calculate power spectra for custom source functions.
# For example, let's split up the standard temperature result into the various sub-terms,
# and see how they contribute to the total

early_ISW = sympy.Piecewise( (ISW, 1/a-1> 30),(0, True))  #redshift > 30
late_ISW = ISW - early_ISW

names = ['mon','ISW','eISW','LISW','dop', 'Q']
source_names =names)

data= camb.get_results(pars)
dic = data.get_cmb_unlensed_scalar_array_dict(CMB_unit='muK')

In [30]:
ls =np.arange(dic['TxT'].shape[0])
plt.figure(figsize=(8,5))
plt.semilogx(ls,dic['monxmon'], color='C0')
plt.semilogx(ls,dic['LISWxLISW'], color='C1')
plt.semilogx(ls,dic['eISWxeISW'], ls='--', color='C1')
plt.semilogx(ls,dic['dopxdop'],color='C2')
plt.semilogx(ls,dic['QxQ'], color='C3')
plt.semilogx(ls,dic['TxT'], lw=2, color='k')
plt.xlabel('$\ell$')
plt.ylabel(cl_label)
plt.xlim(2, ls[-1])
plt.legend(['Monopole','Late ISW','Early ISW','Doppler','Quadrupole', 'Total'], loc = 'upper left');

In [31]:
# The monopole sources can also be decomposed in various different ways,
# e.g. in terms of comoving-frame quantities:

# a comoving photon density term (important on sub-horizon scales),
dg = make_frame_invariant(Delta_g/4*visibility, frame='comoving')

# plus a sum of potential and comoving curvature terms, and ISW
rest=make_frame_invariant(monopole_source - dg, frame=q) +ISW

names = ['dg','rest','mon']
camb.set_custom_scalar_sources([dg, rest, monopole_source+ISW], source_names =names)
data= camb.get_results(pars)
mondic = data.get_cmb_unlensed_scalar_array_dict(CMB_unit='muK')

ls =np.arange(dic['TxT'].shape[0])

plt.figure(figsize=(8,5))
plt.semilogx(ls,mondic['dgxdg'])
plt.semilogx(ls,mondic['restxrest'], color='r')
plt.semilogx(ls,mondic['restxdg'], color='m')
plt.semilogx(ls,mondic['monxmon'], color='k', lw=2)
plt.xlabel('$\ell$')
plt.ylabel(cl_label)
plt.axhline(0)

plt.xlim(2, ls[-1])
plt.legend([r'$\bar{\Delta}_\gamma$',r'$2\phi+\bar{\eta}/2 + {\rm ISW}$',
'Cross term', 'Monopole total']);


We now demonstrate how to calculate non-standard sources, e.g. as needed to calculate corrections to the standard lensing result from emission angle and time delay effects (see arXiv:1706.02673, as now implemeted in camb in the camb.emission_angle module).

In [32]:
#Various sources for emission angle and time delay, from recombination and reionization

chi = tau0-t #assume flat

emission_sources = {
'vperp' :  -(sigma+v_b)*visibility/k/chi,
'emit'  : 15*diff(polter*visibility,t)/8/k**2/chi,
'delay' : 15*diff(polter*visibility,t)/8/k**2/chi**2 *(tau0 - tau_maxvis),
'E'     : scalar_E_source}

sources ={}
for key, source in list(emission_sources.items()):
sources[key+'1'] = sympy.Piecewise((source,1/a-1>30),(0, True))  #recombination
sources[key+'2'] =  sympy.Piecewise((source,1/a-1<=30),(0, True))  #reionization

camb.set_custom_scalar_sources(sources, source_ell_scales={'E1':2,'E2':2})
data= camb.get_results(pars)
dic = data.get_cmb_unlensed_scalar_array_dict(CMB_unit='muK')

In [33]:
#Plot polariazation potentials

ls = np.arange(dic['ExE'].shape[0])
lfac = np.sqrt(ls*(ls+1))
plt.loglog(dic['ExE'])
plt.loglog(dic['emit1xemit1']*ls**2)
plt.loglog(dic['emit1xE1']*ls, color='C2')
plt.loglog(-dic['emit1xE1']*ls,ls ='--',color='C2')

plt.xlim(2,2000)
plt.ylim([1e-5,80])
plt.legend(['$C_\ell^E$',r'$\ell^2 C_\ell^{\psi_\zeta}$',r'$\ell C_\ell^{E\psi_\zeta}$'],
loc='upper left', frameon=False)
plt.xlabel(r'$\ell$')
plt.ylabel(cl_label);

In [34]:
#temperature velocity potentials

plt.plot(ls,ls**2*dic['vperp1xvperp1'],color='C0')
plt.plot(ls,ls*dic['vperp1xT'],color='C1')
plt.plot(ls,dic['TxT'],color='C2',ls='--')
#plt.plot(ls,np.exp(-2*0.0581)*lsampvel**2*vel_deltavis_cl_sp(lsampvel)*norm,ls=':',color='C0')
plt.axhline(0,color='k')
plt.xlim([0,2000])
plt.xlabel(r'$\ell$')
plt.ylabel(cl_label)
plt.legend([r'$\ell^2 C^{\psi_v}_\ell$',r'$\ell C^{T\psi_v}_\ell$',r'$C^{T}_\ell$'],
frameon=False);

In [35]:
#The camb.emission_angle module uses the above to calculate BB from emission angle and time delay
#e.g. compare BB from lensing, field rotation and emission angle/time delay
from camb import emission_angle
from camb import postborn
%time BB = emission_angle.get_emission_delay_BB(pars, lmax=3500)
%time Bom=postborn.get_field_rotation_BB(pars, lmax=3500)

Wall time: 35.3 s
Wall time: 1min 3s

In [36]:
plt.loglog(ls, data.get_lensed_scalar_cls(2500, CMB_unit = 'muK')[:,2])
plt.loglog(ls,Bom(ls))
plt.loglog(ls,BB(ls))
plt.xlim([10, 2000])
plt.ylabel(cl_label)
plt.xlabel(r'$\ell$')
plt.legend([r'$\kappa\kappa$', r'$\omega\omega$','emission + delay']);
plt.title('BB power spectra');