{VERSION 4 0 "IBM INTEL NT" "4.0" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 
1 0 0 0 0 1 }{CSTYLE "" -1 256 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 257 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 
258 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 259 "" 0 1 0 0 
0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 1 }{CSTYLE "" -1 261 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 262 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 
263 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 264 "" 0 1 0 0 
0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 1 }{CSTYLE "" -1 266 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 267 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 
268 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 269 "" 0 1 0 0 
0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 270 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 1 }{CSTYLE "" -1 271 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 272 "" 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 
273 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 274 "" 0 1 0 0 
0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 275 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 1 }{CSTYLE "" -1 276 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 277 "" 1 18 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 
278 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 279 "" 0 1 0 0 
0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 280 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 1 }{CSTYLE "" -1 281 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 282 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 
283 "" 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }{PSTYLE "Normal" -1 0 1 
{CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }0 0 0 -1 -1 -1 0 
0 0 0 0 0 -1 0 }{PSTYLE "Heading 1" 0 3 1 {CSTYLE "" -1 -1 "" 1 18 0 
0 0 0 0 1 0 0 0 0 0 0 0 1 }1 0 0 0 8 4 0 0 0 0 0 0 -1 0 }}
{SECT 0 {EXCHG {PARA 0 "" 0 "" {TEXT 277 15 "Magnetic Bottle" }}{PARA 
0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 586 "In this worksh
eet we investigate charged-particle trajectories in magnetic fields. T
he magnetic field of a current loop was calculated in the worksheet Cu
rrentloop.mws. The magnetic field lines were calculated using a solver
 for systems of first-order differential equations. Then the current l
oop results were extended by combining the fields to form a Helmholtz \+
pair, and a solenoid. Here we investigate charged-particle trajectorie
s by integrating Newton's equations with the Lorentz force. We begin w
ith the study of trajectories in a magnetic field caused by a single c
urrent loop." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 
-1 240 "Subsequently we assemble a magnetic field made up of four curr
ent loops. Two of the loops are larger and form a Helmholtz coil. Two \+
loops that are smaller are added to generate a bottleneck. This will b
e a simple model of a magnetic bottle." }}{PARA 0 "" 0 "" {TEXT -1 0 "
" }}{PARA 0 "" 0 "" {TEXT -1 89 "In principle, we work in MKSA (SI) un
its. We do, however use a charge-over-mass ratio of " }{TEXT 282 1 "e
" }{TEXT -1 1 "/" }{TEXT 281 1 "m" }{TEXT -1 231 "=1, i.e., we do not \+
specify the mass as 1 kg, and the charge unit as 1.6 * 10^(-19) C. For
 real-life applications this would have to be corrected. For the quali
tative analysis of charged-particle trajectories this is not a problem
." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 0 {PARA 3 "
" 0 "" {TEXT 272 60 "Charged-particle trajectories in the field of a c
urrent loop" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 391 "We carry out the s
olution of Newton's equation for a charged particle with the Lorentz f
orce. The usual textbook discussions deal only with homogeneous fields
, where the constant magnetic field imposes circular motion in a plane
 perpendicular to the magnetic field direction. In general this leads \+
to spiralling motion for particles which possess a velocity component \+
aligned with the field." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 
0 "" {TEXT -1 697 "The motion of charged particles in inhomogeneous ma
gnetic fields is much more complicated. We are particularly interested
 in a demonstration that charged particles which move in the direction
 where the magnetic field becomes stronger are slowed down and turn ar
ound, i.e., that they prefer to travel in the direction where the magn
etic field becomes weaker. This has implications for cosmic charged pa
rticles (picking up energy in galactic magnetic fields), and charged p
articles travelling from the Sun to the Earth, which spiral towards th
e poles (and turn around where the field becomes too strong), and coll
ide with each other to produce optical airglow emissions (aurora, or n
orthern lights)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 19 "When we calculated " }{TEXT 261 1 "B" }{TEXT -1 1 "_" }
{TEXT 260 1 "x" }{TEXT -1 101 " we really obtained the transverse comp
onent of the magnetic field, as we worked specifically in the " }
{TEXT 263 1 "x" }{TEXT -1 1 "-" }{TEXT 262 1 "z" }{TEXT -1 135 " plane
. For a solution to Newton's equations we should generate three Cartes
ian components to be calculated anywhere in space (for any " }{TEXT 
266 1 "x" }{TEXT -1 1 "," }{TEXT 265 1 "y" }{TEXT -1 1 "," }{TEXT 264 
1 "z" }{TEXT -1 192 "). Given how the Lorentz force acts (at right ang
les to the velocity and to the B field) we cannot expect the motion to
 be confined to a plane for the case of a non-homogeneous magnetic fie
ld." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 208 "W
e define a procedure to calculate the Cartesian field components for t
he particular case of a current loop defined in Bxa,Bza. We start over
 and define the functions needed for the magnetic field components." }
}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "restart; with(linalg): wit
h(plots):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "Qp:=(Q+R)^2+a^
2: Qm:=(Q-R)^2+a^2:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "Qr:=2*sqrt(Q
*R/Qp);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 83 "Bxex:=i_L*a/(5*1
0^(6))/(sqrt(Qp)*Q)*((R^2+Q^2+a^2)*EllipticE(Qr)/Qm-EllipticK(Qr));" }
}{PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "Bxa:=unapply(subs(i_L=1000,R=0.005
,Bxex),Q,a):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 77 "Bzex:=i_L/(
5*10^(6))/sqrt(Qp)*(EllipticE(Qr)*(R^2-Q^2-a^2)/Qm+EllipticK(Qr));" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "Bza:=unapply(subs(i_L=1000,R=0.005,
Bzex),Q,a):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 120 "To avoid the prob
lems associated with taking the limit at x=0,z=0, we add a tiny consta
nt to the cylindrical coordinate " }{TEXT 276 1 "r" }{TEXT -1 1 "." }}
}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "eta:=10^(-Digits):" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 71 "Bfield:=proc(x,y,z) local r,
Btr,Bx,By,Bz,cphi,sphi; global Bxa,Bza,eta;" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 21 "r:=sqrt(x^2+y^2)+eta;" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 42 "Bz:=evalf(Bza(z,r)); Btr:=evalf(Bxa(z,r));" }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 35 "cphi:=evalf(x/r); sphi:=evalf(y/r);" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 27 "By:=Btr*sphi; Bx:=Btr*cphi;" }}{PARA 0 "> " 0 
"" {MPLTEXT 1 0 16 "[Bx,By,Bz]; end:" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 26 "Bfield(0.002,0.003,0.004);" }}}{EXCHG {PARA 0 "> " 0 
"" {MPLTEXT 1 0 23 "Bfield(0.0,0.00,0.004);" }}}{EXCHG {PARA 0 "> " 0 
"" {MPLTEXT 1 0 307 "NE:=diff(x(t),t)=p_x(t),diff(p_x(t),t)=p_y(t)*Bfi
eld(x(t),y(t),z(t))[3]-p_z(t)*Bfield(x(t),y(t),z(t))[2],diff(y(t),t)=p
_y(t),diff(p_y(t),t)=p_z(t)*Bfield(x(t),y(t),z(t))[1]-p_x(t)*Bfield(x(
t),y(t),z(t))[3],diff(z(t),t)=p_z(t),diff(p_z(t),t)=p_x(t)*Bfield(x(t)
,y(t),z(t))[2]-p_y(t)*Bfield(x(t),y(t),z(t))[1]:" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT -1 29 "We start a trajectory in the " }{TEXT 268 1 "x" }
{TEXT -1 1 "-" }{TEXT 267 1 "z" }{TEXT -1 31 " plane with no momentum \+
in the " }{TEXT 271 1 "x" }{TEXT -1 48 "-direction, and non-vanishing \+
components in the " }{TEXT 270 1 "y" }{TEXT -1 6 "- and " }{TEXT 269 
1 "z" }{TEXT -1 12 "-directions." }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 63 "IC:=x(0)=0.02,p_x(0)=0,y(0)=0,p_y(0)=0.001,z(0)=0,p_z
(0)=0.001;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 88 "sol:=dsolve(
\{NE,IC\},\{x(t),p_x(t),y(t),p_y(t),z(t),p_z(t)\},numeric,output=listp
rocedure):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 33 "This problem is sol
ved in Maple6." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "sol(1);" }
}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 131 "For faster calculations and to \+
be able to work in Maple5 we introduce our own Runge-Kutta procedure u
sed before for the fieldlines." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 64 "RKn:=proc(f,ic::list,N,x1,h) local xn,yn,ynew,m1,m2,m3,m4,ii,n
n;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "xn:=op(1,ic):" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 7 "yn:=[]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "f
or nn from 1 to N do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "yn:=[op(yn
),op(1+nn,ic)]: od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "for ii from \+
1 to round((x1-xn)/h) do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "m1:=[]:
" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "for nn from 1 to N do:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "m1:=[op(m1),evalf(f(xn,yn)[nn])]; o
d:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "m2:=[]:" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 22 "for nn from 1 to N do:" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 55 "m2:=[op(m2),evalf(f(xn+h/2,evalm(yn+h/2*m1))[nn])]; od:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "m3:=[]:" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 22 "for nn from 1 to N do:" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 55 "m3:=[op(m3),evalf(f(xn+h/2,evalm(yn+h/2*m2))[nn])]; od:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "xn:=xn+h;" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 7 "m4:=[]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "for nn f
rom 1 to N do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "m4:=[op(m4),evalf
(f(xn,evalm(yn+h*m3))[nn])]; od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "
ynew:=[]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "for nn from 1 to N do:
" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 83 "ynew:=[op(ynew),evalf(yn[nn]+h/
6*(m1[nn]+2*(m2[nn]+m3[nn])+m4[nn]))]; od: yn:=ynew;" }}{PARA 0 "> " 
0 "" {MPLTEXT 1 0 3 "od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "[xn,op(
ynew)];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end:" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT -1 42 "We define a RHS procedure for the problem." }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "fn:=proc(t,yn) local x,y,z,p
_x,p_y,p_z,BF; global Bfield;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 84 "x:
=yn[1]: p_x:=yn[2]: y:=yn[3]: p_y:=yn[4]: z:=yn[5]: p_z:=yn[6]: BF:=Bf
ield(x,y,z);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 74 "[p_x,p_y*BF[3]-p_z*
BF[2],p_y,p_z*BF[1]-p_x*BF[3],p_z,p_x*BF[2]-p_y*BF[1]];" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 4 "end:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 279 "T
he argument list contains the external RHS, a list that starts with th
e initial time (independent variable), and the initial values, the num
ber of dependent variables, the final value of the time (independent v
ariable), and the time-step (step-size of the independent variable)." 
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "
" {TEXT -1 59 "It seems that we can use a large step. (at least intial
ly)." }}{PARA 0 "" 0 "" {TEXT -1 94 "We prepare for a graph of the cur
rent loop which causes the magnetic field and the trajectory:" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 68 "PL4:=spacecurve([0.005*cos(f
),0.005*sin(f),0],f=0..2*Pi,color=blue):" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 113 "We pick initial conditions for a charged particle locate
d above the current loop without initial momentum in the " }{TEXT 273 
1 "z" }{TEXT -1 48 "-direction, and a small amount of motion in the " 
}{TEXT 275 1 "x" }{TEXT -1 1 "-" }{TEXT 274 1 "y" }{TEXT -1 115 " plan
e. The reader is encouraged to play with the initial conditions patien
tly in order to understand the outcomes." }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 108 "Rt:=[]: dt:=0.5: x0:=0.002: px0:=-0.0: y0:=eta: py0:
=-0.00005: z0:=0.005: pz0:=0.0: for i from 1 to 1000 do:" }}{PARA 0 ">
 " 0 "" {MPLTEXT 1 0 168 "t0:=(i-1)*dt; sol:=RKn(fn,[t0,x0,px0,y0,py0,
z0,pz0],6,t0+dt,dt); x0:=sol[2]: px0:=sol[3]: y0:=sol[4]: py0:=sol[5]:
 z0:=sol[6]: pz0:=sol[7]: Rt:=[op(Rt),[x0,y0,z0]]: od:" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "PL5:=pointplot3d(Rt,axes=boxed,colo
r=green,labels=[x,y,z]): " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
17 "display(PL4,PL5);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "Rt
(nops(Rt));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 388 "An interesting fe
ature of Loretnz-force trajectories is that they are not symmetric und
er time reversal. Under time reversal velocities are reversed, and one
 ought to be able to run through the trajectory from the end to the be
ginning. This is usually a check for the accuracy of the integration s
cheme. To accomplish this with the Lorentz force, we have to change th
e sign of the charge." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "q_p
:=-1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "fn:=proc(t,yn) loc
al x,y,z,p_x,p_y,p_z,BF; global Bfield,q_p;" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 84 "x:=yn[1]: p_x:=yn[2]: y:=yn[3]: p_y:=yn[4]: z:=yn[5]:
 p_z:=yn[6]: BF:=Bfield(x,y,z);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 92 "
[p_x,q_p*(p_y*BF[3]-p_z*BF[2]),p_y,q_p*(p_z*BF[1]-p_x*BF[3]),p_z,q_p*(
p_x*BF[2]-p_y*BF[1])];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end:" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 74 "Rt:=[]: dt:=0.5: px0:=-px0: \+
py0:=-py0: pz0:=-pz0: for i from 1 to 1000 do:" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 168 "t0:=(i-1)*dt; sol:=RKn(fn,[t0,x0,px0,y0,py0,z0,pz0],
6,t0+dt,dt); x0:=sol[2]: px0:=sol[3]: y0:=sol[4]: py0:=sol[5]: z0:=sol
[6]: pz0:=sol[7]: Rt:=[op(Rt),[x0,y0,z0]]: od:" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 60 "PL5:=pointplot3d(Rt,axes=boxed,color=green,lab
els=[x,y,z]): " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "display(P
L4,PL5);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "Rt[nops(Rt)];" 
}}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 97 "Our algorithm for the reverse p
ropagation indeed took us back to the beginning of the trajectory." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 119 "Let us p
ropagate further to see where the particle goes. We accumulate the tra
jectory points into the existing list Rt." }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 24 "for i from 1 to 1000 do:" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 168 "t0:=(i-1)*dt; sol:=RKn(fn,[t0,x0,px0,y0,py0,z0,pz0],
6,t0+dt,dt); x0:=sol[2]: px0:=sol[3]: y0:=sol[4]: py0:=sol[5]: z0:=sol
[6]: pz0:=sol[7]: Rt:=[op(Rt),[x0,y0,z0]]: od:" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 59 "PL5:=pointplot3d(Rt,axes=boxed,color=green,lab
els=[x,y,z]):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "display(PL
4,PL5);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "Rt[nops(Rt)];" }
}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 71 "The charged particle spiralled o
ut again, but in a different direction." }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 271 "We observe qualitatively some spiralling of the charged-
particle orbits, and preferred motion in the z direction towards weake
r fields. This can be used for trapping of charged particles (magnetic
 bottle, confinement in fusion reactors), which is discussed further b
