Contents

Answers

0. Some setup

Import numpy for numerics and pyplot for visualisation

import numpy as np
import matplotlib.pyplot as plt

Define some useful constants

echarge = 1.60217662e-19  # electron charge
emass = 9.10938356e-31    # electron mass
vperm = 1.25663706212e-6  # vacumm permittivity

1. Defining some vector products

Define some vectors for testing

v1 = np.array([1.0, 2.0, 3.0])
v2 = np.array([4.0 ,5.0, 6.0])
v3 = np.array([7.0 ,8.0, 9.0])

1.a. Dot product

def dot(v1, v2):
    return (v1[0]*v2[0]) + (v1[1]*v2[1]) + (v1[2]*v2[2]) 

1.b. Testing the dot product

check1 = 'Passed' if (dot(v1, v2) == 32.0) else 'Failed'
print(f'Check 1: {check1}')
check2 = 'Passed' if (dot(v1, v3) == 50.0) else 'Failed'
print(f'Check 2: {check2}')
check3 = 'Passed' if (dot(v2, v3) == 122.0) else 'Failed'
print(f'Check 3: {check3}')
Check 1: Passed
Check 2: Passed
Check 3: Passed

1.c. Cross product

def cross(v1, v2):
    x = (v1[1]*v2[2]) - (v1[2]*v2[1])
    y = (v1[2]*v2[0]) - (v1[0]*v2[2])
    z = (v1[0]*v2[1]) - (v1[1]*v2[0])
    return np.array([x, y, z])

1.d. Testing the cross product

check1 = 'Passed' if (cross(v1, v2) == np.array([-3, 6, -3])).all() else 'Failed'
print(f'Check 1: {check1}')
check2 = 'Passed' if (cross(v1, v3) == np.array([-6, 12, -6])).all() else 'Failed'
print(f'Check 2: {check2}')
check3 = 'Passed' if (cross(v2, v3) == np.array([-3, 6, -3])).all() else 'Failed'
print(f'Check 3: {check3}')
Check 1: Passed
Check 2: Passed
Check 3: Passed

2. Visualising some magnetic fields

2.a. Defining some field functions

def Bfield_constant(pos, t):
    return np.array([0.0, 0.0, 1.0])

def Bfield_spacevarying(pos, t):
    return np.array([0.0, 0.0, 1/np.sqrt(dot(pos, pos))])

def Bfield_timevarying(pos, t):
    return np.array([0.0, 0.0, np.cos(t*(2*np.pi))])

2.b. Defining a function to generate a grid of points parrallel to yz-plane

def gen_grid(x, y0, z0, ymax, zmax, dyz):
    num_ys = round((ymax-y0)/dyz)
    ymax = y0 + num_ys*dyz  # adjust ymax if (ymax - ymin) incompatible with dyz
    num_zs = round((zmax-z0)/dyz)
    zmax = z0 + num_zs*dyz  # adjust ymax if (zmax - zmin) incompatible with dyz
    points = []
    for i in range(num_ys):
        for j in range(num_zs):
            points.append(np.array([x,
                                    y0 + (dyz/2.0) + (i*dyz), 
                                    z0 + (dyz/2.0) + (j*dyz)]))
            
    return points

Making a few plots to test gen_grid

def plot_grid(xyz):
    x, y, z = zip(*tuple(xyz))
    plt.figure(figsize=(9,6))
    plt.scatter(y, z, s=10, c='b') # s for size, c for color
    plt.xlabel(r'$y$', fontsize=20)
    plt.ylabel(r'$z$', fontsize=20)
    plt.show()
    
grid = gen_grid(0, -2, -2, 2, 2, 0.2)
plot_grid(grid)
grid = gen_grid(0, 0, 0, 5, 10, 0.2)
plot_grid(grid)
_images/mag_23_0.png _images/mag_23_1.png

2.c. Visualizing the magnetic fields

