%matplotlib inline
import numpy
import matplotlib.pyplot as mpl
from mpl_toolkits.mplot3d import Axes3D
# Parameters for plot attributes
mpl.rc("xtick", labelsize="large")
mpl.rc("ytick", labelsize="large")
mpl.rc("axes", labelsize="xx-large")
mpl.rc("axes", titlesize="xx-large")
mpl.rc("figure", figsize=(8,8))
# define key constants
m_p = 1.67E-27 # mass of proton: kg
qe = 1.602E-19 # charge of proton: C
# now, setting up an alpha particle
m = 4.0*m_p
q = 2.0*qe
QoverM = q/m
dt = 1.0E-8
t = numpy.arange(0.0, 0.002, dt)
rp = numpy.zeros((len(t), 3))
vp = numpy.zeros((len(t), 3))
v0 = numpy.sqrt(2.0*QoverM*10.0)
# negative ion
mn = 1.0 * m_p
qn = -1.0*qe
QoverMn = qn/mn
rn = numpy.zeros((len(t), 3))
vn = numpy.zeros((len(t), 3))
# Strength of Magnetic Field
B0 = 1.0E-4
expected_R = m*v0/(q*B0)
expected_T = 2.0*numpy.pi*expected_R / v0
print("v0 = ", v0)
print("Expected trajectory radius = ", expected_R)
print("Expected cyclotron period = ", expected_T)
# initial condition
rp[0,:] = numpy.array([0.0, 0.0, 0.0])
vp[0,:] = numpy.array([v0, 0.0, 0.0])
# initial condition for negative ion
rn[0,:] = numpy.array([0.0, 0.0, 0.0])
vn[0,:] = numpy.array([v0, 0.0, 0.0])
# Euler time steps
for it in range(0, len(t)-1):
Bp = numpy.array([0.0, 0.0, B0])
Ap = QoverM * numpy.cross(vp[it,:], Bp)
vp[it+1] = vp[it] + dt*Ap
rp[it+1] = rp[it] + dt*vp[it]
An = QoverMn * numpy.cross(vn[it,:], Bp)
vn[it+1] = vn[it] + dt*An
rn[it+1] = rn[it] + dt*vn[it]
# Plot the particle's trajectory, in the xy-plane
mpl.plot(rp[:,0], rp[:,1], label="Alpha particle")
mpl.plot(rn[:,0], rn[:,1], linestyle = "--", label="H$^-$ ion")
mpl.xlabel("$x$")
mpl.ylabel("$y$")
mpl.title("Particle Trajectories in Uniform Magnetic Field")
mpl.legend(fontsize="x-large")
mpl.xlim(-13,13)
mpl.ylim(-13,13)