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

lmin = 0
lmax = 1800

plt.close()

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

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) *( y) * 4.0 * sk[x] * special.j0(y * theta)
		x = x + delta2
        return Cgltheta
x = Cgl(0) - Cgl(theta)

plt.plot(theta,x,'b-',label='sigma^2')
plt.show()