def plot_Bfield(Bfield_function, xyz, t, title=None):
    x, y, z = zip(*tuple(xyz))
    Bvecs = [Bfield_function(pos, t) for pos in xyz]
    Bx, By, Bz = zip(*tuple(Bvecs))
    plt.figure(figsize=(9,6))
    plt.quiver(y, z, By, Bz, scale=40.0)
    plt.xlabel(r'$y$', fontsize=20)
    plt.ylabel(r'$z$', fontsize=20)
    if not title is None:
        plt.title(title, fontsize=25)
    plt.show()
    
grid = gen_grid(0, -2, -2, 2, 2, 0.2)
plot_Bfield(Bfield_constant, grid, 0, 'constant')
plot_Bfield(Bfield_spacevarying, grid, 0, 'space-varying')
plot_Bfield(Bfield_timevarying, grid, 0.0, 'time-varying, $t=0.0$')
plot_Bfield(Bfield_timevarying, grid, 0.25, 'time-varying, $t=0.25$')
plot_Bfield(Bfield_timevarying, grid, 0.5, 'time-varying, $t=0.5$')
_images/mag_25_0.png _images/mag_25_1.png _images/mag_25_2.png _images/mag_25_3.png _images/mag_25_4.png

3. Writing and testing the integrator for the vector differential equation

3.a. Writing the integrator

def euler(vec_init, t_init, t_final, dt, diff_func, *args, **kwargs):
    del_t = t_final - t_init # total duration to solve for
    n_dt = round(del_t/dt)
    dt = del_t/n_dt # adjust dt if (t_final - t_init) incompatible with dt
    
    t = t_init
    vec = vec_init
    
    ts, vecs = [t], [vec] # initialise the output lists with the initial values as first element
    for i in range(n_dt):
        dvec = diff_func(vec, t, *args, **kwargs) # get the change in the vector at this step
        vec = (vec + (dvec*dt)) # update vec according to the change calculated above
        t += dt
        
        ts.append(t)
        vecs.append(vec) # add the new values to the output lists
        
    return ts, vecs

3.b. Testing with 1st order differential equation

def diff_func(vec, t, const):  # simple first order case
    return const*vec

const = 2.0
t_init = 0.0
t_final = 2.0
dt = 0.01
vec_init = np.array([1.0])  # boundary condition - the function takes value 1.0 at t_init

ts, vecs = euler(vec_init, t_init, t_final, dt, diff_func, const=const)

y, = zip(*tuple(vecs)) # restructure output for plotting

### plotting ###
plt.figure(figsize=(9,6))
plt.plot(ts, y, c='b', label='numeric')
plt.plot(ts, np.exp(const*np.array(ts)), c='r', label='exact')
plt.xlabel('$t$', fontsize=20)
plt.ylabel('$exp(kt)$', fontsize=20)
plt.title('First order, $k=2.0$', fontsize=25)
plt.legend(loc='best', fontsize=15)
plt.show()
_images/mag_30_0.png

3.c. Testing with 2nd order differential equation

