#HarmonicOscillator1
#This program obtains the spectrum of the harmonic oscillator using the one-run approach 
#of the bootstrap method.

from mpmath import *
mp.dps = 50 

#making a grid
def MakingGrid(start,end,arr,step):
	x = start
	while x <= end:
		arr.append(x)
		x = fadd(x,step)
	return(arr)

#recursion for creating a moment sequence
def RecursionRelation(y,x,K): 
	x =[1,0,mpf(y),0]
	for s in range(4,2*K-1):
		if s%2 == 0:
			x.append(fadd(fdiv(fprod([2,y,(s-1),x[s-2]]),s), 
			fdiv(fprod([s-1,s-2,s-3,x[s-4]]),fmul(4,s))))
		else:
			x.append(0)
	return x
	
#creating a Hankel matrix from the moment sequence
def HankelMatrix(x):
	d = int((len(x) + 1)/2)
	H = zeros(d)
	for i in range(0,d):
		for j in range(0,d):
			if fmul(x[i+1],x[j+1]) != 0:
				H[i,j] = fdiv(x[i+j],fmul(x[i+1],x[j+1]))
			else:
				H[i,j] = x[i+j]
	return H

#creating a search space
N = power(10,4)
K = 50
E_min, E_max = 0, 10
length = fsub(E_max,E_min)
S = MakingGrid(E_min,E_max,[],fdiv(length,N))
S_reduced = MakingGrid(E_min,E_max,[],fdiv(length,N))

#reducing the search space
for i in range(0,len(S)):
	O = RecursionRelation(S[i],[],K)
	try:
		M = HankelMatrix(O)
		L = cholesky(M)
	except ValueError:
		S_reduced[i] = 0
print("\n Reduced search space: \n",S_reduced)

bounds = []
for i in range(1,len(S_reduced)-1):
	if i==1 and S_reduced[0]!=0:
		bounds.append(S_reduced[0])
	if i==len(S_reduced)-2 and S_reduced[-1]!=0:
		bounds.append(S_reduced[-1])	
	if ((S_reduced[i-1]==0) and (S_reduced[i]!=0)): 
		bounds.append(S_reduced[i])
	if ((S_reduced[i+1]==0) and (S_reduced[i]!=0)):
		bounds.append(S_reduced[i])				
print("\n Bounds for accepted regions: \n",bounds)

results = []
for i in range(0,int(len(bounds)/2)):
	if bounds[2*i]!=bounds[2*i+1]:
		results.append([fdiv(fadd(bounds[2*i],bounds[2*i+1]),2),
		fdiv(fsub(bounds2[2*i+1],bounds2[2*i]),2)])
	else:
		results.append([fdiv(fadd(bounds[2*i],bounds[2*i+1]),2),fdiv(length,N)])
print("\n Accepted energies are: \n",results)

#results
#[[mpf('0.49999999999999999999999999999999999999999999999993786'),
#mpf('0.0010000000000000000000000000000000000000000000000000012')],
#[mpf('1.5000000000000000000000000000000000000000000000003207'),
#mpf('0.0010000000000000000000000000000000000000000000000000012')],
#[mpf('2.5000000000000000000000000000000000000000000000007056'),
#mpf('0.0010000000000000000000000000000000000000000000000000012')],
#[mpf('3.5000000000000000000000000000000000000000000000010905'),
#mpf('0.0010000000000000000000000000000000000000000000000000012')],
#[mpf('4.4999999999999999999999999999999999999999999999988026'),
#mpf('0.0010000000000000000000000000000000000000000000000000012')],
#[mpf('5.499999999999999999999999999999999999999999999993842'),
#mpf('0.0010000000000000000000000000000000000000000000000000012')],
#[mpf('6.4999999999999999999999999999999999999999999999888813'),
#mpf('0.0010000000000000000000000000000000000000000000000000012')],
#[mpf('7.4999999999999999999999999999999999999999999999839206'),
#mpf('0.0010000000000000000000000000000000000000000000000000012')],
#[mpf('8.49999999999999999999999999999999999999999999997896'),
#mpf('0.0010000000000000000000000000000000000000000000000000012')],
#[mpf('9.4999999999999999999999999999999999999999999999739993'),
#mpf('0.0010000000000000000000000000000000000000000000000000012')]]

#HarmonicOscillator2
#This program obtains the spectrum of the harmonic oscillator using the subsequent 
#approach of the bootstrap method. Also, it tracks the convergence of the method.

import numpy as np
import scipy as sp
from scipy import odr
import matplotlib.pyplot as plt
from mpmath import *
mp.dps = 50 
plt.rcParams.update({
    "text.usetex": True,
    "font.family": "cmr"
})

#making a grid
def MakingGrid(start,end,arr,step):
	x = start
	while x <= end:
		arr.append(x)
		x = fadd(x,step)
	return(arr)

#recursion for creating a moment sequence
def RecursionRelation(y,x,K): 
	x =[1,0,mpf(y),0]
	for s in range(4,2*K-1):
		if s%2 == 0:
			x.append(fadd(fdiv(fprod([2,y,(s-1),x[s-2]]),s), 
			fdiv(fprod([s-1,s-2,s-3,x[s-4]]),fmul(4,s))))
		else:
			x.append(0)
	return x
	
#creating a Hankel matrix from the moment sequence 	
def HankelMatrix(x):
	d = int((len(x) + 1)/2)
	H = zeros(d)
	for i in range(0,d):
		for j in range(0,d):
			if fmul(x[i+1],x[j+1]) != 0:
				H[i,j] = fdiv(x[i+j],fmul(x[i+1],x[j+1]))
			else:
				H[i,j] = x[i+j]
	return H

#defining a fitting function	
def fitting_function(p,x):
	a,b,c = p
	return a*np.exp(-b*(x-10)**c)

#creating a search space
N = power(10,3)
K = 10
K_list = []
E_min, E_max = 0, 10
length = fsub(E_max,E_min)
S = [MakingGrid(E_min,E_max,[],fdiv(length,N))]
S_reduced = [MakingGrid(E_min,E_max,[],fdiv(length,N))]
results = []
convergence = []

while K<=50:
	#reducing the search space
	for j in range(0,len(S)):
		for i in range(0,len(S[j])):
			O = RecursionRelation(S[j][i],[],K)
			try:
				M = HankelMatrix(O)
				L = cholesky(M)
			except ValueError:
				S_reduced[j][i] = 0
	print(f"\n Reduced search space for K = {K}: \n",S_reduced)

	bounds1 = []
	for j in range(0,len(S_reduced)):
		for i in range(1,len(S_reduced[j])-1):
			if i==1 and S_reduced[j][0]!=0:
				bounds1.append(S_reduced[j][0])
			if i==len(S_reduced[j])-2 and S_reduced[j][-1]!=0:
				bounds1.append(S_reduced[j][-1])	
			if ((S_reduced[j][i-1]==0) and (S_reduced[j][i]!=0)): 
				bounds1.append(S_reduced[j][i])
			if ((S_reduced[j][i+1]==0) and (S_reduced[j][i]!=0)):
				bounds1.append(S_reduced[j][i])
			if ((S_reduced[j][i-1]==0) and (S_reduced[j][i]!=0) and 
			(S_reduced[j][i+1]==0)):
				results.append([S_reduced[j][i],fdiv(fsub(S[j][i+1],S[j][i-1]),2)])								
	print(f"\n Bounds for accepted regions for K = {K}: \n",bounds1)
	
	bounds2 = []
	for i in range(0,len(bounds1)-1,2):
		if bounds1[i] != bounds1[i+1]:
			bounds2.append(bounds1[i])
			bounds2.append(bounds1[i+1])
	print(f"\n Modified bounds for accepted regions for K = {K}: \n",bounds2)
	
	totwidth = 0
	if len(bounds2)>=2:
		for i in range(0,len(bounds2)-1,2):
			totwidth = fadd(totwidth,fsub(bounds2[i+1],bounds2[i]))
	convergence.append(totwidth)
		
	S = []; S_reduced = []
	for i in range(0,int(len(bounds2)/2)):
		S.append(MakingGrid(bounds2[2*i],bounds2[2*i+1],[],
		fdiv(fsub(bounds2[2*i+1],bounds2[2*i]),N)))
		S_reduced.append(MakingGrid(bounds2[2*i],bounds2[2*i+1],[],
		fdiv(fsub(bounds2[2*i+1],bounds2[2*i]),N)))
	K_list.append(K)	
	K += 5
	
for i in range(0,int(len(bounds2)/2)):
	results.append([fdiv(fadd(bounds2[2*i],bounds2[2*i+1]),2),
	fdiv(fsub(bounds2[2*i+1],bounds2[2*i]),2)])
print("\n Accepted energies are: \n",results)
print("\n Total width of accepted energies for given K: \n",convergence)

#fitting data for convergence
for i in range(0,len(convergence)):
	convergence[i] = float(str(convergence[i]))
	
odr_model = sp.odr.Model(fitting_function)
data = odr.Data(K_list,convergence)
fit = sp.odr.ODR(data,odr_model,beta0=[7,0.0001,3])
out = fit.run()
beta = out.beta
std = out.sd_beta
print("\n Parameters of the fit: ",beta,
"\n Standart deviation of the parameters of the fit: ",std)

#plotting
plt.plot(K_list,convergence,color = "black",marker = "o",linestyle = "None")
plt.plot(np.linspace(K_list[0],K_list[-1],50),
beta[0]*np.exp(-beta[1]*(np.linspace(K_list[0],K_list[-1],50)-10)**beta[2]), 
color = "red", linestyle = "--", linewidth = 1)
plt.legend(["Data","Fitted curve $ae^{-b(K-10)^c}$"],loc = "upper right", fontsize=14)
plt.title("Convergence of total interval width",fontsize=16)
plt.xlabel("$K$",fontsize=14)
plt.ylabel("Total interval width",fontsize=14)
plt.grid(linestyle = '-', linewidth = 0.5)
plt.savefig("ConvergenceHO.pdf", dpi=300)

