from pylab import *    # needs Matplotlib

r = 1
res = 500	# grid resolution (100 seems to be sufficient here and will create a smaller SVG file)
figure(figsize=(6, 6))
axes((0,0,1,1), frameon = False)

def dist3(a,b,c,d,e,f):
	return maximum(r, sqrt((a-d)**2 + (b-e)**2 + (c-f)**2))

x = linspace(-150, 150, res)
y = linspace(-150, 150, res)
X,Y = meshgrid(x,y)
F = zeros((res,res,3))

# loop over two coils
for coils in 1,-1:
	# sum field contributions from coil in 10-degree steps
	for p in range(0, 360, 10):
		xc,yc,zc = (100*sin(pi*p/180),50*coils,100*cos(pi*p/180))
		MAG = 1/((r+dist3(X,Y,0,xc,yc,zc))**3)
		# (We leave out the necessary constants that would be required
		# to get proper units because only scaling behavior will be shown
		# in the plot. This is also why a sum instead of an integral
		# can be used.)
		Z = cross((X[:,:,newaxis]-xc,Y[:,:,newaxis]-yc,-zc),(-zc,0,xc))
		F += Z[0,:,:]*MAG[:,:,newaxis]

B = sqrt(F[:,:,0]**2+F[:,:,1]**2+F[:,:,2]**2)
# scale field strength by value at center:
B = B / B[res/2,res/2]

contour(x,y,B, levels = (.5,.8,.9,.95,.99,1.01,1.05,1.1))

# add wire symbols:
scatter((100,100,-100,-100), (50,-50,50,-50), s = 400, color = "black")

axis((-130, 130, -130, 130))
xticks([])
yticks([])

savefig("Helmholtz_coil,_B_magnitude_cross_section.svg")

show()