from math import sqrt
from math import pi
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate

def sn_sq(vn,vf):
	# sn und vn beginne bei v1 und s1
	return 1/2*((vn/(2*pi*vf)+1)/np.sqrt((vn/(2*pi*vf)+1)**2-(vn/(2*pi*vf))**2)-1)

def calc_cm(vn,vf):
	"""
	cms beginnen bei c0 und gehen so weit wie möglich hoch.
	"""
	sn_sq_arr = sn_sq(vn,vf)
	cms = np.array([1,])
	for m in range(1,len(sn_sq_arr)+1):
		cms = np.append(cms,np.sum(cms[::-1]*sn_sq_arr[:m])/m)
	return cms

def calc_al(vn,vf):
	"""
	als beginnen bei 0 und gehen positiv soweit wie möglich hoch
	In der Summe muss mindestens 10 cms vorkommen
	"""
	cms = calc_cm(vn,vf)
	sns = sn_sq(vn,vf)
	factor = np.exp(-2*np.sum(sns/(np.arange(len(sns))+1)))
	als = [np.sum(cms[:]*cms[:]),]
	cmin = 10
	# Hiernochmal richtig cms betrachten.
	for l in range(1,len(cms)-cmin):
		als.append(np.sum(cms[:-l]*cms[l:]))
	return np.array(als)*factor

def calc_bl(vn,vf,ls):
	# l ist relativ zur fermi nf als state n_(nf+l)
	# Schnellermachen indem wir die summe vorher berechnen und dann subtrahieren.
	als = calc_al(vn,vf)
	ns = []
	for l in ls:
		n = 0
		for i in range(len(als)-l):
			n += als[abs(l+i)]
		ns.append(n)
	return np.array(ns)

def vn_rechteck(v0,vf,nc,length):
	print("anomalous dimension: ",2*sn_sq(v0,vf))
	return np.append(np.ones(nc)*v0,np.zeros(length-nc))

def vn_noninteract(length):
	return np.zeros(length)

def vn_exp(v0,vf,nc,length):
	print("anomalous dimension: ",2*sn_sq(v0,vf))
	return np.exp(-np.arange(length)/nc)*v0

def thermodynlim(u,l,alpha):
	return -np.sin(l*u)/u/(2*pi)*(4/(4+u**2))**(alpha/2)

def integrate50(a, b, N,l,alpha):
    x = b-np.linspace(a, b, N,endpoint=False)
    fx = thermodynlim(x,l,alpha)
    area = np.sum(fx)*(b-a)/N
    return area

#/(1-np.exp(1j*u))/(2*pi)
if __name__ == "__main__":

	plot_range = 4

	nc = 1000
	length = 10000
	vf = 1
	m = plot_range*nc

	vn = vn_noninteract(length)
	vn2 = vn_rechteck(10,vf,nc,length)
	vn3 = vn_rechteck(80,vf,nc,length)

	U = set(np.arange(-m,m+1)[::nc//4]) | set(np.arange(-m+nc,m-nc+1)[::nc//8]) | set(np.arange(-m+2*nc,m-2*nc+1)[::nc//32]) | set(np.arange(-m+3*nc,m-3*nc+1)[::nc//32]) | set(np.arange(-nc//8,nc//8))
	x = np.array(sorted(U))
	y = calc_bl(vn,vf,x)
	y2 = calc_bl(vn2,vf,x)
	y3 = calc_bl(vn3,vf,x)

	#y4 = np.array([2*integrate50(0, 1000, 100001,l,alpha) for l in x/nc])

	plt.figure()
	plt.plot(x/nc,y)
	plt.plot(x/nc,y2)
	plt.plot(x/nc,y3)

	plt.xlabel("$\\frac{k-k_f}{k_c}$",size=15)
	plt.ylabel("$<n_k>$",size=15)
	plt.grid(linestyle='--',linewidth=0.3)
	plt.tick_params(axis='x', labelsize=15)
	plt.tick_params(axis='y', labelsize=15)
	plt.tight_layout()
	plt.minorticks_on()
	plt.tick_params(axis='both',which='minor',length=3,width=1)
	plt.show()