#results
#[[mpf('0.49999999999999999999999999999999999999999999999999399'),
#mpf('0.010000000000000000000000000000000000000000000000000174')],
#[mpf('1.4999999999999999958400000000000000000000000000004858'),
#mpf('8.6929920000000000000000000000000011779974046389434164e-18')],
#[mpf('2.500000000000001156181811199999999999999999999999476'),
#mpf('0.0000000000000027174445055999999999999999999999999937265956758702698')],
#[mpf('3.4999999999998113703243775999999999999999999999967355'),
#mpf('0.00000000000037584655257599999999999999999999999998761753563479507')],
#[mpf('4.5000000000128640947149000000000000000000000000035599'),
#mpf('0.000000000027535164427980000000000000000000000000015047608097755')],
#[mpf('5.4999999994353146819329023999999999999999999999984901'),
#mpf('0.000000001317173035830220800000000000000000000000001378724333')],
#[mpf('6.5000000116722248883561663999999999999999999999972256'),
#mpf('0.000000043185842402698195199999999999999999999999996779702612')],
#[mpf('7.4999995904551134172627417599999999999999999999931808'),
#mpf('0.0000010149147534330808819200000000000000000000000563129298')],
#[mpf('8.5000055637842406778576959999999999999999999999938953'),
#mpf('0.000015477452348826660191999999999999999999999999996942449')],
#[mpf('9.4999348724498841208245759999999999999999999999769981'),
#mpf('0.00019103136524200255887359999999999999999999999972475485')]]

#convergence 
#[mpf('8.4299999999999999999999999999999999999999999999998982'),
#mpf('6.850359999999999999999999999999999999999999999997984'),
#mpf('5.683739599999999999999999999999999999999999999995754'),
#mpf('4.1545833661199999999999999999999999999999999999967668'),
#mpf('2.6920863477472799999999999999999999999999999999948565'),
#mpf('0.96519348255029567999999999999999999999999999999490432'),
#mpf('0.12174611420385617280000000000000000000000000000100281'),
#mpf('0.0089499965904711390019199999999999999999999999989644344'),
#mpf('0.00041513652654687589283423999999999999999999999954747678')]	

#fit
#[8.04903532e+00 3.80749926e-03 1.93571082e+00]
#[0.33395913 0.00286182 0.24733784]

#DoubleWell1.py
#This program obtains the picture of the reduced search space for the double-well
#potential using one–run approach of the bootstrap method.

import matplotlib.pyplot as plt
from mpmath import *
mp.dps = 50 
plt.rcParams.update({
    "text.usetex": True,
    "font.family": "cmr"
})

#making a grid
def MakingGrid(start,end,arr,step):
	x = start
	while x <= end:
		arr.append(x)
		x = fadd(x,step)
	return(arr)

#recursion for creating a moment sequence
def RecursionRelation(x,y,z,w,K): 
	x =[1,0,mpf(z),0,fdiv(fadd(fsub(fmul(2,mpf(y)),fdiv(1,fmul(16,mpf(w)))),mpf(z)),
	fmul(3,mpf(w))),0]
	for s in range(6,2*K-1):
		if s%2 == 0:
			x.append(fsum([fdiv(fprod([fsub(fmul(2,mpf(y)),fdiv(1,fmul(16,mpf(w)))),
			s-3,x[s-4]]),fmul(mpf(w),(s-1))),fdiv(fmul(s-2,x[s-2]),fprod([2,mpf(w),s-1])),
			fdiv(fprod([s-3,s-4,s-5,x[s-6]]),fprod([4,mpf(w),(s-1)]))]))
		else:
			x.append(0)
	return x
	
#creating a Hankel matrix from the moment sequence
def HankelMatrix(x):
	d = int((len(x) + 1)/2)
	H = zeros(d)
	for i in range(0,d):
		for j in range(0,d):
			if fmul(x[i+1],x[j+1]) != 0:
				H[i,j] = fdiv(x[i+j],fmul(x[i+1],x[j+1]))
			else:
				H[i,j] = x[i+j]
	return H

#classical limit
g = 0.05
N = 600
E_min, E_max = 0, 3
length1 = fsub(E_max,E_min)
S1 = MakingGrid(E_min,E_max,[],fdiv(length1,N))
mom = []
for i in range(1,len(S1)):
	if S1[i]<=(1/(32*g)):
		x1 = sqrt((1/(2*g)-sqrt(8*S1[i]/g))/2)
		x2 = sqrt((1/(2*g)+sqrt(8*S1[i]/g))/2)
	else:
		x1 = -sqrt((1/(2*g)+sqrt(8*S1[i]/g))/2)
		x2 = sqrt((1/(2*g)+sqrt(8*S1[i]/g))/2)	
	integrand1 = lambda x: 1/(sqrt(2*S1[i]-g*(x**2-1/(4*g))**2))
	integrand2 = lambda x: x**2/(sqrt(2*S1[i]-g*(x**2-1/(4*g))**2))
	T = 2*re(quad(integrand1,[x1,x2]))
	mom.append(2*re(quad(integrand2,[x1,x2])/T))
plt.plot(S1[1:],mom,color= "black",linestyle="--")

#creating a search space
K = 16
mom_min, mom_max = 0, mom[-1]
length2 = fsub(mom_max,mom_min)
S2 = MakingGrid(mom_min,mom_max,[],fdiv(length2,N))
S1_reduced = []
for i in range(0,len(S2)):
	S1_reduced.append(MakingGrid(E_min,E_max,[],fdiv(length1,N))) 

#reducing the search space
for i in range(0,len(S2)):
	for j in range(0,len(S1)):
		O = RecursionRelation([],S1[j],S2[i],g,K)
		try:
			M = HankelMatrix(O)
			L = cholesky(M)
		except ValueError:
			S1_reduced[i].remove(S1[j])		
print("\n Reduced search space: \n",S1_reduced)

#plotting
for i in range(0,len(S2)):
	plt.plot(S1_reduced[i],[S2[i]]*len(S1_reduced[i]),
	color = "red",marker = "o",linestyle = "None")
	
plt.legend(["Classical","Allowed regions"], loc="lower right",fontsize=14)
plt.title("Bootstrap at $g = 0.05$",fontsize=16)
plt.xlabel("$E$",fontsize=14)
plt.ylabel(r"$\langle x^2 \rangle$",fontsize=14)
plt.grid(linestyle = '-', linewidth = 0.5)
plt.savefig("S_reduced_DoubleWell1.pdf", dpi=300)

#DoubleWell2.py
#This program obtains the spectrum of the double-well using the one-run approach 
#of the bootstrap method.

import time
from mpmath import *
mp.dps = 50 

start = time.time()

#making a grid
def MakingGrid(start,end,arr,step):
	x = start
	while x <= end:
		arr.append(x)
		x = fadd(x,step)
	return(arr)

#recursion for creating a moment sequence
def RecursionRelation(x,y,z,w,K): 
	x =[1,0,mpf(z),0,fdiv(fadd(fsub(fmul(2,mpf(y)),fdiv(1,fmul(16,mpf(w)))),mpf(z)),
	fmul(3,mpf(w))),0]
	for s in range(6,2*K-1):
		if s%2 == 0:
			x.append(fsum([fdiv(fprod([fsub(fmul(2,mpf(y)),fdiv(1,fmul(16,mpf(w)))),
			s-3,x[s-4]]),fmul(mpf(w),(s-1))),fdiv(fmul(s-2,x[s-2]),fprod([2,mpf(w),s-1])),
			fdiv(fprod([s-3,s-4,s-5,x[s-6]]),fprod([4,mpf(w),(s-1)]))]))
		else:
			x.append(0)
	return x
	
#creating a Hankel matrix from the moment sequence
def HankelMatrix(x):
	d = int((len(x) + 1)/2)
	H = zeros(d)
	for i in range(0,d):
		for j in range(0,d):
			if fmul(x[i+1],x[j+1]) != 0:
				H[i,j] = fdiv(x[i+j],fmul(x[i+1],x[j+1]))
			else:
				H[i,j] = x[i+j]
	return H

#creating a search space
g = 0.037
N = 900
K = 18
E_min, E_max = 0, 1/(32*g)
length1 = fsub(E_max,E_min)
S1 = MakingGrid(E_min,E_max,[],fdiv(length1,N))
x1 = sqrt((1/(2*g)-sqrt(8*S1[1]/g))/2)
x2 = sqrt((1/(2*g)+sqrt(8*S1[1]/g))/2)	
integrand1 = lambda x: 1/(sqrt(2*S1[1]-g*(x**2-1/(4*g))**2))
integrand2 = lambda x: x**2/(sqrt(2*S1[1]-g*(x**2-1/(4*g))**2))
T = 2*re(quad(integrand1,[x1,x2]))
mom_max = 2*re(quad(integrand2,[x1,x2])/T); mom_min = 0
	
length2 = fsub(mom_max,mom_min)
S2 = MakingGrid(mom_min,mom_max,[],fdiv(length2,N))
S1_reduced = MakingGrid(E_min,E_max,[],fdiv(length1,N))
S2_reduced = []
for i in range(0,len(S1)):
	S2_reduced.append(MakingGrid(mom_min,mom_max,[],fdiv(length2,N))) 

#reducing the search space
for i in range(0,len(S1)):
	for j in range(0,len(S2)):
		O = RecursionRelation([],S1[i],S2[j],g,K)
		try:
			M = HankelMatrix(O)
			L = cholesky(M)
		except ValueError:
			S2_reduced[i].remove(S2[j])
	if S2_reduced[i] == []:
		S1_reduced[i] = 0		
