# copyright 2000 Oregon State University              Capstones in
# Physics: Electromagnetism                  Philip J. Siemens
# 
# MAGNETIC FIELD OF A RING
# 
# This Maplesheet is about the magnetic field of a localized static
# current distribution.
# 
# In the class notes, we calculate the field far away from the sources.
# 
# Here we use numerical methods to view the pattern of the field
# everywhere.
# 
> restart: with (linalg): with (plots):
# 
# BIOT SAVART LAW
# 
# The simplest finite current distribution is a loop of current  I
# circulating around a circle of radius R.
# 
# From the principle of superposition we know that the magnetic field
# for the loop is just the vector sum of the magnetic fields from each
# part of the ring.
# 
# For a current I, the constant of proportionality for the magnetic
# field is  k = mu[0]*I/4/pi.
# 
# Then the magnetic field dB at a position r due to a segment dL  of the
# ring at position R is
# 
#     dB = k (dL x  (r-R))/(r-R)^3   
> 
# 
# SYMMETRY CONSIDERATIONS
# 
# By symmetry, 
#  
#          the total magnetic field B depends only on 
# 
#  the distance r  of the observer from the center of the ring     
# and 
#  the angle theta  of the observer from the ring's axis of symmetry
# 
#          the direction of  B is in the plane which contains the axis
# and r
# 
# We will need to use these symmetry considerations to simplify the
# computation. 
#   
# Otherwise, 
# 
#  the math expressions would get really complicated              
# and
#  it would take too much time for Maple to calculate the fields
# 
# Due to the symmetry of the ring, it is convenient to choose a
# coordinate system in which 
# 
#  the center of the ring is the origin
# and
#  the axis of the ring is the z-axis    
# 
# We will assume that, in the first quadrant,  the current is flowing
# from the x-axis toward the y-axis.
# 
# 
# It will be convenient to switch back and forth between cartesian and
# spherical coordinates.
# 
# 
# We can start by finding the field at a point r  in the  x - z   plane.
#  
# Then we will find the field elsewhere by applying the symmetry
# arguments.
# 
# To set up the integral we choose cartesian coordinates, 
#        because then the unit vectors do not depend on the coordinates.
# 
# The observer's position is
# 
> r:=[x,0,z];
# 
# The position of the current segment on the loop is given by its
# angular coordinate phi:
# 
> R:=[a*cos(phi),a*sin(phi),0];
# 
# 
# With these substitutions, the square of the distance between the
# source and observer is 
# 
> Rr2:=dotprod(r-R,r-R);
# 
# This can be simplified:
# 
> Rr2:=simplify(Rr2);
# 
# 
# The differential segment dL  on the loop is proportional to the
# differential dphi ; 
# the ratio of proportionality is the vector
# 
> dLoverdphi:=a*[-sin(phi),cos(phi),0];
# 
# 
# The numerator in the Biot-Savart law is, except for a factor   dphi  ,
# 
> numerator:=crossprod(dLoverdphi,(r-R));
# 
# This can be simplified: 
# 
> numerator:=[numerator[1],numerator[2],simplify(numerator[3])];
# 
# 
# We are now ready to do the integrals. 
#  
# 
# NUMERICAL INTEGRATIONS
# 
# The integrals can, after some manipulation, be converted to elliptic
# integrals.
# 
# Instead, we will do them numerically.
# 
# 
# We first define the integrands.  
# 
> dBx:=k*numerator[1]/Rr2^(3/2);
# 
# 
> dBz:=k*numerator[3]/Rr2^(3/2);
# 
# There is no point in doing the integral for B[y]  , because it is zero
# by symmetry.
# 
# 
# We need to give values to the ring's radius and current.  Let's take
# a=1,k=1
# 
> dBx:=subs(a=1,k=1,dBx); dBz:=subs(a=1,k=1,dBz);
# 
# 
# 
# Next, we set up the integrals, using the integration routine int. 
#  
# If we use int in its usual format, Maple will try to do the integral
# analytically.
# 
# We can save it this step by capitalizing the first letter.
# 
# We can also save time by halving the range of integration from 0 to pi
#  , then doubling the answer by doubling the integrand:
# 
# 
> Bx:=Int(2*dBx,phi=0..Pi); Bz:=Int(2*dBz,phi=0..Pi);
# 
# These integrals will take some time for Maple to evaluate numerically.
# 
# To save time, we will plot them in the x - z  plane.  We can do this
# because:
# 
#  For y=0, the magnetic field only has x and z components. 
# 
#  We can find the fields elsewhere by symmetry.                  
# 
# 
# 
# To plot the fields, we have to introduce procedures defining the
# integrals.
# 
# We will also save time by limiting the integration procedure to 0.1% =
# 10^(-3)
# 
#   
#  -> 
# 
> Bxp:=(x,z) -> 
> evalf(Int(2*cos(phi)*z/(-2*cos(phi)*x+x^2+1+z^2)^(3/2),phi = 0 ..
> Pi),3);
# 
> Bzp:=(x,z) -> 
> evalf(Int(2*(1-cos(phi)*x)/(-2*cos(phi)*x+x^2+1+z^2)^(3/2),phi = 0 ..
> Pi),3);
# 
# 
# 
# We can save time by plotting the result in just one quadrant.  The
# others follow from symmetry.
# 
> fieldplot([Bxp,Bzp],0.5..2.5,0.5..2.5,grid=[5,5],arrows=slim);
# 
# 
> 
# 
# QUESTIONS FOR HOMEWORK
# 
# Question 1.
#   
# Sketch the extension of the fieldplot to all 4 quadrants of the xz
# plane.
# 
# Question 2. 
#  
# What goes wrong if we change the limits of the plot to include the x
# and z axes, for example to 0..2,0..2  ?  Explain, both physically and
# numerically.
# 
# Question 3. 
# 
# Choose a larger scale, 0.5..5.5,0.5..5.5 .  
# Multiply the fields by an appropriate power of r so that the
# magnitudes of the field vectors aren't falling off with increasing r. 
# What is the power you need?  Explain.  Show an output graph which
# illustrates your answer.  What is happening at small r? Why is this
# okay?
#  
# 
# 
# 
# version 17 October 2001