# 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
import math 

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=100)

delta = 0.01


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

	for number in su:
		TTu = TTu + delta * (2.0 + i) * su[i] * special.j0((2.0 + i) * theta) / (2.0 * np.pi)	
		i = i + delta
	return TTu


j = 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




k = Cl(theta) / Cl(0)

data = np.loadtxt('spectrum.data')
l = data[:,0]
c = data[:,1]
aux = c * l * (l+1.0)/2.0/np.pi
sk = c
delta2 = 0.01



def Cgl(theta):
	
	x = 0
	sigma = 0
	C2 = 0.0
	Cgltheta = 0.0

	for number in sk:
		y = x  + 10.0 
		#C2 = C2 + delta * y**3 * su[x] * special.jv(2,y*r) / (2 * np.pi)
		Cgltheta = Cgltheta + delta2 * (1 / 2.0 * np.pi) * (1/y) * 4 * sk[x] * special.j0(y * theta)
		x = x + delta
        return Cgltheta


def Cl2(theta):
	x = 0
        TTl = 0.0
        #TTl2=0
        
	for number in su:
		y = x + 2.0
		TTl = TTl + delta * y * su[x] * np.exp((-1) * y**2 * (Cgl(0)-Cgl(theta)) / 2.0) * special.j0(theta * y) / (2.0 * np.pi)  
		#TTl2 = TTl2 + delta * y * sl[x] * np.exp((-1)*y**2 * (Cgl0-Cgltheta) / 2) * (special.j0(y * theta) + (y**2 * C2 * special.jv(2,theta*y) / 2)) / (2 * np.pi)
		x = x + delta
	return TTl


z = Cl2(theta)/Cl2(0)

h = np.pi/1

mmax = 200.0
mmin = -mmax

plt.close()


cov = np.array([[Cu(0),Cu(h)],[Cu(h),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(h)],[Cl(h),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()




cov3 = np.array([[Cl2(0),Cl2(h)],[Cl2(h),Cl2(0)]])

ico3 = np.linalg.inv(cov3)

r3 = cov3[1,0] / np.sqrt(cov3[0,0] * cov3[1,1])
print(r3)

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

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()


mmax = 200
mmin = -mmax

x = y = np.linspace(mmin,mmax,101)
X,Y = np.meshgrid(x,y)
Z4 =(10**6)*((Z-Z2)**2)**0.5
mm = np.linspace(1,3,3)
levels = -2.0 * np.log(special.erfc(mm/np.sqrt(2.0)))
print(Z4)
plt.contourf(X,Y,Z4,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)
Z5 =(10**8)*((Z3-Z)**2)**0.5
mm = np.linspace(1,3,3)
levels = -2.0 * np.log(special.erfc(mm/np.sqrt(2.0)))
print(Z5)
plt.contourf(X,Y,Z5,levels)
plt.axis('equal')

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

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