print("\n Reduced search space: \n",S1_reduced)

bounds = []
for i in range(1,len(S1_reduced)-1):
	if i==1 and S1_reduced[0]!=0:
		bounds.append(S1_reduced[0])
	if i==len(S1_reduced)-2 and S1_reduced[-1]!=0:
		bounds.append(S1_reduced[-1])
	if ((S1_reduced[i-1]==0) and (S1_reduced[i]!=0)): 
		bounds.append(S1_reduced[i])
	if ((S1_reduced[i+1]==0) and (S1_reduced[i]!=0)):
		bounds.append(S1_reduced[i])				
print("\n Bounds for accepted regions: \n",bounds)

results = []
for i in range(0,int(len(bounds)/2)):
	results.append(fdiv(fadd(bounds[2*i],bounds[2*i+1]),2))	
print("\n Accepted energies are: \n",results)

end = time.time()
print("\n",(end - start)/60)

#DoubleWell3.py
#This program obtains the graph of energies of the ground state and the first excited 
#state versus the coupling constant.

import time
import matplotlib.pyplot as plt
from mpmath import *
mp.dps = 50
plt.rcParams.update({
    "text.usetex": True,
    "font.family": "cmr"
})

start = time.time()

#making a grid
def MakingGrid(start,end,arr,step):
	x = start
	while x <= end:
		arr.append(x)
		x = fadd(x,step)
	return(arr)

#recursion for creating a moment sequence
def RecursionRelation(x,y,z,w,K): 
	x =[1,0,mpf(z),0,fdiv(fadd(fsub(fmul(2,mpf(y)),fdiv(1,fmul(16,mpf(w)))),mpf(z)),
	fmul(3,mpf(w))),0]
	for s in range(6,2*K-1):
		if s%2 == 0:
			x.append(fsum([fdiv(fprod([fsub(fmul(2,mpf(y)),fdiv(1,fmul(16,mpf(w)))),
			s-3,x[s-4]]),fmul(mpf(w),(s-1))),fdiv(fmul(s-2,x[s-2]),fprod([2,mpf(w),s-1])),
			fdiv(fprod([s-3,s-4,s-5,x[s-6]]),fprod([4,mpf(w),(s-1)]))]))
		else:
			x.append(0)
	return x
	
#creating a Hankel matrix from the moment sequence
def HankelMatrix(x):
	d = int((len(x) + 1)/2)
	H = zeros(d)
	for i in range(0,d):
		for j in range(0,d):
			if fmul(x[i+1],x[j+1]) != 0:
				H[i,j] = fdiv(x[i+j],fmul(x[i+1],x[j+1]))
			else:
				H[i,j] = x[i+j]
	return H

#creating a search space
N = 900
K = 18
g_min, g_max = 0.037, 0.050
g = g_min; g_step = 0.001
g_list = []; GS = []
while g<(g_max+g_step):
	g_list.append(g)
	E_min, E_max = 0, 1/(32*g)
	length1 = fsub(E_max,E_min)
	S1 = MakingGrid(E_min,E_max,[],fdiv(length1,N))
	x1 = sqrt((1/(2*g)-sqrt(8*S1[1]/g))/2)
	x2 = sqrt((1/(2*g)+sqrt(8*S1[1]/g))/2)	
	integrand1 = lambda x: 1/(sqrt(2*S1[1]-g*(x**2-1/(4*g))**2))
	integrand2 = lambda x: x**2/(sqrt(2*S1[1]-g*(x**2-1/(4*g))**2))
	T = 2*re(quad(integrand1,[x1,x2]))
	mom_max = 2*re(quad(integrand2,[x1,x2])/T); mom_min = 0
		
	length2 = fsub(mom_max,mom_min)
	S2 = MakingGrid(mom_min,mom_max,[],fdiv(length2,N))
	S1_reduced = MakingGrid(E_min,E_max,[],fdiv(length1,N))
	S2_reduced = []
	for i in range(0,len(S1)):
		S2_reduced.append(MakingGrid(mom_min,mom_max,[],fdiv(length2,N)))

	#reducing the search space
	for i in range(0,len(S1)):
		for j in range(0,len(S2)):
			O = RecursionRelation([],S1[i],S2[j],g,K)
			try:
				M = HankelMatrix(O)
				L = cholesky(M)
			except ValueError:
				S2_reduced[i].remove(S2[j])
		if S2_reduced[i] == []:
			S1_reduced[i] = 0	
	print(f"\n Reduced search space for g = {g}: \n",S1_reduced)

	bounds = []
	for i in range(1,len(S1_reduced)-1):
		if i==1 and S1_reduced[0]!=0:
			bounds.append(S1_reduced[0])
		if i==len(S1_reduced)-2 and S1_reduced[-1]!=0:
			bounds.append(S1_reduced[-1])
		if ((S1_reduced[i-1]==0) and (S1_reduced[i]!=0)): 
			bounds.append(S1_reduced[i])
		if ((S1_reduced[i+1]==0) and (S1_reduced[i]!=0)):
			bounds.append(S1_reduced[i])		
	print(f"\n Bounds for accepted regions for g = {g}: \n",bounds)

	results = []
	for i in range(0,int(len(bounds)/2)):
		results.append(fdiv(fadd(bounds[2*i],bounds[2*i+1]),2))
	
	print(f"\n Accepted energies for g = {g} are: \n",results)
	GS.append(results)
	g += g_step
print("\n Accepted energies for given g are: \n",GS)

#plotting
GS1 = []; GS2 = []
for i in range(len(GS)):
	if len(GS[i])==1:
		GS1.append(GS[i][0])
		GS2.append(GS[i][0])
	else:
		GS1.append(GS[i][0])
		GS2.append(GS[i][1])
		
plt.plot(g_list,GS1,color = "red",marker = "o",linestyle = "--")
plt.plot(g_list,GS2,color = "blue",marker = "o",linestyle = "--")
plt.legend(["Lower ground state","Higher ground state"], loc="lower left",fontsize=14)
plt.title("Ground state energy vs. $g$",fontsize=16)	
plt.xlabel("$g$",fontsize=14)
plt.ylabel("$E$",fontsize=14)
plt.grid(linestyle = '-', linewidth = 0.5)
plt.savefig("GroundState.pdf", dpi=300)

end = time.time()
print("\n",(end - start)/3600)

#DoubleWell4.py
#This program obtains the spectrum of the double-well using the subsequent 
#approach of the bootstrap method. Also, it tracks the convergence of the method.

import numpy as np
import scipy as sp
from scipy import odr
import matplotlib.pyplot as plt
from mpmath import *
mp.dps = 50
plt.rcParams.update({
    "text.usetex": True,
    "font.family": "cmr"
})

#making a grid
def MakingGrid(start,end,arr,step):
	x = start
	while x <= end:
		arr.append(x)
		x = fadd(x,step)
	return(arr)

#recursion for creating a moment sequence
def RecursionRelation(x,y,z,w,K): 
	x =[1,0,mpf(z),0,fdiv(fadd(fsub(fmul(2,mpf(y)),fdiv(1,fmul(16,mpf(w)))),mpf(z)),
	fmul(3,mpf(w))),0]
	for s in range(6,2*K-1):
		if s%2 == 0:
			x.append(fsum([fdiv(fprod([fsub(fmul(2,mpf(y)),fdiv(1,fmul(16,mpf(w)))),
			s-3,x[s-4]]),fmul(mpf(w),(s-1))),fdiv(fmul(s-2,x[s-2]),fprod([2,mpf(w),s-1])),
			fdiv(fprod([s-3,s-4,s-5,x[s-6]]),fprod([4,mpf(w),(s-1)]))]))
		else:
			x.append(0)
	return x
	
#creating a Hankel matrix from the moment sequence
def HankelMatrix(x):
	d = int((len(x) + 1)/2)
	H = zeros(d)
	for i in range(0,d):
		for j in range(0,d):
			if fmul(x[i+1],x[j+1]) != 0:
				H[i,j] = fdiv(x[i+j],fmul(x[i+1],x[j+1]))
			else:
				H[i,j] = x[i+j]
	return H
	
#defining a fitting function	
def fitting_function(p,x):
	a,b,c = p
	return a*np.exp(-b*(x-10)**c)

#creating a search space
g = 0.05
N = 300
K = 10; K_max = 18 
K_list = []
E_min, E_max = 0, 3
if E_max<=(1/(32*g)):
	x1 = sqrt((1/(2*g)-sqrt(8*E_max/g))/2)
	x2 = sqrt((1/(2*g)+sqrt(8*E_max/g))/2)
else:
	x1 = -sqrt((1/(2*g)+sqrt(8*E_max/g))/2)
	x2 = sqrt((1/(2*g)+sqrt(8*E_max/g))/2)	
integrand1 = lambda x: 1/(sqrt(2*E_max-g*(x**2-1/(4*g))**2))
integrand2 = lambda x: x**2/(sqrt(2*E_max-g*(x**2-1/(4*g))**2))
T = 2*re(quad(integrand1,[x1,x2]))
mom_max = 2*re(quad(integrand2,[x1,x2])/T); mom_min = 0
	
length1 = fsub(E_max,E_min)
length2 = fsub(mom_max,mom_min)
S1 = [MakingGrid(E_min,E_max,[],fdiv(length1,N))]
S2 = MakingGrid(mom_min,mom_max,[],fdiv(length2,N))
S1_reduced = [MakingGrid(E_min,E_max,[],fdiv(length1,N))]
S2_reduced = MakingGrid(mom_min,mom_max,[],fdiv(length2,N))
results = []
convergence = []

