! compute the magnetic field from a circular current loop
! Program to accompany "Computational Physics" by N. Giordano
! Copyright Prentice Hall 1997
program current_loop
   option nolet
   library "sglib*","sgfunc*"
   dim field(0,0,0)      ! magnetic field at grid location i,j,k
   call initialize(field,max)
   call calculate(field,max)
   call display(field,max)
end
! initialize variables
! True basic initializes all variables to zero, so no need to zero out 
! field array        max determines number of grid steps
sub initialize(field(,,),max)
   input prompt "number of grid steps -> ": max
   mat redim field(3,-max to max,-max to max)  ! resize field array
end sub
! sweep through the lattice here  
sub calculate(field(,,),max)
   for i = -max to max        ! compute field in x-z plane
      for k = -max to max
         call calculate_field(bx,by,bz,max,i/max,0,k/max)
         field(1,i,k) = bx    ! vector components of the field
         field(2,i,k) = by
         field(3,i,k) = bz
      next k
   next i
end sub
! calculate the field components here
! bx,by,bz are the field components at site x,y,z
sub calculate_field(bx,by,bz,max,x,y,z)
   radius = 0.5    ! radius of circular current loop
   dtheta = 2 * pi / 20    ! step size for integration along theta
   bx = 0
   by = 0
   bz = 0
   for theta = 0 to 2*pi-dtheta step dtheta
      dlx = -radius * dtheta * sin(theta)    ! components of dr
      dly = radius * dtheta * cos(theta)      
      rx = x - radius*cos(theta)             ! components of r
      ry = y - radius*sin(theta)
      rz = z 
      r = sqr(rx^2 + ry^2 + rz^2)
      if r > 1 / max then                    ! avoid the wire itself
         bx = bx + dly * rz / r^3            ! Biot-Savart law (cross product)
         by = by - dlx * rz / r^3
         bz = bz + (dlx * ry - dly * rx) / r^3
      end if
   next theta
end sub
! display the results
sub display(field(,,),max)
   dim x(0),z(0),bx(0),by(0),bz(0)  ! arrays used for plotting
   max2 = (2 * max + 1)^2
   mat redim x(max2),z(max2),bx(max2),by(max2),bz(max2)
   m = 0
   if max > 5 then     ! choose a pleasing number of elements to display
      n = 2
   else
      n = 1
   end if
   for i = -max to max step n
      for k = -max to max step n
         m = m + 1
         x(m) = i/max
         z(m) = k /max
         bx(m) = field(1,i,k)
         bz(m) = field(3,i,k)
      next k
   next i
   mat redim x(m),z(m),bx(m),bz(m)
   call settitle("Magnetic field from a current loop")
   call sethlabel("X")
   call setvlabel("Z")
   call vectorgraph(x,z,bx,bz,.010,3,"black")
end sub