# plot_cmb.py - plots the output from CAMB
# by Bjoern Malte Schaefer, bjoern.malte.schaefer@uni-heidelberg.de

import numpy as np
import pylab as plt
import scipy.special as special
 

lmin = 0
lmax = 1800

plt.close()


data = np.loadtxt('camb_cmb.data')

l = data[lmin:lmax,0]
cu = data[lmin:lmax,1]

su = cu / l / (l+1.0) * 2.0 * np.pi


tmin = 0
tmax = 4


 
theta = np.linspace(tmin,tmax,num=1000)

delta = 0.01


def Cu(theta):
	
	Tu = 0
	TTu = 0
	i = 0


	for number in su:

		TTu = TTu + delta * (2 + i) * su[i] * special.j0((2 + i) * theta) / (2 * np.pi)	
		i = i + delta
	return TTu




x = Cu(theta)/Cu(0)



	




data = np.loadtxt('camb_lensed.data')

l = data[lmin:lmax,0]
cl = data[lmin:lmax,1]

sl = cl / l / (l+1.0) * 2.0 * np.pi


def Cl(theta):
	
	Tl = 0
	TTl = 0
	i = 0


	for number in sl:

		TTl = TTl + delta * (2 + i) * sl[i] * special.j0((2 + i) * theta) / (2 * np.pi)	
		i = i + delta
		
	return TTl




y = Cl(theta) / Cl(0)


plt.plot(theta,x,'b-',label='unlensed')
plt.plot(theta,y,'r-',label='lensed')

plt.legend(loc='upper right')
plt.xlabel('Theta')
plt.ylabel('Corellation of theta')

plt.show()

Delta = Cu(theta)/Cu(0)  - Cl(theta)/Cl(0) 


plt.plot(theta,Delta,'b-',label='differenz')


plt.show()



mmax = 200.0
mmin = -mmax

plt.close()




cov = np.array([[Cu(0),Cu(0.218)],[Cu(0.218),Cu(0)]])

ico = np.linalg.inv(cov)

r = cov[1,0] / np.sqrt(cov[0,0] * cov[1,1])
print(r)

x = y = np.linspace(mmin,mmax,101)
X,Y = np.meshgrid(x,y)
Z = X**2 * ico[0,0] + 2.0 * X * Y * ico[1,0] + Y**2 * ico[1,1]
mm = np.linspace(1,3,3)
levels = -2.0 * np.log(special.erfc(mm/np.sqrt(2.0)))

plt.contourf(X,Y,Z,levels)

plt.axis('equal')

plt.xlabel('$x$-axis')
plt.ylabel('$y$-axis')

plt.xlim([mmin,mmax])
plt.ylim([mmin,mmax])
plt.show()


cov2 = np.array([[Cl(0),Cl(0.218)],[Cl(0.218),Cl(0)]])

ico2 = np.linalg.inv(cov2)

r2 = cov2[1,0] / np.sqrt(cov2[0,0] * cov2[1,1])
print(r2)

x = y = np.linspace(mmin,mmax,101)
X,Y = np.meshgrid(x,y)
Z2 = X**2 * ico2[0,0] + 2.0 * X * Y * ico2[1,0] + Y**2 * ico2[1,1]

mm = np.linspace(1,3,3)
levels = -2.0 * np.log(special.erfc(mm/np.sqrt(2.0)))

plt.contourf(X,Y,Z2,levels)
plt.axis('equal')

plt.xlabel('$x$-axis')
plt.ylabel('$y$-axis')

plt.xlim([mmin,mmax])
plt.ylim([mmin,mmax])
plt.show()




mmax = 200
mmin = -mmax

x = y = np.linspace(mmin,mmax,101)
X,Y = np.meshgrid(x,y)
Z3 =10**5*((Z2-Z)**2)**0.5
mm = np.linspace(1,3,3)
levels = -2.0 * np.log(special.erfc(mm/np.sqrt(2.0)))

plt.contourf(X,Y,Z3,levels)
plt.axis('equal')

plt.xlabel('$x$-axis')
plt.ylabel('$y$-axis')

plt.xlim([mmin,mmax])
plt.ylim([mmin,mmax])
plt.show()