# Decompose a shear map into E and B modesΒΆ

from lenstools import ShearMap

import numpy as np
import matplotlib.pyplot as plt

from astropy.io import fits
from astropy.units import deg

shear_file_1 = fits.open(filename1)
gamma = shear_file_1[0].data.astype(np.float)
shear_file_1.close()

shear_file_2 = fits.open(filename2)
assert shear_file_2[0].header["ANGLE"] == angle
gamma = np.array((gamma,shear_file_2[0].data.astype(np.float)))
shear_file_2.close()

return angle*deg,gamma

l_edges = np.arange(200.0,50000.0,200.0)

l,EE,BB,EB = test_map.decompose(l_edges,keep_fourier=True)

assert l.shape == EE.shape == BB.shape == EB.shape

fig,ax = plt.subplots()
ax.plot(l,l*(l+1)*EE/(2.0*np.pi),label=r"$P_{EE}$")
ax.plot(l,l*(l+1)*BB/(2.0*np.pi),label=r"$P_{BB}$")
ax.plot(l,l*(l+1)*np.abs(EB)/(2.0*np.pi),label=r"$\vert P_{EB}\vert$")

ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel(r"$l$")
ax.set_ylabel(r"$l(l+1)P_l/2\pi$")

ax.legend(loc="upper left")

plt.savefig("EB.png")
plt.clf()

fig,ax = plt.subplots()
ax.plot(l,np.abs(EB)/np.sqrt(EE*BB))
ax.set_xlabel(r"$l$")
ax.set_ylabel(r"$P_{EB}/\sqrt{P_{EE}P_{BB}}$")

plt.savefig("EB_corr.png")


These are the $$P_{EE}(l),P_{BB}(l),P_{EB}(l)$$ power spectra

and this is their cross correlation