while K<=K_max:
	#reducing the search space
	for k in range(0,len(S1)):
		for i in range(0,len(S1[k])):
			for j in range(0,len(S2)):
				O = RecursionRelation([],S1[k][i],S2[j],g,K)
				try:
					M = HankelMatrix(O)
					L = cholesky(M)
				except ValueError:
					S2_reduced.remove(S2[j])
			if S2_reduced == []:
				S1_reduced[k][i] = 0
			S2_reduced = MakingGrid(mom_min,mom_max,[],fdiv(length2,N))	
	print(f"\n Reduced search space for K = {K}: \n",S1_reduced)
	
	bounds1 = []
	for j in range(0,len(S1_reduced)):
		for i in range(1,len(S1_reduced[j])-1):
			if i==1 and S1_reduced[j][0]!=0:
				bounds1.append(S1_reduced[j][0])
			if i==len(S1_reduced[j])-2 and S1_reduced[j][-1]!=0:
				bounds1.append(S1_reduced[j][-1])	
			if ((S1_reduced[j][i-1]==0) and (S1_reduced[j][i]!=0)): 
				bounds1.append(S1_reduced[j][i])
			if ((S1_reduced[j][i+1]==0) and (S1_reduced[j][i]!=0)):
				bounds1.append(S1_reduced[j][i])
			if ((S1_reduced[j][i-1]==0) and (S1_reduced[j][i]!=0) and 
			(S1_reduced[j][i+1]==0)):
				results.append([S1_reduced[j][i],fdiv(fsub(S1[j][i+1],S1[j][i-1]),2)])								
	print(f"\n Bounds for accepted regions for K = {K}: \n",bounds1)
	
	bounds2 = []
	for i in range(0,len(bounds1)-1,2):
		if bounds1[i] != bounds1[i+1]:
			bounds2.append(bounds1[i])
			bounds2.append(bounds1[i+1])
	print(f"\n Modified bounds for accepted regions for K = {K}: \n",bounds2)
	
	totwidth = 0
	if len(bounds2)>=2:
		for i in range(0,len(bounds2)-1,2):
			totwidth = fadd(totwidth,fsub(bounds2[i+1],bounds2[i]))
	convergence.append(totwidth)
		
	S1 = []; S1_reduced = []
	for i in range(0,int(len(bounds2)/2)):
		S1.append(MakingGrid(bounds2[2*i],bounds2[2*i+1],[],
		fdiv(fsub(bounds2[2*i+1],bounds2[2*i]),N)))
		S1_reduced.append(MakingGrid(bounds2[2*i],bounds2[2*i+1],[],
		fdiv(fsub(bounds2[2*i+1],bounds2[2*i]),N)))
	K_list.append(K)
	K += 1
	
for i in range(0,int(len(bounds2)/2)):
	results.append([fdiv(fadd(bounds2[2*i],bounds2[2*i+1]),2),
	fdiv(fsub(bounds2[2*i+1],bounds2[2*i]),2)])
results.sort()
print("\n Accepted energies are: \n",results)
print("\n Total width of accepted energies for given K: \n",convergence)

#fitting data for convergence
for i in range(0,len(convergence)):
	convergence[i] = float(str(convergence[i]))
	
odr_model = sp.odr.Model(fitting_function)
data = odr.Data(K_list,convergence)
fit = sp.odr.ODR(data,odr_model,beta0=[2,0.02,3])
out = fit.run()
beta = out.beta
std = out.sd_beta
print("\n Parameters of the fit: ",beta,
"\n Standart deviation of the parameters of the fit: ",std)

#plotting
plt.plot(K_list,convergence,color = "black",marker = "o",linestyle = "None")
plt.plot(np.linspace(K_list[0],K_list[-1],50),
beta[0]*np.exp(-beta[1]*(np.linspace(K_list[0],K_list[-1],50)-10)**beta[2]), 
color = "red", linestyle = "--", linewidth = 1)
plt.legend(["Data","Fitted curve $ae^{-b(K-10)^c}$"],loc = "upper right", fontsize=14)
plt.title("Convergence of total interval width",fontsize=16)
plt.xlabel("$K$",fontsize=14)
plt.ylabel("Total interval width",fontsize=14)
plt.grid(linestyle = '-', linewidth = 0.5)
plt.savefig("ConvergenceDoubleWell.pdf", dpi=300)

#results
#[[mpf('0.38500731825598968888888888888888888888888888888905774'),
#mpf('0.00018427105329339259259259259259259259259259259263156421')],
#[mpf('0.46004632296296296296296296296296296296296296296301074'),
#mpf('0.0007083585185185185185185185185185185185185185185187336')],
#[mpf('0.49711127185266958683127572016460905349794238683136109'),
#mpf('0.00019051459532440164609053497942386831275720164610727058')],
#[mpf('1.0643332315934011522633744855967078189300411522632037'),
#mpf('0.00078962488095824417009602194787379972565157750340581069')],
#[mpf('1.598485666831874865514403292181069958847736625514162'),
#mpf('0.0031721236683901576954732510288065843621399176953393079')],
#[mpf('1.8078714192592592592592592592592592592592592592590193'),
#mpf('0.0034735383703703703703703703703703703703703703703688972')], 
#[mpf('1.8287126494814814814814814814814814814814814814812327'),
#mpf('0.0034735383703703703703703703703703703703703703703688972')],
#[mpf('1.8460803413333333333333333333333333333333333333330772'),
#mpf('0.0034735383703703703703703703703703703703703703703688972')], 
#[mpf('1.884289263407407407407407407407407407407407407407135'),
#mpf('0.0034735383703703703703703703703703703703703703703688972')],
#[mpf('2.2435506794396049382716049382716049382716049382704013'),
#mpf('0.001984741229958847736625514403292181069958847736247743')],
#[mpf('2.2633389120625107753086419753086419753086419753078807'),
#mpf('0.0024279323063266502057613168724279835390946502058009403')],
#[mpf('2.2827222743983558452674897119341563786008230452668841'),
#mpf('0.0019216589364299193415637860082304526748971193417973891')], 
#[mpf('2.7841231355191784287517146776406035665294924554181078'),
#mpf('0.0000086764112673176406035665294924554183813443074175045279')],
#[mpf('2.8077900770460365432098765432098765432098765432098045'),
#mpf('0.00024579915354864197530864197530864197530864197520797301')],
#[mpf('2.8316262696529876543209876543209876543209876543217444'),
#mpf('0.00040405118663374485596707818930041152263374485617643594')],
#[mpf('2.855615735537332137174211248285322359396433470507773'),
#mpf('0.00051578517062453333333333333333333333333333333353064338')],
#[mpf('2.8797538475264077958847736625514403292181069958847435'),
#mpf('0.00058664109322808888888888888888888888888888888882271157')],
#[mpf('2.9040314383107259193415637860082304526748971193413168'),
#mpf('0.00063526258013235226337448559670781893004115226323443868')],
#[mpf('2.928447906570000197530864197530864197530864197531168'),
#mpf('0.00065352393358103703703703703703703703703703703707270418')],
#[mpf('2.9529981761743965322359396433470507544581618655697338'),
#mpf('0.00065364369199771742112482853223593964334705075471313141')],
#[mpf('2.9776929392531144647462277091906721536351165980794469'),
#mpf('0.00064564071304849821673525377229080932784636488322655824')]]

#convergence
#[mpf('2.3599999999999999999999999999999999999999999999999622'),
#mpf('2.2952000000000000000000000000000000000000000000000377'),
#mpf('1.9677040000000000000000000000000000000000000000000239'),
#mpf('1.3831647674074074074074074074074074074074074074132735'),
#mpf('1.0153060478222222222222222222222222222222222222217852'),
#mpf('0.78367976025053497942386831275720164609053497942630204'),
#mpf('0.21419701285336581618655692729766803840877914952360425'),
#mpf('0.064994164480024625514403292181069958847736625513805028'),
#mpf('0.030039781209487090041152263374485596707818930041464253')]

#fit
#[2.3669436  0.04273358 2.159403  ]
#[0.07796334 0.01635894 0.23888899]

#DoubleWell5.py
#This program obtains the graph of the ground–state energy splitting versus the 
#coupling constant, together with estimates of approximation methods.

import time
import matplotlib.pyplot as plt
from mpmath import *
mp.dps = 50
plt.rcParams.update({
    "text.usetex": True,
    "font.family": "cmr"
})

start = time.time()

#making a grid
def MakingGrid(start,end,arr,step):
	x = start
	while x <= end:
		arr.append(x)
		x = fadd(x,step)
	return(arr)

#recursion for creating a moment sequence
def RecursionRelation(x,y,z,w,K): 
	x =[1,0,mpf(z),0,fdiv(fadd(fsub(fmul(2,mpf(y)),fdiv(1,fmul(16,mpf(w)))),mpf(z)),
	fmul(3,mpf(w))),0]
	for s in range(6,2*K-1):
		if s%2 == 0:
			x.append(fsum([fdiv(fprod([fsub(fmul(2,mpf(y)),fdiv(1,fmul(16,mpf(w)))),
			s-3,x[s-4]]),fmul(mpf(w),(s-1))),fdiv(fmul(s-2,x[s-2]),fprod([2,mpf(w),s-1])),
			fdiv(fprod([s-3,s-4,s-5,x[s-6]]),fprod([4,mpf(w),(s-1)]))]))
		else:
			x.append(0)
	return x
	
#creating a Hankel matrix from the moment sequence
def HankelMatrix(x):
	d = int((len(x) + 1)/2)
	H = zeros(d)
	for i in range(0,d):
		for j in range(0,d):
			if fmul(x[i+1],x[j+1]) != 0:
				H[i,j] = fdiv(x[i+j],fmul(x[i+1],x[j+1]))
			else:
				H[i,j] = x[i+j]
	return H

