1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
|
# Magnetic field of a toroid, using Biot-Savart law
# Ruth Chabay, Carnegie Mellon University, 2000-05-30
from visual import *
scene.width=800
scene.height=800
scene.title="Magnetic Field of a Toroid"
scene.minscale=.95
bigradius = 0.25
smallradius=0.03
I=1.
muthing=1e-7
nturns=75. # number of loops in helix
nsteps=20. # number of straight segments to approximate 1 turn
Bideal=muthing*2.*nturns*I/bigradius # Ampere's law, inside helix
bscale = bigradius/(Bideal*4.)
# create toroid and arrows at observation locations
toroid = curve(color=(1.,.7,.2), radius=0.003)
tr = vector(bigradius, 0, 0)
dtheta=2*pi/(nsteps*nturns)
phi= 2*pi/nsteps
# use invisible arrows here because vector.rotate illegal at present
a=arrow(pos=(0,0,0), axis=tr, visible=0)
b=arrow(pos=a.axis, axis=(smallradius,0,0), visible=0)
zaxis=vector(0,0,1)
ocount=(nsteps*nturns/10.) # in order to add zeroth observation loc
Barrows=[]
for t in arange(0, nsteps*nturns+1):
a.rotate(angle=dtheta, axis=(0,0,1))
b.rotate(angle=phi, axis=cross(zaxis,a.axis))
toroid.append(pos=a.axis+b.axis)
ocount=ocount+1
if ocount > (nsteps*nturns/10.):
ocount=0
# outside toroid (note that B is *very* small here)
c=arrow(pos=1.2*a.axis, axis=(0,0,0), color=(0,.6, 1), shaftwidth=0.005)
# inside the wire helix (B proportional to 1/r)
Barrows.append(c)
c=arrow(pos=1.07*a.axis, axis=(0,0,0), color=(0,.6, 1), shaftwidth=0.005)
Barrows.append(c)
c=arrow(pos=a.axis, axis=(0,0,0), color=(0,.6, 1), shaftwidth=0.005)
Barrows.append(c)
c=arrow(pos=0.93*a.axis, axis=(0,0,0), color=(0,.6, 1), shaftwidth=0.005)
Barrows.append(c)
# inside toroid (small, in +z direction, due to loop current)
c=arrow(pos=0.8*a.axis, axis=(0,0,0), color=(0,.6, 1), shaftwidth=0.005)
Barrows.append(c)
# uncomment to visualize r
#rarrow=arrow(pos=(0,0,0), axis=(0,0,0), color=(1,0,1), shaftwidth=0.01)
# visualize dl as integration progresses
dlarrow=arrow(pos=(0,0,0), axis=(0,0,0), color=(0,1,0), shaftwidth=0.005)
oldpt = toroid.pos[0]
for pt in toroid.pos:
dl=vector(pt-oldpt)
source=oldpt+dl/2.
dlarrow.axis=dl*2
dlarrow.pos=source
for b in Barrows:
r=vector(b.pos-source)
# uncomment to visualize r
# rarrow.pos=source
# rarrow.axis=r
# rate(30)
dB=muthing*I*cross(dl,r)/r.mag**3
b.axis=b.axis+dB*bscale
oldpt=pt
dlarrow.visible=0
print "computed value at R: ",Barrows[2].axis.mag/bscale
print "Ampere's law value: ", Bideal
|
