A charged particle moving in a magnetic field will experience a magnetic force governed by Lorentz force law. This force is always perpendicular to the direction of motion of the particle. If the particle moves in uniform magnetic field and the velocity of the particle is normal to the direction of the field, the particle end up moving in a circle. Now, if the direction of the field is not normal to the velocity, the particle starts moving in a helix, unless its motion is completely parallel to the field in which the there is no Lorentz force and the particle moves in a perfectly straight line. Additionally, with a non-uniform field, you may be able to trap the helical moving particle in a so called magnetic mirror. This may happen if the particle is traveling from a region of weak magnetic field into a region with strong magnetic field.
In this example, we will construct a non-uniform magnetic field and trap the motion of the charged particle in a magnetic mirror. If you are familiar with our notebook Uniform magnetic field, you may skip some of the theory section, as its the same in both notebooks.
A particle with charge $Q$ and mass $m_p$ moving with a velocity $\mathbf{v}$ in a magnetic field $\mathbf{B}$ experiences a force
\begin{equation} \mathbf{F}_{mag} = Q(\mathbf{v} \times \mathbf{B}) \label{Lorentz} \end{equation}known as the Lorentz force [1].
The Lorentz force can be written as a second order ordinary differential equation (ODE)
\begin{equation} m_p\frac{d^2\mathbf{x}}{dt^2} = Q\big(\frac{d\mathbf{x}}{dt} \times \mathbf{B}\big), \label{LorentzODE} \end{equation}where the velocity components are written as the time derivative of the positional components $v_i = dx_i/dt$. By decomposing the vectors, \eqref{LorentzODE} can be reduced into three sets of two first order ODE's each.
\begin{equation} \begin{aligned} &\frac{dx}{dt} = v_x\\[4pt] &\frac{dv_x}{dt} = \frac{Q}{m_{p}}(v_yB_z - v_zB_y)\\[10pt] &\frac{dy}{dt} = v_y\\[4pt] &\frac{dv_y}{dt}= \frac{Q}{m_{p}}(v_zB_x - v_xB_z)\\[10pt] &\frac{dz}{dt} = v_z\\[4pt] &\frac{dv_z}{dt} = \frac{Q}{m_{p}}(v_xB_y - v_yB_x) . \label{LorentzComponents} \end{aligned} \end{equation}Using a numerical method such as the Euler method, Runge-Kutta fourth order, or an embedded Runge-Kutta pair, these equations can be numerically integrated using a specified time step, $\Delta t$, and the motion of a charged particle in a magnetic field can be found.
To to trap a charged particle in a mirror we want a magnetic field which mainly points, and whose magnitude varies, in the $z$-direction. A quadratic function of the direction of the field, plus a radial component is typical [3]
\begin{equation} \mathbf{B}(x, y, z) = B_r\hat{\mathbf{r}} + B_0 \big(1 + \frac{z^2}{L^2} \big)\hat{\mathbf{z}} \; . \label{field} \end{equation}Here the mangetic field is on cylindrical form and $L$ is a parameter controlling the mirror length in $z$-direction. This field is axisymmetric (cylindrical symmetry), i.e. the azimuthal component $B_\theta = 0$ and $\partial/\partial\theta = 0$. From Maxwell's equations $\nabla \cdot \mathbf{B} = 0$ we derive (using the nabla operator corresponding to cylindrical coordinate frames)
\begin{equation} \begin{aligned} \nabla \cdot \mathbf{B} &= \frac{1}{r}\frac{\partial(rB_r)}{\partial r} + \frac{1}{r}\frac{\partial B_\theta}{\partial \theta} + \frac{\partial B_z}{\partial z} \\[4pt] &= \frac{1}{r}\frac{\partial(rB_r)}{\partial r} + \frac{\partial B_z}{\partial z} = 0 \end{aligned} \end{equation}Carrying out the partial derivation of $B_z$ we get
\begin{equation} \begin{aligned} \frac{1}{r}\frac{\partial(rB_r)}{\partial r} &= -2B_0\frac{z}{L^2} \\[4pt] rB_r &= -\int_0^r 2B_0\frac{z}{L^2} r dr \\[4pt] \implies B_r &= -rB_0\frac{z}{L^2} , \end{aligned} \end{equation}leading to the magnetic field
\begin{equation} \mathbf{B}(x, y, z) = -xB_0\frac{z}{L^2}\hat{\mathbf{x}} -yB_0\frac{z}{L^2}\hat{\mathbf{y}} + B_0 \big(1 + \frac{z^2}{L^2} \big)\hat{\mathbf{z}} , \label{BCartesian} \end{equation}on cartesian form.
The usual definition for the magnetic moment of a current loop with a current $I$ and area $A$ is [1]
$$ \mu = IA.$$An electron moving around in a circle of radius $r$, completing a loop with a frequency $f = \omega/2\pi$, generates a current $I = e\omega/2\pi$. Here $\omega$ is the angular frequency of the movement. With the area $A = \pi r^2 = \pi v^2\omega^2$, where $v$ is the velocity of the electron, the magnetic moment can be written as
\begin{equation} \mu = \frac{1}{2}\frac{v^2e}{\omega} \label{magneticMoment} \end{equation}As stated in the introduction, a charged particle will move in a circle if it has a velocity component perpendicular to an external magnetic field. This is called Larmor gyration, and the gyroradius (or Larmor radius), the radius of the circle, is given by [2]
\begin{equation} r_L = \frac{mv_{\perp}}{|Q| B} , \label{larmorRadius} \end{equation}where $v_{\perp}$ is the velocity of the particle perpendicular to the magnetic field with magnitude $B$.
import numpy as np
import matplotlib.pyplot as plt
import sys
from mpl_toolkits.mplot3d import Axes3D
plt.style.use("bmh")
%matplotlib inline
import progressbar
import sys
import time
figsize = (6, 6)
dpi = 600
To begin with, we plot the magnetic field itself. To plot the streamlines of the field in the different planes, a function which returns the magnetic field vector at a given point is needed.
def magneticMirror(x, y, z, B_0, L):
"""Calculates and returns the magnetic field components
given the cartesian position of the particle and B_0 and L."""
# Component field strength in cartesian coordinates
B_x = -x*B_0*z/L**2
B_y = -y*B_0*z/L**2
B_z = B_0*(1 + z**2/L**2)
# Magnetic vector field matrix
B = np.array([B_x, B_y, B_z])
return B
Using a meshgrid, this function will return the magnetic vector field in any plane of interest.
B_0 = 30 # The magnetic field magnitude at origin
L = 0.4 # Mirror length parameter
a = 2*L # To be used for grid creation
# Define meshgrid limits
[xmin, xmax, ymin, ymax, zmin, zmax] = [-a, a, -a, a, -a, a]
# Grid resolution
n = 100j
# 2D meshgrid for creating streamplot
y, x = np.mgrid[xmin:xmax:n, ymin:ymax:n]
def plotMagneticFieldInPlane(subplot, plane, xlabel, ylabel, abscissa, ordinate, B_a, B_o):
"""Function for streamplotting (in a specified subplot)the magnetic
field in a given plane with given grid and field components.
Parameters:
subplot int; the position of the subplot
plane string; which plane is plottet? to be used in
in the title.
xlabel string; label of the subplot's abscissa.
ylabel string; label of the subplot's ordinate.
abscissa nxn array; the "x-values" of the mesh grid points.
ordinate nxn array; the "y-values" of the mesh grid points.
B_a nxn array; the magnetic field components in the direction
of the abscissa.
B_o nxn array; the magnetic field components in the direction
of the ordinate."""
plt.subplot(subplot)
plt.title("Magnetic field in the %s-plane"%(plane))
# Adding linewidth
# Thicker streamline represent higher magnetic field magnitude
magnitude = np.sqrt(B_a**2 + B_o**2) # Absolute magnitude of the magnetic field
lw = 5*magnitude/np.amax(magnitude)
plt.streamplot(abscissa, ordinate, B_a, B_o, linewidth=lw)
plt.plot(0, 0, 'o', color="blue")
plt.xlabel(xlabel)
plt.ylabel(ylabel)
plt.figure(figsize=(12,10), dpi=200)
# xy-plane
B = magneticMirror(x, y, x*0, B_0, L)
plotMagneticFieldInPlane(221, 'xy', r"$x$", r"$y$", x, y, B[0], B[1])
# yz-plane
# In case of unequal axis limits, re-create meshgrid
z, y = np.mgrid[ymin:ymax:n, zmin:zmax:n]
B = magneticMirror(y*0, y, z, B_0, L)
plotMagneticFieldInPlane(222, 'yz', r"$y$", r"$z$", y, z, B[1], B[2])
# xz-plane
z, x = np.mgrid[xmin:xmax:n, zmin:zmax:n]
B = magneticMirror(x, x*0, z, B_0, L)
plotMagneticFieldInPlane(223, 'xz', r"$x$", r"$z$", x, z, B[0], B[2])
# xy-plane with z=1
y, x = np.mgrid[xmin:xmax:n, ymin:ymax:n]
z = np.ones(np.shape(x))
B = magneticMirror(x, y, z, B_0, L)
plotMagneticFieldInPlane(224, 'xy', r"$x$", r"$y$", x, y, B[0], B[1])
plt.show()
/Users/eilifso/anaconda3/envs/numfys/lib/python3.7/site-packages/ipykernel_launcher.py:34: RuntimeWarning: invalid value encountered in true_divide
As expected, there is no magnetic field component in the $xy$-plane at $z=0$. From the other three subplots, it's obvious that the magnetic field magnitude increases away from the origin, as wanted. To plot the streamlines in 3D, the function mlab.quiver3d
from the package mayavi
is used.
from mayavi import mlab
# Enabling 3D plot from mlab in jupyter notebook
mlab.init_notebook()
Notebook initialized with x3d backend.
# 3D meshgrid for calculating magnetic vector field
Z, Y, X = np.mgrid[xmin:xmax:n, ymin:ymax:n, zmin:zmax:n]
# Get magnetic field at each grid point
B = magneticMirror(X, Y, Z, B_0, L)
# Get values relative to max field strength
B = B/np.amax(B)
Bx, By, Bz = B
# Initialise figure
fig = mlab.figure(size=(1000,1000))
# Plot quivers (arrows)
field = mlab.quiver3d(X, Y, Z, Bx, By, Bz)
mlab.axes(figure=fig, x_axis_visibility=True, xlabel="x", y_axis_visibility=True,
ylabel="y", z_axis_visibility=True, zlabel="z")
# Changing background color and viewing angle
fig.scene.background = (0,0,0)
fig.scene.y_minus_view()
fig.scene.camera.roll(90)
fig.scene.camera.zoom(1.5)
# The field strength far from origin is relatively very strong.
# Thus, we mask out those, and showing only points with a strength
# within the specified interval
vectors = fig.children[0].children[0].children[0]
vectors.glyph.glyph.range = np.array([0, .7])
# Setting max amount of data points
vectors.glyph.mask_points.maximum_number_of_points = 3000
vectors.glyph.mask_input_points = True
# Displaying the figure using `fig` (does not work for nbviewer)
#fig
# Save the 3D plot
mlab.savefig('particleTrajectory.jpg' )
# We will here display a saved screenshot
Since the Lorentz force is a cross product between the particle's velocity and the magnetic field, its direction is always perpendicular to the velocity of the particle. Consequently, the Lorentz force does no work on the particle, and the particle's kinetic energy is conserved. A test for numerical validity is then to check the deviation in initial and final absolute velocity.
The velocity component perpendicular to the magnetic field, $v_{\perp}$, will cause a Larmor gyration of radius $r_L$ that creates a magnetic moment $\mu$. Realising that for circular motion $\omega = v/r$ we get (from \eqref{magneticMoment})
$$ \mu = \frac{1}{2}\frac{v^2e}{\omega} = \frac{ve}{2}r. $$Substituting $v \rightarrow v_{\perp}$, $r \rightarrow r_L$ and $e \rightarrow Q$, we can insert the expression for Larmor radius (equation \eqref{larmorRadius})
\begin{equation} \mu = \frac{ve}{2}r \rightarrow \frac{v_\perp Q}{2}r_L = \frac{v_{\perp}|Q|}{2} \frac{mv_{\perp}}{|Q| B} = \frac{1}{2}\frac{mv_{\perp}^2}{B} \label{mu} \end{equation}This is the magnetic moment of a gyrating particle. As the particle moves into a region of stronger magnetic field, its perpendicular velocity will increase accordingly, such that the magnetic moment stays invariant
$$ \frac{d\mu}{dt} = 0.$$See for [2] proof.
Hence, the non-invariance of the particle's magnetic moment can be used as an indication of numerical error.
As discussed in the uniform field section, magnetic fields do no work, and the total energy is invariant. As $v_\perp$ increases for stronger magnetic fields, the velocity component parallel to the magnetic field, $v_\parallel$, must decrease, and if the magnetic field reaches a certain magnitude, $v_\parallel$ becomes zero and the particle is reflected towards the weak-field region.
The midplane of the magnetic mirror is (using \eqref{BCartesian}) where the magnetic field $B=B_0$. At this plane $v_\perp = v_{\perp 0}$ and $v_\parallel = v_{\parallel 0}$. Let the field be $B_{max}$ at the reflection point. There $v_\parallel = 0$ and $v_\perp = v_\perp^\prime$. Conservation of magnetic moment yields
\begin{equation} \begin{aligned} \frac{1}{2}\frac{mv_{\perp 0}^2}{B_0} &= \frac{1}{2}\frac{mv_{\perp}^{\prime 2}}{B_{max}}\\[4pt] \frac{B_0}{B_{max}} &= \frac{v_{\perp 0}^2}{v_{\perp}^{\prime 2}}. \end{aligned} \end{equation}The total energy is invariant at both points
$$ v_0^2 = v_{\perp 0}^2 + v_{\parallel 0}^2 = v_\perp^{\prime 2} $$so that
\begin{equation} \frac{B_0}{B_{max}} = \sin^2\theta \label{lossCone} \end{equation}where $\theta$ is the angle between $\mathbf{v}_0$ and $\mathbf{v}_{\parallel 0}$. If the particle never reaches a region of field magnitude $B_{max}$, given a $B_0$ and $\theta$, it will not be reflected and will escape the mirror. Thus, for a given $B_0$ and $B_{max}$, a magnetic mirror has a loss cone defined by \eqref{lossCone}. On the other hand, the magnetic field equation used in this notebook diverges with $z$ so that every possible particle is reflected at some point, as long as they have a magnetic moment (i.e. $v_{\perp 0}$ is nonzero).
By initialising the particle at the midplane, we set $\theta$. If $B_0$ is known, one can find $B_{max}$ using \eqref{lossCone}. Setting \eqref{BCartesian} equal to $B_{max}$ yields the mirror length
\begin{equation} z = \pm L \sqrt{\frac{B_{max}}{B_0} -1}. \label{mirrorLength} \end{equation}The radial component of the magnetic field generates a force in the azimuthal direction. This is easily seen from writing the Lorentz force in cylindrical components. Remember, $B_\theta$ is zero.
\begin{equation} \begin{aligned} F_r &= Qv_\theta B_z \\[4pt] F_\theta &= Q(v_z B_r - v_r B_z) \\[4pt] F_z &= Qv_\theta B_r \label{lorentzForceCylindrical} \end{aligned} \end{equation}$F_r$ and the second term in $F_\theta$ give rise to the Larmor gyration. $F_z$ causes the particle to slow down and reflect. The last azimuthal force term causes a drift in the radial direction. In essence, it causes the the guiding centre to follow the lines of field [1]. Thus, if the particle is initialised at the Larmor radius $r_L$ such that the initial radial force points towards the $z$-axis, then the guiding centre is at the z-axis, at which the field is parallel, and the guiding centre will not drift.
In this non-uniform magnetic field we place a charged particle. The particle has an initial velocity component parallel and perpendicular to the magnetic field in its initial position. We solve the Lorentz force law to calculate the trajectory. The magnetic field has to be calculated using the function magneticMirror
during each integration step. To avoid diverging error, the adaptive Runge-Kutta method is used instead of the Euler method which is used in Uniform magnetic field. At the point of reflection $v_\perp$ is large relative to the velocity components when the particle is at $z=0$. The error in a integration step increases with increasing particle velocity and acceleration. Thus, to minimise computation time while maintaining an acceptable truncation error, it would be ideal to have a short time step close to the mirror edge, and a longer time step in the middle. An embedded Runge-Kutta pair will estimate the the local truncation error and adjust the time step accordingly, almost without extra computation time. In this problem we will implement a second (Heun's method) and third order embedded Runge-Kutta pair.
We start with an initial time step $h_0$, and will in each integration step $i$ calculate the difference between the approximation by Heun's method and the third order Runge-Kutta method, denoted as $e_i$. If $e_i$ is larger than a given tolerance, $TOL$, we will create a new, ideal, reduced time step, $\tilde{h}$, according to the formula
\begin{equation} \tilde{h} = \left(\frac{TOL |w_i|}{e_i}\right)^{1/p+1}h_i, \label{newTimeStep} \end{equation}where $w_i$ is the approximation of the lower order Runge-Kutta method (in this case Heun's method), $p$ is the order of this method, and $h_i$ is the current time step. For the derivation of this expression and other relevant embedded Runge-Kutta pairs, check out the Adaptive Runge-Kutta notebook. If $e_i$ is too small we will double the time step to reduce computation time. Morover, we will in this notebook use a safety factor of $0.8$ when calculating the ideal time step, such that $h_{i+1} = 0.8 \tilde{h}$. We will also calculate the relative error as
$$ e_{\text{rel}} = \frac{e_i}{\text{max}(w_i, \theta)}, $$where $\theta > 0$, to avoid dividing by a tiny $w$.
To implement the new method, we define a function for the right hand side of the equations in \eqref{LorentzComponents}, an embedded Runge-Kutta pair integration step function, and a function that loops the integration steps and calculates new time steps.
def RHS(M, RHSargs):
""" Returns the right hand side of the sets of ordinary
differential equations.
Parameters:
M 1x6 array containing the position and
velocity of the particle in cartesian
coordinates.
RHSargs 1x5 array containing the arguments/variables
of the right hand side of the Lorentz force.
Returns:
The time derivative of the position and velocity.
"""
x, y, z, vx, vy, vz = M
Q, mass, Bx, By, Bz = RHSargs
ddt_vx = Q/mass*(vy*Bz - vz*By)
ddt_vy = Q/mass*(vz*Bx - vx*Bz)
ddt_vz = Q/mass*(vx*By - vy*Bx)
ddt_x = vx
ddt_y = vy
ddt_z = vz
return ddt_x, ddt_y, ddt_z, ddt_vx, ddt_vy, ddt_vz
def ode23(M, dt, RHS, RHSargs):
""" Performs one step of a second order Runge-Kutta
method (Heun's method), and one step of a third
order Runge-Kutta method RK3. Both steps are returned.
Parameters:
M 1x6 array containing the position and
velocity of the particle in cartesian
coordinates.
dt Integration time step.
RHS Function that returns the right hand side of the
first order ODE's, i.e. their time derivatives.
RHSargs 1x5 array containing the arguments/variables
of the right hand side of the Lorentz force.
Returns:
Heun Next integration step (position and velocity)
calculated using Heun's method.
RK3 Next integration step (position and velocity)
calculated using a third order Runge-Kutta method.
"""
k1 = np.array(RHS(M, RHSargs))
k2 = np.array(RHS(M + k1*dt, RHSargs))
k3 = np.array(RHS(M + (k1 + k2)*dt/4, RHSargs))
Heun = M + dt/2*(k1 + k2)
RK3 = M + dt/6*(k1 + k2 + 4*k3)
return Heun, RK3
Then, a function for calculating particle propagation based on initial conditions is defined.
def adaptivePair(M, mu, times, RHSargs, B_0, L, timeStep, odePair, TOL, tmax, theta):
""" Integrates the particle position in the magnetic
mirror field using an adaptive Runge-Kutta pair. The time
step of the integration step is reduced until the relative
error of the lower powered RK-method is below the tolerance.
Every integration step is stored in M, which contain the initial
velocity components and position. Returns a renewed M and the
magnetic moment at each point.
Parameters:
M 6xNdataPoints array; position and velocity
mu 1xNdataPoints array; magnetic moment
RHSargs 1x5 array; arguments for the RHS of the ODE's
B_0 float; the magnetic field magnitude at origin
L float; mirror length parameter
timeStep float; the initial integration time step
odePair function; ODE solver
TOL float; error tolerance
tmax int/float; the maximum (end) time of the simulation.
theta float; a small number to protect against dividing by
an very small0 lower order approximation (see RKlow
in the function) when calculating the relative error
between the approximations.
returns:
M 6xNdataPoints array; position and velocity (filled)
mu 1xNdataPoints array; magnetic moment (filled)
"""
# Create an 6x1 array of theta to be able to compare with
# approximations in x, y, z, vx, vy, and vz at the same time.
thetas = np.zeros((len(M[:, 0])))
thetas[:] = theta
i = 0 # index keeping track of the iterations
t = 0 # time corresponding to iteration #i
h = timeStep # time step of the integration (dynamic)
x0, y0, z0 = M[0:3, 0] # initial position of the particle
Bx, By, Bz = magneticMirror(x0, y0, z0, B_0, L) # magnetic field components
# corresponding to the initial position
while t < tmax:
x, y, z = M[0:3, i] # Current position
# Store magnetic field components
RHSargs[2], RHSargs[3], RHSargs[4] = Bx, By, Bz
# Calculate next step by calling ode-solver
RKlow, RKhigh = odePair(M[:, i], timeStep, RHS, RHSargs)
# Get difference in position and velocity approximation between the
# two Runge-Kutta methods
error = np.abs(RKhigh - RKlow)
# Get the largest error relative to the lower order approximations
# Using theta if approximation is too small
largestRelError = np.amax(error/np.amax([thetas, RKlow], axis=0))
# In the case the largest relatice error larger than our tolerance
if largestRelError > TOL:
# New reduced time step according to formula
h = 0.8*(TOL*largestRelError)**(1/3)*timeStep
# Get new approximations with new time step
RKlow, RKhigh = odePair(M[:, i], h, RHS, RHSargs)
# Calculate new relative error
largestRelError = np.amax(np.abs(RKhigh - RKlow)/np.amax([thetas, RKlow], axis=0))
# If the error is too big still, half the time step until satisfaction
while largestRelError > TOL:
h = h/2.0 # New smaller time step
RKlow, RKhigh = odePair(M[:, i], h, RHS, RHSargs)
largestRelError = np.amax(np.abs(RKhigh - RKlow)/np.amax([thetas, RKlow], axis=0))
# In the case the relative error is too small, if so, double step size
elif largestRelError < 0.1*TOL:
h = h*2.0
# Assign new time step
timeStep = h
# Insert new accepted position and velocity into M-matrix
M = np.insert(M, i+1, RKhigh, axis=1)
# Extract position, velocity and alculate magnetic field
# to calculate magnetic moment
vx, vy, vz = M[3, i+1], M[4, i+1], M[5, i+1] # velocity
x, y, z = M[0, i+1], M[1, i+1], M[2, i+1] # position
Bx, By, Bz = magneticMirror(x, y, z, B_0, L) # magnetic field in position
B_squared = Bx**2 + By**2 + Bz**2 # Magnetic field magnitude squared
v_squared = vx**2 + vy**2 + vz**2 # Velocity magnitude squared
# Velocity perpendicular to the magnetic field
v_perp_squared = v_squared - (vx*Bx + vy*By + vz*Bz)**2/B_squared
# Calculate magnetic moment using formula
mu = np.insert(mu, i+1, 0.5*mass*v_perp_squared/np.sqrt(B_squared))
# Insert current time into time array
times = np.insert(times, i+1, times[i] + timeStep)
# Increment iterations
i += 1
# increment current time by time step
t += timeStep
print('Number of iterations: %1.0i' %(i))
return M, mu, times
Finally, we initialise the particle, integrate and plot the trajectory.
TOL = 1e-3 # Error tolerance betweeen Runge-Kutta pair
theta = 1e-2 # Constant to prevent to large relative error
timeStep = 1e-3 # Initial time step
simulationDuration = 6
# Setting parameters
Q = 1 # Particle charge
mass = 1 # Particle mass
B_0 = 30 # Magnetic field parameter
L = 0.4 # Magnetic field parameter
RHSargs = np.array([Q, mass, 0.0, 0.0, 0.0])
vx0 = 0 # Initial velocity components
vy0 = 1
vz0 = 2
v_perp_0 = np.sqrt(vx0**2 +vy0**2)
r_L = mass*v_perp_0/(Q*B_0) # The Larmor radius given vx0 and vy0
# Initial v_perp in positive y-direction. Right hand rule
# reveals radial force in positive x-direction. Hence,
# initial negative x-coordinate to avoid guiding centre drift
x0 = -r_L
y0 = 0
z0 = 0
mu_0 = 0.5*mass*v_perp_0**2/B_0 # The initial magnetic moment
# Initialising arrays
M = np.zeros((6, 1))
mu = np.zeros((1))
times = np.zeros((1))
mu[0] = mu_0
M[:, 0] = np.array([x0, y0, z0, vx0, vy0, vz0])
times[0] = 0
M, mu, times = adaptivePair(M, mu, times, RHSargs, B_0, L, timeStep, ode23, TOL, simulationDuration, theta)
# Calculating deviation in total energy and magnetic moment
startVel = np.absolute(np.linalg.norm(M[3:6, 0], axis=0)) # Initial velocity magnitude
endVel = np.absolute(np.linalg.norm(M[3:6, -1], axis=0)) # Final velocity magnitude
velocityDeviation = np.abs(endVel - startVel)/startVel*100 # Magnitude deviation in percent
print("Initial to final velocity deviation:\t\t%.3f percent"%(velocityDeviation))
muDeviation = np.abs(mu[0] - mu[-1])/mu[0]*100 # Deviation of initial and final
# magnetic moment in percent
print("Initial to final magnetic moment deviation:\t%.3f percent"%(muDeviation))
# Plotting x and y position of the particle
plt.figure(figsize=(10,5), dpi = 600)
plt.subplot(311)
plt.plot(times, M[0])
plt.ylabel("x, m")
plt.xlabel('Time, s')
plt.title(r"$x$-position of particle")
plt.subplot(312)
plt.plot(times, M[1])
plt.ylabel("y, m")
plt.xlabel('Time, s')
plt.title(r"$y$-position of particle")
plt.subplots_adjust(hspace=1)
# Plotting the z-position of the particle
#times = np.arange(0, simulationDuration, timeStep)
plt.figure(figsize=(10,3.4), dpi = 600)
plt.plot(times, M[2])
plt.ylabel("z, m")
plt.xlabel('Time, s')
plt.title(r"$z$-position of particle")
plt.show()
# Evaluating analytical mirror length
B_max = B_0*(vx0**2 + vy0**2 + vz0**2)/v_perp_0
mirrorLength = L*np.sqrt(B_max/B_0 - 1)
print("Analytical mirror length:\t%1.2f m"%(mirrorLength))
Number of iterations: 20408 Initial to final velocity deviation: 0.280 percent Initial to final magnetic moment deviation: 2.113 percent
Analytical mirror length: 0.80 m
From the result we see that the particle is mirrored five times during the simulation, and the analytical mirror length matches very good with the numerical result. In addition, both the magnetic moment and the total energy remained fairly constant. The guiding centre remains along the $z$-axis, seen from the pattern of the $x$ and $y$ coordinates.
# Changing parameters to change mirror characteristics
L = 2 # Extending L
x0 = 0.1 # Changing start position away from Larmor Radius
simulationDuration *= 2
# Re-initialising arrays
M = np.zeros((6, 1))
mu = np.zeros((1))
times = np.zeros((1))
mu[0] = mu_0
M[:, 0] = np.array([x0, y0, z0, vx0, vy0, vz0])
times[0] = 0
M, mu, times = adaptivePair(M, mu, times, RHSargs, B_0, L, timeStep, ode23, TOL, simulationDuration, theta)
# Plotting the cartesian coordinate of the particle
plt.figure(figsize=(10,9), dpi = 600)
# Time evolution of x
plt.subplot(5, 2, (1, 2))
plt.plot(times, M[0])
plt.ylabel(r"$x$, m")
plt.xlabel('Time, s')
plt.title(r"$x$-position of particle")
# Time evolution of y
plt.subplot(5, 2, (3, 4))
plt.plot(times, M[1])
plt.ylabel(r"$y$, m")
plt.xlabel('Time, s')
plt.title(r"$y$-position of particle")
# Evaluating analytical mirror length
B_max = B_0*(vx0**2 + vy0**2 + vz0**2)/v_perp_0
mirrorLength = L*np.sqrt(B_max/B_0 - 1)
# Time evolution of z
plt.subplot(5, 2, (5, 6))
plt.plot(times, M[2])
plt.ylabel(r"$z$, m")
plt.xlabel('Time, s')
plt.title(r"$z$-position of particle")
plt.yticks((mirrorLength, -mirrorLength))
# Trajectory seen from xy-plane
plt.subplot(5, 2, (7, 9))
plt.plot(M[0], M[1])
plt.xlabel(r"$x$, m")
plt.ylabel(r"$y$, m")
plt.title(r"Particle motion in $xy$-plane")
# Trajectory seen from xz-plane
plt.subplot(5, 2, (8, 10))
plt.plot(M[0], M[2])
plt.xlabel(r"$x$, m")
plt.ylabel(r"$z$, m")
plt.title(r"Particle motion in $xz$-plane")
plt.subplots_adjust(hspace=1.2)
plt.show()
print("Analytical mirror length:\t%1.2f m"%(mirrorLength))
Number of iterations: 41567
Analytical mirror length: 4.00 m
As expected, the mirror length is proportional to $L$. In both examples we have plottet until now, the mirror length has been $2L$. This is because the loss cone $\theta$ in equation \eqref{lossCone} have remained the same, giving a constant $B_{max}/B_0 = 5$. Try changing the initial velocity components to get $L \neq 2z$.
We see also from the upper and two lower subplots that the guiding centre follows an arc, such as in the field-plots, when the particle is initialised at origo.
The particle's trajectory is plotted in 3D using mayavi
below.
L = 0.8
# Re-initialising arrays
M = np.zeros((6, 1))
mu = np.zeros((1))
times = np.zeros((1))
mu[0] = mu_0
M[:, 0] = np.array([x0, y0, z0, vx0, vy0, vz0])
times[0] = 0
M, mu, times = adaptivePair(M, mu, times, RHSargs, B_0, L, timeStep, ode23, TOL, simulationDuration, theta)
# Calculating deviation in total energy and magnetic moment
startVel = np.absolute(np.linalg.norm(M[3:6, 0], axis=0)) # Initial velocity magnitude
endVel = np.absolute(np.linalg.norm(M[3:6, -1], axis=0)) # Final velocity magnitude
velocityDeviation = np.abs(endVel - startVel)/startVel*100 # Magnitude deviation in percent
print("Initial to final velocity deviation:\t\t%.3f percent"%(velocityDeviation))
muDeviation = np.abs(mu[0] - mu[-1])/mu[0]*100 # Deviation of initial and final
# magnetic moment in percent
print("Initial to final magnetic moment deviation:\t%.3f percent"%(muDeviation))
x, y, z = M[0], M[1], M[2]
fig = mlab.figure()
field = mlab.plot3d(x, y, z, line_width=0.1, tube_radius=0.001)
# Changing background color and viewing angle
fig.scene.background = (0,0,0)
fig.scene.y_minus_view()
# Displaying the figure using `fig` (does not work for nbviewer)
#fig
# Save the 3D plot
mlab.savefig('particleTrajectory.jpg' )
# We will here display a saved screenshot
Number of iterations: 38167 Initial to final velocity deviation: 0.569 percent Initial to final magnetic moment deviation: 4.787 percent
[1]: Chen, F. 1984, Introduction to plasma physics and controlled fusion. Volume 1: Plasma physics (New York: Plenum Press)
[2]: Bittencourt, J. 2004, Fundamentals of plasma physics (New York: Springer-Verlag)
[3]: Ripperda, B. et al. 2019, A comprehensive comparison of relativistic particle integrators. ApJS.