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




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

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 + delta2
        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)  
		x = x + delta
	return TTl


z = Cl2(theta)/Cl2(0)




plt.plot(theta,j,'b-',label='unlensed')
plt.plot(theta,k,'r-',label='lensed')

plt.legend(loc='upper right')
plt.xlabel('Theta [radian]')
plt.ylabel('Correlation of Theta')

plt.show()



plt.plot(theta,j,'b-',label='unlensed')
plt.plot(theta,z,'g-',label='lensed_appr')

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

plt.show()


plt.plot(theta,j,'b-',label='unlensed')
plt.plot(theta,k,'r-',label='lensed')
plt.plot(theta,z,'g-',label='lensed_appr')

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

plt.show()



plt.plot(theta,k,'r-',label='lensed')
plt.plot(theta,z,'g-',label='lensed_appr')

plt.legend(loc='upper right')
plt.xlabel('Theta [radian]')
plt.ylabel('Corellation of Theta')

plt.show()



Delta = (Cu(theta)/Cu(0))  - (Cl(theta)/Cl(0))
plt.plot(theta,Delta,'b-',label='difference_unlensed_lensed')

plt.legend(loc='upper right')
plt.xlabel('Theta [radian]')
plt.ylabel('Correlation of Theta')

plt.show()


Delta4 = Delta
Delta2 = (Cu(theta) / Cu(0))- (Cl2(theta) / Cl2(0))

plt.plot(theta,Delta4,'b-',label='difference_unlensed-lensed')
plt.plot(theta,Delta2,'g-',label='difference_unlensed-lensed_appr')

plt.legend(loc='upper right')
plt.xlabel('Theta [radian]')
plt.ylabel('Correlation of Theta')

plt.show()



Delta2 = (Cu(theta) / Cu(0)) - (Cl2(theta) / Cl2(0))

plt.plot(theta,Delta2,'b-',label='difference_unlensed-lensed_appr')

plt.legend(loc='upper right')
plt.xlabel('Theta [radian]')
plt.ylabel('Correlation of Theta')

plt.show()



Delta3 = (Cl(theta)/Cl(0)) - (Cl2(theta) / Cl2(0))
plt.plot(theta,Delta3,'b-',label='difference_lensed-lensed_appr')

plt.legend(loc='upper right')
plt.xlabel('Theta [radian]')
plt.ylabel('Correlation of Theta')

plt.show()