From Werner KRAUTH
import math, cmath, random
random.seed('Werner')
N = 8 # take even
beta = math.sqrt(2.0)
PhiR = []
for k in range(N):
PhiR.append(random.uniform(-1.0, 1.0)) # this is an example
rvalues = range(N)
kvalues = [2.0 * math.pi / N * float(l) for l in range(N)]
epskvalues = [4.0 * math.sin( k / 2.0) **2 for k in kvalues]
Energy_real = -2.0 * N + (PhiR[N - 1] - PhiR[0]) ** 2
for r in range(N - 1):
Energy_real += (PhiR[r + 1] - PhiR[r]) ** 2
print Energy_real, 'Energy of the conf. in real representation'
#
# Fourier transform (by hand, no FFT)
#
PhiK = []
for fourierk in kvalues:
dummy = 0.0
for i in range(N):
term = cmath.exp(- 1j * fourierk * rvalues[i]) * PhiR[i] / math.sqrt(N)
dummy += term
PhiK.append(dummy)
#
# real-valued Fourier coefficients, don't forget to count PsiK[N/2] = PhiK[N/2]
#
PsiKpos = [0.0]
PsiKneg = [0.0]
for k in range(1,N / 2):
PsiKpos.append((PhiK[k] + PhiK[N - k]) / math.sqrt(2.0))
PsiKneg.append((PhiK[k] - PhiK[N - k]) / (math.sqrt(2.0) * 1j))
#
# Express the energy in Fourier modes
#
Energy_fourier = - 2.0 * N + 4.0 * PhiK[N / 2] ** 2
for k in range(1, N / 2):
Energy_fourier += 2.0 * epskvalues[k] * PhiK[k] * PhiK[N - k]
print Energy_fourier, 'energy of the conf. in Fourier representation'
#
# Energy in Fourier representation with real-valued terms
#
Energy_real = - 2.0 * N + 4.0 * PhiK[N / 2]**2
for k in range(1, N / 2):
Energy_real += (PsiKneg[k] ** 2 + PsiKpos[k] ** 2) * epskvalues[k]
print Energy_real, ' energy of th conf. with real-valued Fourier coefficients'
#
# Correlation between site 0 and site k
#
print 'check correlation in real and (real Fourier)'
for r in range(N):
dummy = PhiK[N / 2] * (1.0 - (-1) ** r) / math.sqrt(N)
for k in range(1, N / 2):
term = PsiKpos[k] * (1.0 - math.cos(kvalues[k] * r)) + \
PsiKneg[k] * math.sin(kvalues[k] * r)
dummy += term * math.sqrt(2.0 / N)
print r, dummy.real, PhiR[0] - PhiR[r], ' r, Fourier term, direct term'
xdata = []
ydata = []
for r in range(0, N / 2 + 1):
summe = (1.0 - (-1) ** r) ** 2 / (32.0)
for k in range(1, N / 2):
summe += math.sin(math.pi * k * r / N) ** 2 / (4.0 * math.sin(math.pi *
k / N) ** 2)
ydata.append( math.exp(- 2.0 * summe / N / beta))
xdata.append(r)
print ydata