#creating a search space
N = 900
K = 18
g_min, g_max = 0.037, 0.050
g = g_min; g_step = 0.001
g_list = []; GS = []; OneLoop = []; TwoLoop = []; ThreeLoop = []
while g<(g_max+g_step):
	g_list.append(g)
	OneLoop.append(2*exp(-1/(6*g))/(sqrt(pi*g)))
	TwoLoop.append(2*exp(-1/(6*g)-71*g/12)/(sqrt(pi*g)))
	ThreeLoop.append(2*exp(-1/(6*g))*(1-71*g/12-21.8713*g**2)/(sqrt(pi*g)))
	E_min, E_max = 0, 1/(32*g)
	length1 = fsub(E_max,E_min)
	S1 = MakingGrid(E_min,E_max,[],fdiv(length1,N))
	x1 = sqrt((1/(2*g)-sqrt(8*S1[1]/g))/2)
	x2 = sqrt((1/(2*g)+sqrt(8*S1[1]/g))/2)	
	integrand1 = lambda x: 1/(sqrt(2*S1[1]-g*(x**2-1/(4*g))**2))
	integrand2 = lambda x: x**2/(sqrt(2*S1[1]-g*(x**2-1/(4*g))**2))
	T = 2*re(quad(integrand1,[x1,x2]))
	mom_max = 2*re(quad(integrand2,[x1,x2])/T); mom_min = 0
	
	length2 = fsub(mom_max,mom_min)
	S2 = MakingGrid(mom_min,mom_max,[],fdiv(length2,N))
	S1_reduced = MakingGrid(E_min,E_max,[],fdiv(length1,N))
	S2_reduced = []
	for i in range(0,len(S1)):
		S2_reduced.append(MakingGrid(mom_min,mom_max,[],fdiv(length2,N)))

	#reducing the search space
	for i in range(0,len(S1)):
		for j in range(0,len(S2)):
			O = RecursionRelation([],S1[i],S2[j],g,K)
			try:
				M = HankelMatrix(O)
				L = cholesky(M)
			except ValueError:
				S2_reduced[i].remove(S2[j])
		if S2_reduced[i] == []:
			S1_reduced[i] = 0	
	print(f"\n Reduced search space for g = {g}: \n",S1_reduced)

	bounds = []
	for i in range(1,len(S1_reduced)-1):
		if i==1 and S1_reduced[0]!=0:
			bounds.append(S1_reduced[0])
		if i==len(S1_reduced)-2 and S1_reduced[-1]!=0:
			bounds.append(S1_reduced[-1])
		if ((S1_reduced[i-1]==0) and (S1_reduced[i]!=0)): 
			bounds.append(S1_reduced[i])
		if ((S1_reduced[i+1]==0) and (S1_reduced[i]!=0)):
			bounds.append(S1_reduced[i])		
	print(f"\n Bounds for accepted regions for g = {g}: \n",bounds)

	results = []
	for i in range(0,int(len(bounds)/2)):
		results.append(fdiv(fadd(bounds[2*i],bounds[2*i+1]),2))
	
	print(f"\n Accepted energies for g = {g} are: \n",results)
	GS.append(results)
	g += g_step
print("\n Accepted energies for given g are: \n",GS)

#plotting
GS1 = []; GS2 = []; DeltaGS = []
for i in range(len(GS)):
	if len(GS[i])==1:
		GS1.append(GS[i][0])
		GS2.append(GS[i][0])
	else:
		GS1.append(GS[i][0])
		GS2.append(GS[i][1])
for i in range(len(GS1)):
	DeltaGS.append(GS2[i]-GS1[i])
				
plt.plot(g_list,DeltaGS,color = "black",marker = "o",linestyle = "None")
plt.plot(g_list,OneLoop,color = "red",marker = "None",linestyle = "--",linewidth = 1)
plt.plot(g_list,TwoLoop,color = "blue",marker = "None",linestyle = "--",linewidth = 1)
plt.plot(g_list,ThreeLoop,color = "green",marker = "None",linestyle = "--",linewidth = 1)
plt.legend(["GS splitting","One loop approximation","Two loop approximation",
"Three loop approximation"], fontsize=14)
plt.title("Ground state energy splitting vs. $g$",fontsize=16)	
plt.xlabel("$g$",fontsize=14)
plt.ylabel("$\Delta E$",fontsize=14)
plt.grid(linestyle = '-', linewidth = 0.5)
plt.savefig("GroundStateSplitting.pdf", dpi=300)

end = time.time()
print("\n",(end - start)/3600)

#DoubleWell6.py
#This program obtains the picture of the reduced search space for the double-well
#potential using the subsequent approach of the bootstrap method (we are using the data 
#from the DoubleWell4.py.).

import numpy as np
import scipy as sp
from scipy import odr
import matplotlib.pyplot as plt
from mpmath import *
mp.dps = 50
plt.rcParams.update({
    "text.usetex": True,
    "font.family": "cmr"
})

#making a grid
def MakingGrid(start,end,arr,step):
	x = start
	while x <= end:
		arr.append(x)
		x = fadd(x,step)
	return(arr)

#recursion for creating a moment sequence
def RecursionRelation(x,y,z,w,K): 
	x =[1,0,mpf(z),0,fdiv(fadd(fsub(fmul(2,mpf(y)),fdiv(1,fmul(16,mpf(w)))),mpf(z)),
	fmul(3,mpf(w))),0]
	for s in range(6,2*K-1):
		if s%2 == 0:
			x.append(fsum([fdiv(fprod([fsub(fmul(2,mpf(y)),fdiv(1,fmul(16,mpf(w)))),
			s-3,x[s-4]]),fmul(mpf(w),(s-1))),fdiv(fmul(s-2,x[s-2]),fprod([2,mpf(w),s-1])),
			fdiv(fprod([s-3,s-4,s-5,x[s-6]]),fprod([4,mpf(w),(s-1)]))]))
		else:
			x.append(0)
	return x
	
#creating a Hankel matrix from the moment sequence
def HankelMatrix(x):
	d = int((len(x) + 1)/2)
	H = zeros(d)
	for i in range(0,d):
		for j in range(0,d):
			if fmul(x[i+1],x[j+1]) != 0:
				H[i,j] = fdiv(x[i+j],fmul(x[i+1],x[j+1]))
			else:
				H[i,j] = x[i+j]
	return H

