from visual import *

print """
Right button drag to rotate view.
Left button drag up or down to move in or out.
"""
scene.title = "Magnetic field due to a coil"
scene.center = vector(0,0,0)
scene.lights = [vector(0.5,0.5,0.5)]
scene.ambient = 0.3
#scene.forward = vector(0,1,0)

n = 10.
length = 1.
radius = 0.5
a = 2.*pi*n/length
c = curve(z=arange(-length/2.,length/2.,length/(20.*n)),color=(1,0,0))
c.x = radius*cos(a*c.z)
c.y = radius*sin(a*c.z)

dtheta = 2.*pi/20.
dlz = dtheta*length/(n*2*pi)
max = 5.
scale = 0.1

for i in arange(-max,max+1,1.):
   for j in arange(-max,max+1,1.):
      for k in arange(-2.*max,2.*max+1,1.):
         bx = 0.
         by = 0.
         bz = 0.
         for theta in arange(-pi*n,pi*n,dtheta):
            dlx = -radius*dtheta*sin(theta)
            dly = radius*dtheta*cos(theta)      
            rx = i/max-radius*cos(theta)
            ry = j/max-radius*sin(theta)
            rz = k/max-length*theta/(2.*n*pi)
            r = sqrt(rx**2+ry**2+rz**2)
            bx = bx+(dly*rz-dlz*ry)/r**3
            by = by+(dlz*rx-dlx*rz)/r**3
            bz = bz+(dlx*ry-dly*rx)/r**3
         arrow(pos=(i/max,j/max,k/max),axis=(scale*bx,scale*by,scale*bz),shaftwidth=0.02)