Kronig model for the article.

kronig.py

import matplotlib.pyplot as plt
from scipy.linalg import eigh
import numpy as np
#Number of barriers
nbar = 10
#Distance between barriers
a = 1
#Width of the barriers
b = 1/6
#The height of the constant potential
V0 = 100
#Array of barrier locations
xi = np.zeros(nbar)
#Length of our box
L=10
#Number of basis functions. See what happens with fewer or more of them!
npoints = 100
#Our Hamiltonian
H = np.zeros((npoints,npoints))
#This for-loop calculates the f_nm function as explained in the text
def fnm(x,n,m):
	if n==m:
		return x/L - np.sin(2*np.pi*n*x/L)/(2*np.pi*n)
	else:
		return np.sin((m-n)*np.pi*x/L)/(np.pi*(m-n)) - np.sin((m+n)*np.pi*x/L)/(np.pi*(m+n))
def hnm(x,n,m):
	return fnm(x+b/2,n,m) - fnm(x-b/2,n,m)        
#We set up the locations of the barriers.
for i in range(nbar):
	xi[i] = (-0.5)*a + (i+1) * a
#Creating the hamiltonian. Notice that I've used units where
#hbar^2/2m = 1
for i in range(npoints):
	for j in range(i,npoints):
		for k in range(nbar):
			H[i,j] += V0*hnm(xi[k],i+1,j+1)

		if i == j:
			H[i,j] += np.pi**2*(i+1)**2/(L**2)
		H[j,i] = H[i,j]

e, vec = eigh(H)
x = np.zeros(30)
for i in range(30):
	x[i] = (i+1)*np.pi/L
plt.xlabel("Wave number k of solution")
plt.ylabel("Energy")
plt.scatter(x/np.pi,e[0:30]/np.pi**2)
plt.show()