#data from the DoubleWell4.py
bounds = [[mpf('0.28999999999999999999999999999999999999999999999999735'), mpf('0.56999999999999999999999999999999999999999999999999353'), mpf('0.91999999999999999999999999999999999999999999999998792'), mpf('2.9999999999999999999999999999999999999999999999999519')],
[mpf('0.34133333333333333333333333333333333333333333333334554'), mpf('0.5634666666666666666666666666666666666666666666667425'), mpf('0.91999999999999999999999999999999999999999999999998792'), mpf('2.993066666666666666666666666666666666666666666666628')],
[mpf('0.35095911111111111111111111111111111111111111111112576'), mpf('0.56346666666666666666666666666666666666666666666673582'), mpf('1.0236533333333333333333333333333333333333333333333211'), mpf('2.6268248888888888888888888888888888888888888888888822'), mpf('2.6406453333333333333333333333333333333333333333333234'), mpf('2.6959271111111111111111111111111111111111111111110883'), mpf('2.7097475555555555555555555555555555555555555555555295'), mpf('2.7166577777777777777777777777777777777777777777777501'), mpf('2.7304782222222222222222222222222222222222222222221914'), mpf('2.7442986666666666666666666666666666666666666666666326'), mpf('2.7581191111111111111111111111111111111111111111110738'), mpf('2.7650293333333333333333333333333333333333333333332944'), mpf('2.7788497777777777777777777777777777777777777777777356'), mpf('2.7926702222222222222222222222222222222222222222221769'), mpf('2.8064906666666666666666666666666666666666666666666181'), mpf('2.8134008888888888888888888888888888888888888888888387'), mpf('2.8272213333333333333333333333333333333333333333332799'), mpf('2.8341315555555555555555555555555555555555555555555005'), mpf('2.8548622222222222222222222222222222222222222222221623'), mpf('2.861772444444444444444444444444444444444444444444383'), mpf('2.8755928888888888888888888888888888888888888888888242'), mpf('2.8825031111111111111111111111111111111111111111110448'), mpf('2.9032337777777777777777777777777777777777777777777066'), mpf('2.9101439999999999999999999999999999999999999999999272'), mpf('2.9239644444444444444444444444444444444444444444443684'), mpf('2.9308746666666666666666666666666666666666666666665891'), mpf('2.9516053333333333333333333333333333333333333333332509'), mpf('2.9585155555555555555555555555555555555555555555554715'), mpf('2.9723359999999999999999999999999999999999999999999127'), mpf('2.9792462222222222222222222222222222222222222222221333')],
[mpf('0.35095911111111111111111111111111111111111111111112576'), mpf('0.41258630222222222222222222222222222222222222222225558'), mpf('0.46004632296296296296296296296296296296296296296301074'), mpf('0.46004632296296296296296296296296296296296296296301074'), mpf('0.4614630400000000000000000000000000000000000000000482'), mpf('0.50396455111111111111111111111111111111111111111116888'), mpf('1.0236533333333333333333333333333333333333333333333211'), mpf('1.1038119111111111111111111111111111111111111111110834'), mpf('1.4778852740740740740740740740740740740740740740739741'), mpf('2.5199467851851851851851851851851851851851851851848838'), mpf('2.5306345955555555555555555555555555555555555555552521'), mpf('2.5413224059259259259259259259259259259259259259256204'), mpf('2.5520102162962962962962962962962962962962962962959887'), mpf('2.562698026666666666666666666666666666666666666666357'), mpf('2.5733858370370370370370370370370370370370370370367253'), mpf('2.5840736474074074074074074074074074074074074074070936'), mpf('2.5947614577777777777777777777777777777777777777774619'), mpf('2.6054492681481481481481481481481481481481481481478302'), mpf('2.6268248888888888888888888888888888888888888888885668'), mpf('2.6214809837037037037037037037037037037037037037033827'), mpf('2.6406453333333333333333333333333333333333333333333234'), mpf('2.654834322962962962962962962962962962962962962962851'), mpf('2.6634951348148148148148148148148148148148148148146405'), mpf('2.6774998518518518518518518518518518518518518518515768'), mpf('2.6867134814814814814814814814814814814814814814811401'), mpf('2.6959271111111111111111111111111111111111111111107034'), mpf('2.7100700325925925925925925925925925925925925925926017'), mpf('2.7166347437037037037037037037037037037037037037044273'), mpf('2.7336108562962962962962962962962962962962962962962429'), mpf('2.7442986666666666666666666666666666666666666666665364'), mpf('2.7581191111111111111111111111111111111111111111110738'), mpf('2.7650062992592592592592592592592592592592592592599716'), mpf('2.781014980740740740740740740740740740740740740740683'), mpf('2.7926702222222222222222222222222222222222222222220806'), mpf('2.8064906666666666666666666666666666666666666666666181'), mpf('2.8133778548148148148148148148148148148148148148155158'), mpf('2.8289258548148148148148148148148148148148148148149469'), mpf('2.8341085214814814814814814814814814814814814814821777'), mpf('2.8548622222222222222222222222222222222222222222221623'), mpf('2.8617494103703703703703703703703703703703703703710601'), mpf('2.8773204444444444444444444444444444444444444444445678'), mpf('2.8824800770370370370370370370370370370370370370377219'), mpf('2.9032337777777777777777777777777777777777777777777066'), mpf('2.9101209659259259259259259259259259259259259259266044'), mpf('2.9261987496296296296296296296296296296296296296297968'), mpf('2.9308516325925925925925925925925925925925925925932662'), mpf('2.9516053333333333333333333333333333333333333333332509'), mpf('2.9584925214814814814814814814814814814814814814821486'), mpf('2.9755607703703703703703703703703703703703703703706341'), mpf('2.9792231881481481481481481481481481481481481481488105')],
[mpf('0.38136185872592592592592592592592592592592592592597651'), mpf('0.41114833442962962962962962962962962962962962962971542'), mpf('0.41176460634074074074074074074074074074074074074082725'), mpf('0.41238087825185185185185185185185185185185185185193909'), mpf('0.4614630400000000000000000000000000000000000000000482'), mpf('0.46245474192592592592592592592592592592592592592597256'), mpf('0.46273808533333333333333333333333333333333333333337952'), mpf('0.50382287940740740740740740740740740740740740740740593'), mpf('1.0236533333333333333333333333333333333333333333333211'), mpf('1.0987352011851851851851851851851851851851851851850834'), mpf('1.4778852740740740740740740740740740740740740740739741'), mpf('1.7314535751111111111111111111111111111111111111109036'), mpf('1.7418741902222222222222222222222222222222222222220103'), mpf('1.752294805333333333333333333333333333333333333333117'), mpf('1.7661889588148148148148148148148148148148148148145925'), mpf('1.7696624971851851851851851851851851851851851851849614'), mpf('1.7870301890370370370370370370370370370370370370368059'), mpf('1.7905037274074074074074074074074074074074074074071748'), mpf('1.8078714192592592592592592592592592592592592592590193'), mpf('1.8078714192592592592592592592592592592592592592590193'), mpf('1.8287126494814814814814814814814814814814814814812327'), mpf('1.8287126494814814814814814814814814814814814814812327'), mpf('1.8460803413333333333333333333333333333333333333330772'), mpf('1.8460803413333333333333333333333333333333333333330772'), mpf('1.8634480331851851851851851851851851851851851851849217'), mpf('1.8669215715555555555555555555555555555555555555552906'), mpf('1.884289263407407407407407407407407407407407407407135'), mpf('1.884289263407407407407407407407407407407407407407135'), mpf('1.9016569552592592592592592592592592592592592592589795'), mpf('1.9051304936296296296296296296296296296296296296293484'), mpf('1.919024647111111111111111111111111111111111111110824'), mpf('1.9259717238518518518518518518518518518518518518515618'), mpf('1.9398658773333333333333333333333333333333333333330374'), mpf('1.9468129540740740740740740740740740740740740740737752'), mpf('1.9572335691851851851851851851851851851851851851848819'), mpf('1.9676541842962962962962962962962962962962962962959886'), mpf('1.9780747994074074074074074074074074074074074074070953'), mpf('1.988495414518518518518518518518518518518518518518202'), mpf('1.9954424912592592592592592592592592592592592592589398'), mpf('2.0093366447407407407407407407407407407407407407404153'), mpf('2.0162837214814814814814814814814814814814814814811531'), mpf('2.3289021748148148148148148148148148148148148148143539'), mpf('2.3358492515555555555555555555555555555555555555550917'), mpf('2.3497434050370370370370370370370370370370370370365673'), mpf('2.3566904817777777777777777777777777777777777777773051'), mpf('2.3705846352592592592592592592592592592592592592587807'), mpf('2.3775317119999999999999999999999999999999999999995184'), mpf('2.391425865481481481481481481481481481481481481480994'), mpf('2.3983729422222222222222222222222222222222222222217318'), mpf('2.4122670957037037037037037037037037037037037037032074'), mpf('2.4192141724444444444444444444444444444444444444439452'), mpf('2.4296347875555555555555555555555555555555555555550519'), mpf('2.4400554026666666666666666666666666666666666666661586'), mpf('2.4504760177777777777777777777777777777777777777772653'), mpf('2.4643701712592592592592592592592592592592592592587409'), mpf('2.4713172479999999999999999999999999999999999999994787'), mpf('2.4852114014814814814814814814814814814814814814809543'), mpf('2.495632016592592592592592592592592592592592592592061'), mpf('2.5060526317037037037037037037037037037037037037031676'), mpf('2.5164732468148148148148148148148148148148148148142743'), mpf('2.5306345955555555555555555555555555555555555555552521'), mpf('2.5388642095407407407407407407407407407407407407398512'), mpf('2.5520102162962962962962962962962962962962962962959887'), mpf('2.5606317166617283950617283950617283950617283950608068'), mpf('2.5733858370370370370370370370370370370370370370367253'), mpf('2.5825417279209876543209876543209876543209876543200239'), mpf('2.5947614577777777777777777777777777777777777777774619'), mpf('2.6046654953876543209876543209876543209876543209866331'), mpf('2.6406453333333333333333333333333333333333333333333234'), mpf('2.6494425069037037037037037037037037037037037037038559'), mpf('2.6634951348148148148148148148148148148148148148146405'), mpf('2.6721313769876543209876543209876543209876543209870355'), mpf('2.6867134814814814814814814814814814814814814814811401'), mpf('2.6949443239506172839506172839506172839506172839497391'), mpf('2.7100700325925925925925925925925925925925925925926017'), mpf('2.7166347437037037037037037037037037037037037037040264'), mpf('2.7336108562962962962962962962962962962962962962962429'), mpf('2.741056697520987654320987654320987654320987654320404'), mpf('2.7581191111111111111111111111111111111111111111110738'), mpf('2.7643405377382716049382716049382716049382716049389839'), mpf('2.781014980740740740740740740740740740740740740740683'), mpf('2.7877750207999999999999999999999999999999999999994689'), mpf('2.8064906666666666666666666666666666666666666666666181'), mpf('2.8113805702518518518518518518518518518518518518523924'), mpf('2.8289258548148148148148148148148148148148148148149469'), mpf('2.8341085214814814814814814814814814814814814814817767'), mpf('2.8548622222222222222222222222222222222222222222221623'), mpf('2.8589945351111111111111111111111111111111111111115491'), mpf('2.8773204444444444444444444444444444444444444444445678'), mpf('2.8824800770370370370370370370370370370370370370374012'), mpf('2.9032337777777777777777777777777777777777777777777066'), mpf('2.907159475022222222222222222222222222222222222222624'), mpf('2.9261987496296296296296296296296296296296296296297968'), mpf('2.9308516325925925925925925925925925925925925925931058'), mpf('2.9516053333333333333333333333333333333333333333332509'), mpf('2.9558524326913580246913580246913580246913580246917873'), mpf('2.9755607703703703703703703703703703703703703703706341'), mpf('2.9792109800888888888888888888888888888888888888896298')],
[mpf('0.38136185872592592592592592592592592592592592592597651'), mpf('0.38851061289481481481481481481481481481481481481486823'), mpf('0.49451032608395061728395061728395061728395061728398268'), mpf('0.5036859300938271604938271604938271604938271604938371'), mpf('1.0549374449382716049382716049382716049382716049382031'), mpf('1.0987352011851851851851851851851851851851851851850379'), mpf('1.564098496426666666666666666666666666666666666666561'), mpf('1.730608347440987654320987654320987654320987654320871'), mpf('1.7418741902222222222222222222222222222222222222220103'), mpf('1.7522948053333333333333333333333333333333333333327802'), mpf('1.7661889588148148148148148148148148148148148148145925'), mpf('1.7696509187239506172839506172839506172839506172838814'), mpf('1.7870301890370370370370370370370370370370370370368059'), mpf('1.7904921489461728395061728395061728395061728395060948'), mpf('1.8634480331851851851851851851851851851851851851849217'), mpf('1.8669099930943209876543209876543209876543209876542105'), mpf('1.9016569552592592592592592592592592592592592592589795'), mpf('1.9051189151683950617283950617283950617283950617282684'), mpf('1.919024647111111111111111111111111111111111111110824'), mpf('1.9259485669293827160493827160493827160493827160494017'), mpf('1.9398658773333333333333333333333333333333333333330374'), mpf('1.9467666402291358024691358024691358024691358024691449'), mpf('1.9572335691851851851851851851851851851851851851848819'), mpf('1.9668900058548148148148148148148148148148148148141953'), mpf('1.9780747994074074074074074074074074074074074074070953'), mpf('1.9870365284029629629629629629629629629629629629623574'), mpf('1.9954424912592592592592592592592592592592592592589398'), mpf('2.0072062078735802469135802469135802469135802469135188'), mpf('2.0162837214814814814814814814814814814814814814811531'), mpf('2.0267043365925925925925925925925925925925925925922438'), mpf('2.0329567056592592592592592592592592592592592592588982'), mpf('2.0475455668148148148148148148148148148148148148144251'), mpf('2.0517138128592592592592592592592592592592592592588614'), mpf('2.0673447355259259259259259259259259259259259259254973'), mpf('2.0715129815703703703703703703703703703703703703699336'), mpf('2.0871439042370370370370370370370370370370370370365696'), mpf('2.0913121502814814814814814814814814814814814814810059'), mpf('2.1069430729481481481481481481481481481481481481476418'), mpf('2.110069257481481481481481481481481481481481481480969'), mpf('2.1277843031703703703703703703703703703703703703698231'), mpf('2.1298684261925925925925925925925925925925925925920413'), mpf('2.1475834718814814814814814814814814814814814814808954'), mpf('2.1496675949037037037037037037037037037037037037031135'), mpf('2.1673826405925925925925925925925925925925925925919676'), mpf('2.1705088251259259259259259259259259259259259259252948'), mpf('2.1882238708148148148148148148148148148148148148141489'), mpf('2.1903079938370370370370370370370370370370370370363671'), mpf('2.2080230395259259259259259259259259259259259259252212'), mpf('2.2101071625481481481481481481481481481481481481474393'), mpf('2.2278222082370370370370370370370370370370370370362934'), mpf('2.2309483927703703703703703703703703703703703703696206'), mpf('2.2486634384592592592592592592592592592592592592584747'), mpf('2.2507475614814814814814814814814814814814814814806929'), mpf('2.268462607170370370370370370370370370370370370369547'), mpf('2.2715887917037037037037037037037037037037037037028742'), mpf('2.2893038373925925925925925925925925925925925925917283'), mpf('2.2924300219259259259259259259259259259259259259250555'), mpf('2.3091030061037037037037037037037037037037037037028005'), mpf('2.3132712521481481481481481481481481481481481481472368'), mpf('2.3289021748148148148148148148148148148148148148138728'), mpf('2.3358492515555555555555555555555555555555555555550917'), mpf('2.349697091192098765432098765432098765432098765432247'), mpf('2.3566904817777777777777777777777777777777777777773051'), mpf('2.3705383214143209876543209876543209876543209876544604'), mpf('2.3775317119999999999999999999999999999999999999995184'), mpf('2.3913795516365432098765432098765432098765432098766738'), mpf('2.3983729422222222222222222222222222222222222222217318'), mpf('2.4122207818587654320987654320987654320987654320988872'), mpf('2.4192141724444444444444444444444444444444444444439452'), mpf('2.4296000521718518518518518518518518518518518518518117'), mpf('2.4400554026666666666666666666666666666666666666661586'), mpf('2.4504412823940740740740740740740740740740740740740251'), mpf('2.4643701712592592592592592592592592592592592592587409'), mpf('2.4712940910775308641975308641975308641975308641973186'), mpf('2.4852114014814814814814814814814814814814814814809543'), mpf('2.4955972812088888888888888888888888888888888888888208'), mpf('2.5060526317037037037037037037037037037037037037031676'), mpf('2.5164385114311111111111111111111111111111111111110342'), mpf('2.5306345955555555555555555555555555555555555555552521'), mpf('2.5388642095407407407407407407407407407407407407396908'), mpf('2.5520102162962962962962962962962962962962962962959887'), mpf('2.5606029783271769547325102880658436213991769547318039'), mpf('2.5733858370370370370370370370370370370370370370367253'), mpf('2.5825417279209876543209876543209876543209876543194359'), mpf('2.5947614577777777777777777777777777777777777777774619'), mpf('2.6046654953876543209876543209876543209876543209858848'), mpf('2.6406453333333333333333333333333333333333333333333234'), mpf('2.6494131829918024691358024691358024691358024691362098'), mpf('2.6634951348148148148148148148148148148148148148146405'), mpf('2.6721313769876543209876543209876543209876543209869018'), mpf('2.6868232260477366255144032921810699588477366255140553'), mpf('2.6949168878090534979423868312757201646090534979415531'), mpf('2.7102669739259259259259259259259259259259259259259444'), mpf('2.7166347437037037037037037037037037037037037037040264'), mpf('2.7338342315330370370370370370370370370370370370369886'), mpf('2.7410318780502386831275720164609053497942386831276827'), mpf('2.7581191111111111111111111111111111111111111111110738'), mpf('2.7643405377382716049382716049382716049382716049389678'), mpf('2.7813529827437037037037037037037037037037037037035838'), mpf('2.7877750207999999999999999999999999999999999999986991'), mpf('2.8064906666666666666666666666666666666666666666666181'), mpf('2.8113805702518518518518518518518518518518518518523443'), mpf('2.8293577437037037037037037037037037037037037037038494'), mpf('2.8341085214814814814814814814814814814814814814817767'), mpf('2.8548622222222222222222222222222222222222222222221623'), mpf('2.8589807607348148148148148148148148148148148148155708'), mpf('2.8778536064790123456790123456790123456790123456791605'), mpf('2.8824800770370370370370370370370370370370370370374012'), mpf('2.9032337777777777777777777777777777777777777777777066'), mpf('2.907159475022222222222222222222222222222222222222608'), mpf('2.9268346436345679012345679012345679012345679012347824'), mpf('2.9308516325925925925925925925925925925925925925931058'), mpf('2.9516053333333333333333333333333333333333333333332509'), mpf('2.955852432691358024691358024691358024691358024691226'), mpf('2.976302979679802469135802469135802469135802469136061'), mpf('2.9792109800888888888888888888888888888888888888891273')],
[mpf('0.38312521808758518518518518518518518518518518518524042'), mpf('0.38538899024106666666666666666666666666666666666672788'), mpf('0.49451032608395061728395061728395061728395061728398268'), mpf('0.49790529956760493827160493827160493827160493827162437'), mpf('1.0549374449382716049382716049382716049382716049382031'), mpf('1.065156921395884773662551440329218106995884773662371'), mpf('1.564098496426666666666666666666666666666666666666561'), mpf('1.5652085621000954732510288065843621399176954732509211'), mpf('1.5718689561406683127572016460905349794238683127570815'), mpf('1.6251521084652510288065843621399176954732510288063649'), mpf('2.1848580121339259259259259259259259259259259259252797'), mpf('2.1881648206625185185185185185185185185185185185178687'), mpf('2.2030037765807407407407407407407407407407407407400574'), mpf('2.207963989373629629629629629629629629629629629628941'), mpf('2.2218581428551111111111111111111111111111111111103899'), mpf('2.2277631580847407407407407407407407407407407407400132'), mpf('2.2412231192699259259259259259259259259259259259251653'), mpf('2.2486043883069629629629629629629629629629629629621945'), mpf('2.2608451375241481481481481481481481481481481481473489'), mpf('2.2684035570180740740740740740740740740740740740732668'), mpf('2.2808006154619259259259259259259259259259259259250867'), mpf('2.2892447872402962962962962962962962962962962962954481'), mpf('2.3009888204705185185185185185185185185185185185180264'), mpf('2.3090474294897777777777777777777777777777777777776419'), mpf('2.3213993319348148148148148148148148148148148148138458'), mpf('2.3289021748148148148148148148148148148148148148137926'), mpf('2.3419884604610897119341563786008230452674897119342957'), mpf('2.3496509317266436213991769547325102880658436214000691'), mpf('2.3627835312178567901234567901234567901234567901235827'), mpf('2.3704921619488658436213991769547325102880658436222824'), mpf('2.3837632398364444444444444444444444444444444444445751'), mpf('2.3913333921710880658436213991769547325102880658444958'), mpf('2.4048814268513975308641975308641975308641975308643466'), mpf('2.4121746223933102880658436213991769547325102880667092'), mpf('2.4261727118618074074074074074074074074074074074070546'), mpf('2.4296000521718518518518518518518518518518518518515712'), mpf('2.4476370948676740740740740740740740740740740740737255'), mpf('2.4504412823940740740740740740740740740740740740737846'), mpf('2.4692630745975045267489711934156378600823045267489335'), mpf('2.4712710113448032921810699588477366255144032921812296'), mpf('2.4910274941288296296296296296296296296296296296292248'), mpf('2.4954241832134320987654320987654320987654320987651198'), mpf('2.5129419319228839506172839506172839506172839506168929'), mpf('2.5164385114311111111111111111111111111111111111107936'), mpf('2.5349962909677037037037037037037037037037037037030046'), mpf('2.5380138160956049382716049382716049382716049382706322'), mpf('2.55719451605492762688614540466392318244170096021845'), mpf('2.5595432043433683401920438957475994513031550068573551'), mpf('2.5795508035655637860082304526748971193415637860070838'), mpf('2.5812599031972345679012345679012345679012345678998564'), mpf('2.6020574321503868312757201646090534979423868312743334'), mpf('2.6031138628287736625514403292181069958847736625498985'), mpf('2.6698283790748971193415637860082304526748971193409655'), mpf('2.6702026162357201646090534979423868312757201646084301'), mpf('2.6923269160454320987654320987654320987654320987644812'), mpf('2.6932711765842524005486968449931412894375857338809758'), mpf('2.7150215753599999999999999999999999999999999999996709'), mpf('2.7164224847111111111111111111111111111111111111106795'), mpf('2.7378889057377272976680384087791495198902606310016539'), mpf('2.7397363016771423868312757201646090534979423868317085'), mpf('2.7609187530933333333333333333333333333333333333336261'), mpf('2.763158466679111111111111111111111111111111111111668'), mpf('2.7841144591079111111111111111111111111111111111106903'), mpf('2.7867260879174716049382716049382716049382716049375661'), mpf('2.8074686473837037037037037037037037037037037037037633'), mpf('2.8104025895348148148148148148148148148148148148151991'), mpf('2.8309888440740740740740740740740740740740740740745005'), mpf('2.8340926855555555555555555555555555555555555555565163'), mpf('2.8548622222222222222222222222222222222222222222221623'), mpf('2.8581845099557135802469135802469135802469135802472359'), mpf('2.8785167339256625514403292181069958847736625514404387'), mpf('2.8822641750776625514403292181069958847736625514402196'), mpf('2.9032337777777777777777777777777777777777777777777066'), mpf('2.9064659351757037037037037037037037037037037037040088'), mpf('2.9268346436345679012345679012345679012345679012347824'), mpf('2.9308114627030123456790123456790123456790123456795967'), mpf('2.9516053333333333333333333333333333333333333333332509'), mpf('2.9552719957790946502057613168724279835390946502056361'), mpf('2.976302979679802469135802469135802469135802469136061'), mpf('2.9792109800888888888888888888888888888888888888887157')],
[mpf('0.38482304720269629629629629629629629629629629629642618'), mpf('0.38538144433388839506172839506172839506172839506188282'), mpf('0.49517800420240263374485596707818930041152263374489816'), mpf('0.49520063735896032921810699588477366255144032921814953'), mpf('0.49559671759872000000000000000000000000000000000004855'), mpf('0.49642282781307588477366255144032921810699588477372363'), mpf('0.49692075725734518518518518518518518518518518518525382'), mpf('0.49736210381022024691358024691358024691358024691365558'), mpf('1.0588889758352153635116598079561042524005486968449847'), mpf('1.0651228564743593964334705075445816186556927297668902'), mpf('1.564098496426666666666666666666666666666666666666561'), mpf('1.5652085621000954732510288065843621399176954732507233'), mpf('1.5761316083266349300411522633744855967078189300410344'), mpf('1.5850121337140653827160493827160493827160493827159364'), mpf('1.5953135431634847078189300411522633744855967078188227'), mpf('1.6020627424579318518518518518518518518518518518517481'), mpf('2.1848580121339259259259259259259259259259259259252797'), mpf('2.1881317525772325925925925925925925925925925925919005'), mpf('2.2030037765807407407407407407407407407407407407400574'), mpf('2.2077986489471999999999999999999999999999999999992495'), mpf('2.2218581428551111111111111111111111111111111111103899'), mpf('2.2272907568663703703703703703703703703703703703691516'), mpf('2.2412231192699259259259259259259259259259259259251653'), mpf('2.2466360498970864197530864197530864197530864197518301'), mpf('2.2608451375241481481481481481481481481481481481473489'), mpf('2.2657833049268464197530864197530864197530864197519187'), mpf('2.2808006154619259259259259259259259259259259259250867'), mpf('2.2846567872407150617283950617283950617283950617275542'), mpf('2.7841144591079111111111111111111111111111111111106903'), mpf('2.7841318699666415144032921810699588477366255144028675'), mpf('2.8074686473837037037037037037037037037037037037037633'), mpf('2.8080358761995851851851851851851851851851851851851835'), mpf('2.8309888440740740740740740740740740740740740740745005'), mpf('2.8320338040395061728395061728395061728395061728399181'), mpf('2.8548622222222222222222222222222222222222222222221623'), mpf('2.8561357658533939094650205761316872427983539094652875'), mpf('2.8785167339256625514403292181069958847736625514404387'), mpf('2.8803404886196358847736625514403292181069958847736732'), mpf('2.9032337777777777777777777777777777777777777777777066'), mpf('2.9046667008908582716049382716049382716049382716051927'), mpf('2.9268478996981293827160493827160493827160493827162658'), mpf('2.9291014305035812345679012345679012345679012345684545'), mpf('2.9516053333333333333333333333333333333333333333332509'), mpf('2.9536586643029596707818930041152263374485596707822266'), mpf('2.9763417530185902880658436213991769547325102880660964'), mpf('2.9783385799661629629629629629629629629629629629629193')],
[mpf('0.38482304720269629629629629629629629629629629629642618'), mpf('0.3851915893092830814814814814814814814814814814816893'), mpf('0.49692075725734518518518518518518518518518518518525382'), mpf('0.49730178644799398847736625514403292181069958847746836'), mpf('1.0635436067124429080932784636488340192043895747597979'), mpf('1.0651228564743593964334705075445816186556927297666096'), mpf('1.5953135431634847078189300411522633744855967078188227'), mpf('1.6016577905002650232098765432098765432098765432095013'), mpf('2.2415659382096460905349794238683127572016460905341536'), mpf('2.2455354206695637860082304526748971193415637860066491'), mpf('2.2609109797561841251028806584362139917695473251020824'), mpf('2.2657668443688374255144032921810699588477366255136843'), mpf('2.2808006154619259259259259259259259259259259259250867'), mpf('2.2846439333347857646090534979423868312757201646086814'), mpf('2.7841144591079111111111111111111111111111111111106903'), mpf('2.7841318119304457463923182441700960219478737997255253'), mpf('2.8075442778924879012345679012345679012345679012345965'), mpf('2.8080358761995851851851851851851851851851851851850125'), mpf('2.831222218466353909465020576131687242798353909465568'), mpf('2.8320303208396213991769547325102880658436213991779208'), mpf('2.8550999503667076038408779149519890260631001371742423'), mpf('2.8561315207079566705075445816186556927297668038413036'), mpf('2.8791672064331797069958847736625514403292181069959208'), mpf('2.8803404886196358847736625514403292181069958847735663'), mpf('2.9033961757305935670781893004115226337448559670780823'), mpf('2.9046667008908582716049382716049382716049382716045512'), mpf('2.9277943826364191604938271604938271604938271604940953'), mpf('2.9291014305035812345679012345679012345679012345682407'), mpf('2.9523445324823988148148148148148148148148148148150207'), mpf('2.953651819866394249657064471879286694101508916324447'), mpf('2.9770472985400659665294924554183813443072702331962203'), mpf('2.9783385799661629629629629629629629629629629629626734')]]

