import numpy as np
import itertools

N = 20

#find the unique distances to lattice points from the origin
r = np.empty(((N+1)**3,3),np.float)

m = 0
for i in range(N+1):
	for j in range(N+1):
		for k in range(N+1):
			r[m] = [i,j,k]
			m += 1

r = np.sqrt(np.unique((r**2).sum(axis=1)))

#then loop through them evaluating the unconverged Madelung constant
for rad in r:
	n = np.int(np.ceil(rad))
	M = 0.
	for i in range(-n,n+1):
		for j in range(-n,n+1):
			for k in range(-n,n+1):
				if  np.sqrt(i**2+j**2+k**2) <= rad and not(i==0 and j ==0 and k == 0):
					M += (-1)**(i+j+k+1)/np.sqrt(i**2+j**2+k**2)
	print str(rad)+'\t'+str(M)