elow." }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 35 "A number of open questio
ns remains:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 
-1 232 "1) we have simply followed the mathematical convention for the
 traversal of the current loop. What happens to the magnetic field com
ponents if the current direction (direction in which the curve integra
l is carried out) is reversed?" }}{PARA 0 "" 0 "" {TEXT -1 109 "2) we \+
have considered positively charged particles; what happens if the sign
 of the Lorentz force is changed?" }}{PARA 0 "" 0 "" {TEXT -1 268 "3) \+
most importantly: we have been sloppy as far as the units are concerne
d, i.e., we are just looking qualitatively at the trajectories. For re
al-life applications we would need to understand better what our units
 (m=1, e=1) imply, i.e., how do we convert to SI units?" }}{PARA 0 "" 
0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 0 {PARA 3 "" 0 "" 
{TEXT 283 22 "Simple magnetic bottle" }}{EXCHG {PARA 0 "" 0 "" {TEXT 
-1 173 "We now proceed with the model of a magnetic bottle based on fo
ur current loops. First we restart and define the analytic expressions
 for the magnetic field of a single loop." }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 35 "restart; with(linalg): with(plots):" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 38 "We pick for the Helmholtz coil radius:" }
}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "R_Hh:=0.005;" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "Qp:=(Q+R)^2+a^2: Qm:=(Q-R)^2+a^2:" 
}}{PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "Qr:=2*sqrt(Q*R/Qp);" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 83 "Bxex:=i_L*a/(5*10^(6))/(sqrt(Qp)*Q)
*((R^2+Q^2+a^2)*EllipticE(Qr)/Qm-EllipticK(Qr));" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 45 "Bxa:=unapply(subs(i_L=1000,R=R_Hh,Bxex),Q,a):" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 77 "Bzex:=i_L/(5*10^(6))/sqrt(Qp
)*(EllipticE(Qr)*(R^2-Q^2-a^2)/Qm+EllipticK(Qr));" }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 45 "Bza:=unapply(subs(i_L=1000,R=R_Hh,Bzex),Q,a):" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 250 "We define a second function for t
he loops with smaller radius, which are introduced to confine the char
ged particles. The idea is to create a higher density of magnetic fiel
d lines inside the Helmholtz coils. For the confining field radius we \+
choose:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "R_cf:=0.002;" }}
}{EXCHG {PARA 0 "" 0 "" {TEXT -1 119 "We choose twice the current in t
he confining coils to obtain a large increase in the magnetic field at
 the bottlenecks." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "Bxb:=u
napply(subs(i_L=2000,R=R_cf,Bxex),Q,a):" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 45 "Bzb:=unapply(subs(i_L=2000,R=R_cf,Bzex),Q,a):" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 58 "To avoid the problems associated with tak
ing the limit at " }{TEXT 279 1 "x" }{TEXT -1 3 "=0," }{TEXT 280 2 " z
" }{TEXT -1 57 "=0, we add a tiny constant to the cylindrical coordina
te " }{TEXT 259 1 "r" }{TEXT -1 1 "." }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 18 "eta:=10^(-Digits):" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 71 "Bfield:=proc(x,y,z) local r,Btr,Bx,By,Bz,cphi,sphi; g
lobal Bxa,Bza,eta;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "r:=sqrt(x^2+y
^2)+eta;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 169 "Bz:= evalf(Bza(z+0.5*R
_Hh,r)+Bza(z-0.5*R_Hh,r)+Bzb(z+0.5*R_Hh,r)+Bzb(z-0.5*R_Hh,r)); Btr:=ev
alf(Bxa(z+0.5*R_Hh,r)+Bxa(z-0.5*R_Hh,r)+Bxb(z+0.5*R_Hh,r)+Bxb(z-0.5*R_
Hh,r));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "cphi:=evalf(x/r); sphi:=
evalf(y/r);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "By:=Btr*sphi; Bx:=Bt
r*cphi;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "[Bx,By,Bz]; end:" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "Bfield(0.002,0.003,0.004);" 
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "Bfield(0.0,0.00,0.004);" 
}}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 168 "For fast calculations and to b
e able to work in Maple5 we introduce our own Runge-Kutta procedure us
ed before for the fieldlines rather than relying on dsolve[numeric]." 
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "RKn:=proc(f,ic::list,N,x1
,h) local xn,yn,ynew,m1,m2,m3,m4,ii,nn;" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 13 "xn:=op(1,ic):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "yn:=[]:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "for nn from 1 to N do:" }}{PARA 0 "
> " 0 "" {MPLTEXT 1 0 29 "yn:=[op(yn),op(1+nn,ic)]: od:" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 37 "for ii from 1 to round((x1-xn)/h) do:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "m1:=[]:" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 22 "for nn from 1 to N do:" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 37 "m1:=[op(m1),evalf(f(xn,yn)[nn])]; od:" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 7 "m2:=[]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "for nn f
rom 1 to N do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "m2:=[op(m2),evalf
(f(xn+h/2,evalm(yn+h/2*m1))[nn])]; od:" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 7 "m3:=[]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "for nn from 1 to \+
N do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "m3:=[op(m3),evalf(f(xn+h/2
,evalm(yn+h/2*m2))[nn])]; od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "xn:
=xn+h;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "m4:=[]:" }}{PARA 0 "> " 0 
"" {MPLTEXT 1 0 22 "for nn from 1 to N do:" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 49 "m4:=[op(m4),evalf(f(xn,evalm(yn+h*m3))[nn])]; od:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "ynew:=[]:" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 22 "for nn from 1 to N do:" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 83 "ynew:=[op(ynew),evalf(yn[nn]+h/6*(m1[nn]+2*(m2[nn]+m3[nn])+m4[
nn]))]; od: yn:=ynew;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "[xn,op(ynew)];" }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 4 "end:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 42 "We defi
ne a RHS procedure for the problem." }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 57 "fn:=proc(t,yn) local x,y,z,p_x,p_y,p_z,BF; global Bfi
eld;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 84 "x:=yn[1]: p_x:=yn[2]: y:=yn
[3]: p_y:=yn[4]: z:=yn[5]: p_z:=yn[6]: BF:=Bfield(x,y,z);" }}{PARA 0 "
> " 0 "" {MPLTEXT 1 0 74 "[p_x,p_y*BF[3]-p_z*BF[2],p_y,p_z*BF[1]-p_x*B
F[3],p_z,p_x*BF[2]-p_y*BF[1]];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "en
d:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 279 "The argument list contains
 the external RHS, a list that starts with the initial time (independe
nt variable), and the initial values, the number of dependent variable
s, the final value of the time (independent variable), and the time-st
ep (step-size of the independent variable)." }}}{EXCHG {PARA 0 "" 0 "
" {TEXT -1 220 "It seems that we can use a large step dt=0.5. The read
er should check occasionally whether a reduction in the time step resu
lts in reproduced trajectories. This is the price to be paid for using
 a non-adaptive algorithm." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 
"" 0 "" {TEXT -1 94 "We prepare for a graph of the current loop which \+
causes the magnetic field and the trajectory:" }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 73 "PL1:=spacecurve([R_Hh*cos(f),R_Hh*sin(f),0.5*R_H
h],f=0..2*Pi,color=blue):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 74 "PL2:=s
pacecurve([R_Hh*cos(f),R_Hh*sin(f),-0.5*R_Hh],f=0..2*Pi,color=blue):" 
}}{PARA 0 "> " 0 "" {MPLTEXT 1 0 72 "PL3:=spacecurve([R_cf*cos(f),R_cf
*sin(f),0.5*R_Hh],f=0..2*Pi,color=red):" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 73 "PL4:=spacecurve([R_cf*cos(f),R_cf*sin(f),-0.5*R_Hh],f=0..2*Pi,
color=red):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 121 "We pick initial c
onditions for a charged particle located in between the current loops \+
with some initial momentum in the " }{TEXT 256 1 "z" }{TEXT -1 48 "-di
rection, and a small amount of motion in the " }{TEXT 258 1 "x" }
{TEXT -1 1 "-" }{TEXT 257 1 "y" }{TEXT -1 154 " plane. We start the tr
ajectory in the middle of the bottle (z=0) so that we are in a relativ
ely weaker magnetic field region with our initial conditions." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 107 "The read
er is encouraged to play with the initial conditions patiently in orde
r to understand the outcomes." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 35 "p_long:=-0.00004; p_trans:=0.00002;" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 109 "Rt:=[]: dt:=0.50: x0:=0.0005: px0:=-0.0: y0:=eta: \+
py0:=p_trans: z0:=0.0: pz0:=p_long: for i from 1 to 400 do:" }}{PARA 
0 "> " 0 "" {MPLTEXT 1 0 168 "t0:=(i-1)*dt; sol:=RKn(fn,[t0,x0,px0,y0,
py0,z0,pz0],6,t0+dt,dt); x0:=sol[2]: px0:=sol[3]: y0:=sol[4]: py0:=sol
[5]: z0:=sol[6]: pz0:=sol[7]: Rt:=[op(Rt),[x0,y0,z0]]: od:" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "PL5:=pointplot3d(Rt,axes=boxed,colo
r=green,labels=[x,y,z]): " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
29 "display(PL1,PL2,PL3,PL4,PL5);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 
196 "p_long=-0.00004; p_trans=0.00002 work. Too much longitudinal mome
ntum: violates the condition below. Too much transverse momentum gets \+
the charged particle out of the strong-field region sideways." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 16 "J.D. Jack
son in " }{TEXT 278 25 "Classical Electrodynamics" }{TEXT -1 270 " der
ives an expression for the trapping condition. The ratio of the longit
udinal to transverse initial velocity of the particle has to be bounde
d from above by an expression that involves the ratio of magnetic fiel
d strengths at maximum as compared to the middle region." }}{PARA 0 "
" 0 "" {TEXT -1 31 "v_L(0)/v_T(0) < sqrt(B_max/B-1)" }}{PARA 0 "" 0 "
" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 250 "The above example show
s that the trajectories can be complicated. We can check whether the c
onfinement result is consistent with the above formula. The reader sho
uld change the initial conditions to see that not all of them lead to \+
particle trapping." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "abs(p
_long/p_trans);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "B0:=norm
(Bfield(0.,0.,0.));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "B_m:
=norm(Bfield(eta,eta,0.0024));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 15 "sqrt(B_m/B0-1);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 70 "We can
 make the particle escape by exceeding this ratio substantially:" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "p_long:=-0.00008; p_trans:=0
.00002;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 109 "Rt:=[]: dt:=0.5
0: x0:=0.0000: px0:=-0.0: y0:=eta: py0:=p_trans: z0:=0.0: pz0:=p_long:
 for i from 1 to 400 do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 168 "t0:=(i
-1)*dt; sol:=RKn(fn,[t0,x0,px0,y0,py0,z0,pz0],6,t0+dt,dt); x0:=sol[2]:
 px0:=sol[3]: y0:=sol[4]: py0:=sol[5]: z0:=sol[6]: pz0:=sol[7]: Rt:=[o
p(Rt),[x0,y0,z0]]: od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "P
L5:=pointplot3d(Rt,axes=boxed,color=green,labels=[x,y,z]): " }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "display(PL1,PL2,PL3,PL4,PL5)
;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 ">
 " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "
" }}}}{MARK "2" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 
2 33 1 1 }