#creating a search space and classical limit
g = 0.05
N = 300
K_min = 10; K_max = 18
K = K_min
E_min, E_max = 0, 3
length1 = fsub(E_max,E_min)
S1 = MakingGrid(E_min,E_max,[],fdiv(length1,N))
mom = []
for i in range(1,len(S1)):
	if S1[i]<=(1/(32*g)):
		x1 = sqrt((1/(2*g)-sqrt(8*S1[i]/g))/2)
		x2 = sqrt((1/(2*g)+sqrt(8*S1[i]/g))/2)
	else:
		x1 = -sqrt((1/(2*g)+sqrt(8*S1[i]/g))/2)
		x2 = sqrt((1/(2*g)+sqrt(8*S1[i]/g))/2)	
	integrand1 = lambda x: 1/(sqrt(2*S1[i]-g*(x**2-1/(4*g))**2))
	integrand2 = lambda x: x**2/(sqrt(2*S1[i]-g*(x**2-1/(4*g))**2))
	T = 2*re(quad(integrand1,[x1,x2]))
	mom.append(2*re(quad(integrand2,[x1,x2])/T))
plt.plot(S1[1:],mom,color= "black",linestyle="--")
	
mom_min, mom_max = 0, mom[-1]	
length2 = fsub(mom_max,mom_min)
S2 = MakingGrid(mom_min,mom_max,[],fdiv(length2,N))
results = [[],[],[],[],[],[],[],[],[]]

