2711.sagews

The import commands load the libraries we will need to call certain functions.

import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as la

Define some useful constants

mu = [1.0,0.0,0.0] # magnetic dipole moment of bar magnet
kmag = 1e-7

Define a function for computing the field of a particle with a magnetic dipole

def Bdipole(source,obsloc):
    r = np.subtract(source,obsloc)
    return kmag*(3*np.dot(mu,r/la.norm(r))*r/la.norm(r) - mu)/(np.dot(r,r))^(3/2)

Define a function for computing the field of a magnetic monopole

def Bmono(source,obsloc):
    r = np.subtract(source,obsloc)
    return 100*kmag*r/la.norm(r)/(np.dot(r,r))

We will compute the flux at points along a radius of the loop, separated by a distance dR, and use this to compute the integrated flux over the entire loop

center=[0.2,0,0]   #Center of the loop
Rdisk = 0.3        #Radius of the loop
Bsurface = [] 
dR = 0.05*Rdisk     #Determines how many points to sample to compute the flux
for y in np.arange(0,Rdisk,dR):
    pos=[center[0], y, 0]
    Bsurface.append(pos)   #Add point to list of points

Below is the main loop over time.  At each time step, we move the particle forward according to its speed, compute the fields at each point along the radius of the loop, and from this we compute the flux.  By keeping track of the flux from the previous time step in variables called oldflux, we may approximate the flux derivative as the ratio of the change in flux and the change in time.

dt=0.01  #Time step in seconds
xhat = [1,0,0]  #unit vector for the loop surface
vx=0.1  #Speed of the particle in meters per second
dfluxddt=[]   #lists for the flux derivatives
dfluxmdt=[]
magnetpos=[0.0,0.0,0.0]  #initial position of the particle
time=np.arange(0.0,4.0,dt)  #create a list of time steps
fluxd=0.0
fluxm=0.0
oldfluxd=fluxd
oldfluxm=fluxm
for t in time:
    magnetpos[0] += vx*dt  #Advance the particle
    oldfluxd=fluxd   #Store the previous flux values
    oldfluxm=fluxm
    fluxd=0.0
    fluxm=0.0
    #Perform the integral over loop enclosed area
    for arr in Bsurface:    #For each point along the loop radius
        Bd = Bdipole(magnetpos, arr)  #Compute the dipole field
        Bm = Bmono(magnetpos,arr)     #Compute the monopole field
        fluxd+=np.dot(Bd,xhat)*2.0*np.pi*abs(arr[1])*dR  #compute the dipole flux
        fluxm+=np.dot(Bm,xhat)*2.0*np.pi*abs(arr[1])*dR  #Compute the monopole flux
    dfluxddt.append((fluxd-oldfluxd)/dt)  #Compute the flux derivatives
    dfluxmdt.append((fluxm-oldfluxm)/dt)

The following commands make a nice plot of our results.

plt.cla()
plt.clf()
plt.figure()
#[1:-1] means plot the second to last points.  Since we have computed derivatives as differences, the first point has not been computed correctly.
plt.plot(time[1:-1],dfluxddt[1:-1],'k', time[1:-1], dfluxmdt[1:-1],'k--')
plt.xlabel('Time(s)')
plt.ylabel('Current (arbitrary units)')
plt.title('Current induced in a conducting loop due to a particle')
plt.yticks([])
plt.legend(['Magnetic Dipole','Magnetic Monopole'])
ax=plt.gca()
plt.savefig('blah')
plt.show()