def diff_func(vec, t, const):  # 2nd order case
    vec_length = len(vec)
    return np.concatenate((vec[vec_length//2:], (const**2)*vec[:vec_length//2])) # put function value and 1st derivative together in a single vector

const = 2.0
t_init = 0.0
t_final = 2.0
dt = 0.01
vec_init = np.array([1.0, const])  # boundary condition - the function and 1st derivative take values 1.0 and k at t_init respectively

ts, vecs = euler(vec_init, t_init, t_final, dt, diff_func, const=const)

y, dy = zip(*tuple(vecs)) # restructure output for plotting

### plotting ###
plt.figure(figsize=(9,6))
plt.plot(ts, y, c='b', label='numeric')
plt.plot(ts, np.exp(const*np.array(ts)), c='r', label='exact')
plt.xlabel('$t$', fontsize=20)
plt.ylabel('$exp(kt)$', fontsize=20)
plt.title('2nd order, $k=2.0$', fontsize=25)
plt.legend(loc='best', fontsize=15)
plt.show()
_images/mag_32_0.png

4. Simulate a charged test particle moving in a magnetic field

Representing the particle

def get_particle():
    return {'charge' : echarge,
            'mass' : emass,
            'pos' : np.array([0.0, 0.0, 0.0]),
            'vel' : np.array([1.0e11, 0.0, 0.0])}

4.a. Writing the differential equation

def diff_func(vec, t, mass, charge, Bfield_function):
    vec_length = len(vec)
    pos = vec[:vec_length//2]
    vel = vec[vec_length//2:]
    dpos = vel
    dvel = (charge/mass)*cross(vel, Bfield_function(pos, t))
    return np.concatenate((dpos, dvel))

4.b. and 4.c. Performing some simulations

Define a simulation function

def simulation(Bfield_function, particle, T, dt, integrator, diff_func):
    vec = np.concatenate((particle['pos'], particle['vel']))
    ts, vecs = integrator(vec, 0, T, dt, diff_func, 
                          charge=particle['charge'], 
                          mass=particle['mass'], 
                          Bfield_function=Bfield_function)
    
    pos_xs, pos_ys, pos_zs, vel_xs, vel_ys, vel_zs = zip(*tuple(vecs))
    poss = [np.array(triple) for triple in zip(*(pos_xs, pos_ys, pos_zs))] # restructure integrator output in convenient way
    vels = [np.array(triple) for triple in zip(*(vel_xs, vel_ys, vel_zs))]
    
    particle['pos'] = poss[-1].copy() # update position of particle using final timestep
    particle['vel'] = vels[-1].copy() # update velocity of particle using final timestep
    
    results = {'t':ts, 'pos':poss, 'vel':vels} # package the results up in a nice dict
    return results

Make a 2nd particle with opposite charge but otherwise identical

def get_particle2():
    return {'charge' : -echarge,
            'mass' : emass,
            'pos' : np.array([0.0, 0.0, 0.0]),
            'vel' : np.array([1.0e11, 0.0, 0.0])}

Run the simulations with the 2 particles

t_unit = emass/(1.0*echarge)
dt = 0.001*t_unit
T = (2.0*np.pi)*t_unit

particle = get_particle()
results = simulation(Bfield_constant, particle, T, dt, euler, diff_func)

pos_x, pos_y, pos_z = zip(*tuple(results['pos'])) # restructure positions for plotting

particle2 = get_particle2()
results2 = simulation(Bfield_constant, particle2, T, dt, euler, diff_func)

pos2_x, pos2_y, pos2_z = zip(*tuple(results2['pos'])) # restructure positions for plotting

### plotting ###
plt.figure(figsize=(9,6))
plt.plot(pos_x, pos_y, c='b', label='positive charge')
plt.plot(pos2_x, pos2_y, c='r', label='negative charge')
plt.xlabel('$x$', fontsize=20)
plt.ylabel('$y$', fontsize=20)
plt.legend(loc='best', fontsize=15)
plt.show()
_images/mag_44_0.png

5. Visualising the magnetic field induced by a moving charge

5.a. Writing the induced field function

Function with particle as an argument

def induced_Bfield(pos, t, particle):
    disp = pos - particle['pos'] # displacement vector from particle to position we're looking at
    return (vperm*particle['charge']/(4*np.pi)) * cross(particle['vel'], disp)/dot(disp, disp)

Create a moving particle moving quickly along x axis

def get_xmoving_particle():
    return {'charge' : echarge,
            'mass' : emass,
            'pos' : np.array([0.1, 0.1, 0.1]),
            'vel' : np.array([1.0e26, 0.0, 0.0,])}

Get the field function in its usual form

xmoving_particle = get_xmoving_particle()
def Bfield_particle(pos, t):
    return induced_Bfield(pos, t, xmoving_particle)

5.b. Visualise the induced field

grid = gen_grid(0, -2, -2, 2, 2, 0.2)
plot_Bfield(Bfield_particle, grid, 0)
_images/mag_54_0.png

6. Gauss’s law for magnetism

6.a. and 6.b. Creating the cubic surface and area element vectors

Defining a function to create ‘oriented’ surfaces

def gen_oriented_grid(axis, orientation, coord1, coord2_0, coord3_0, coord2_max, coord3_max, dcoord):
    grid = gen_grid(coord1, coord2_0, coord3_0, coord2_max, coord3_max, dcoord) # get a grid perpendicular to z axis
    axis_i = {'x':0, 'y':1, 'z':2}[axis.lower()]
    grid = [ point[[axis_i, (axis_i+1)%3, (axis_i+2)%3]] for point in grid ] # swap co-ords to make the grid perpendicular to chosen axis
    
    dA_mag = dcoord**2 # get magnitude of area element
    
    dA = np.array([0.0, 0.0, 0.0])
    dA[axis_i] = np.sign(orientation)*dA_mag # get vector form of area element 

    return grid, dA

Define a function to calculate magnetic flux through a regular grid surface

def calculate_flux(Bfield_function, t, surface_points, surface_dA):
    flux = 0.0
    for point in surface_points:
        flux += dot(Bfield_function(point, t), surface_dA)
    return flux

Create a cube centred at the origin

xpos, dA_xpos = gen_oriented_grid('x', 1.0, 0.5, -0.5, -0.5, 0.5, 0.5, 0.001) # create the face on +ve x-axis
xneg, dA_xneg = gen_oriented_grid('x', -1.0, -0.5, -0.5, -0.5, 0.5, 0.5, 0.001) # create the face on -ve x-axis

ypos, dA_ypos = gen_oriented_grid('y', 1.0, 0.5, -0.5, -0.5, 0.5, 0.5, 0.001) # create the face on +ve y-axis
yneg, dA_yneg = gen_oriented_grid('y', -1.0, -0.5, -0.5, -0.5, 0.5, 0.5, 0.001) # create the face on -ve y-axis

zpos, dA_zpos = gen_oriented_grid('z', 1.0, 0.5, -0.5, -0.5, 0.5, 0.5, 0.001) # create the face on +ve z-axis
zneg, dA_zneg = gen_oriented_grid('z', -1.0, -0.5, -0.5, -0.5, 0.5, 0.5, 0.001) # create the face on -ve z-axis

6.c Calculating fluxes through the cube faces

flux_xpos = calculate_flux(Bfield_particle, 0.0, xpos, dA_xpos)
flux_xneg = calculate_flux(Bfield_particle, 0.0, xneg, dA_xneg)

flux_ypos = calculate_flux(Bfield_particle, 0.0, ypos, dA_ypos)
flux_yneg = calculate_flux(Bfield_particle, 0.0, yneg, dA_yneg)

flux_zpos = calculate_flux(Bfield_particle, 0.0, zpos, dA_zpos)
flux_zneg = calculate_flux(Bfield_particle, 0.0, zneg, dA_zneg)

Printing the results

print(f'The flux through the +ve x face is: {flux_xpos}')
print(f'The flux through the -ve x face is: {flux_xneg}')
flux_x = flux_xpos + flux_xneg
print(f'The total flux through the x faces is: {flux_x}')

print(f'The flux through the +ve y face is: {flux_ypos}')
print(f'The flux through the -ve y face is: {flux_yneg}')
flux_y = flux_ypos + flux_yneg
print(f'The total flux through the y faces is: {flux_y}')

print(f'The flux through the +ve z face is: {flux_zpos}')
print(f'The flux through the -ve z face is: {flux_zneg}')
flux_z = flux_zpos + flux_zneg
print(f'The total flux through the x faces is: {flux_z}')

flux = flux_x + flux_y + flux_z
print(f'The total flux through the cube surface is: {flux}')
The flux through the +ve x face is: 0.0
The flux through the -ve x face is: 0.0
The total flux through the x faces is: 0.0
The flux through the +ve y face is: -2.1294852298939033
The flux through the -ve y face is: -1.85747510863246
The total flux through the y faces is: -3.9869603385263632
The flux through the +ve z face is: 2.1294852298939033
The flux through the -ve z face is: 1.85747510863246
The total flux through the x faces is: 3.9869603385263632
The total flux through the cube surface is: 0.0