while K<=K_max:
	S1 = []
	for i in range(0,int(len(bounds[K-K_min])/2)):
		if bounds[K-K_min][2*i]!=bounds[K-K_min][2*i+1]:
			S1.append(MakingGrid(bounds[K-K_min][2*i],bounds[K-K_min][2*i+1],[],
			fdiv(fsub(bounds[K-K_min][2*i+1],bounds[K-K_min][2*i]),N)))
		else:
			S1.append([bounds[K-K_min][2*i]])

	#reducing the search space
	for k in range(0,len(S1)):
		for i in range(0,len(S1[k])):
			for j in range(0,len(S2)): 
				O = RecursionRelation([],S1[k][i],S2[j],g,K)
				try:
					M = HankelMatrix(O)
					L = cholesky(M)
				except ValueError:
					S2[j] = 0
			(results[K-K_min]).append([S1[k][i],S2])
			S2 = MakingGrid(mom_min,mom_max,[],fdiv(length2,N))	
	print(results)
	print(f"K = {K} \n")
	K += 1
	
#plotting
red = 0.2
for i in range(0,len(results)):
	for j in range(0,len(S2)):
		x = []
		for k in range(0,len(results[i])):
			if results[i][k][1][j]!=0:
				x.append(results[i][k][0])
		if x!=[]:
			plt.plot(x,[S2[j]]*len(x),
			color = (red,0,0),marker = "o",linestyle = "None")
	red += 0.1

plt.legend(["Classical","Allowed regions"], loc="lower right",fontsize=14)
plt.title("Bootstrap at $g = 0.05$",fontsize=16)
plt.xlabel("$E$",fontsize=14)
plt.ylabel(r"$\langle x^2 \rangle$",fontsize=14)
plt.grid(linestyle = '-', linewidth = 0.5)
plt.savefig("S_reduced_DoubleWell2.pdf", dpi=300)