#
# ExB_Filter_Exercise_4.py
#
# This file is used to numerically integrate
# the second order linear differential equations
# that describe the trajectory of a charged particle through
# an E x B velocity filter.
#
# Here, it is assumed that the axis of the filter
# is aligned with the z-axis, that the magnetic field
# is along the +x-direction, and that the electric field
# is along the -y-direction.
#
# The numerical integration is done using the built-in
# routine odeint.
#
# Many particles selected from a normal distribution of
# z-velocities are sent through the filter and histograms
# of the z-velocities of the incident and transmitted particles
# are produced.
#
# By:
# Ernest R. Behringer
# Department of Physics and Astronomy
# Eastern Michigan University
# Ypsilanti, MI 48197
# (734) 487-8799 (Office)
# ebehringe@emich.edu
#
# Last updated:
#
# 20160624 ERB
#
from pylab import figure,xlim,xlabel,ylim,ylabel,grid,title,hist,show,text
from numpy import sqrt,array,arange,random,absolute,zeros,linspace
from scipy.integrate import odeint
#
# Initialize parameter values
#
q = 1.60e-19 # particle charge [C]
m = 7.0*1.67e-27 # particle mass [kg]
KE_eV = 100.0 # particle kinetic energy [eV]
Ex = 0.0 # Ex = electric field in the +x direction [N/C]
Ey = -105.0 # Ey = electric field in the +y direction [N/C]
Ez = 0.0 # Ez = electric field in the +z direction [N/C]
Bx = 0.002 # Bx = magnetic field in the +x direction [T]
By = 0.0 # By = magnetic field in the +x direction [T]
Bz = 0.0 # Bz = magnetic field in the +x direction [T]
R_mm = 0.5 # R = radius of the exit aperture [mm]
L = 0.25 # L = length of the crossed field region [mm]
Ntraj = 40000 # number of trajectories
transmitted_v = zeros(Ntraj) # array to save velocities of transmitted particles
n_transmitted = 0 # counter for the number of transmitted particles
# Derived quantities
qoverm = q/m # charge to mass ratio [C/kg]
KE = KE_eV*1.602e-19 # particle kinetic energy [J]
vmag = sqrt(2.0*KE/m) # particle velocity magnitude [m/s]
R = 0.001*R_mm # aperture radius [m]
vzpass = -Ey/Bx # z-velocity for zero deflection [m/s]
# Set up the distribution of incident velocities
mean = vzpass # the mean of the velocity distribution is vzpass
sigma = 0.1*vzpass # the width of the velocity distribution is 0.1*vzpass
vz = mean + sigma*random.randn(Ntraj) # the set of initial velocity magnitudes
scaled_vz = vz/vzpass # the set of scaled initial velocity magnitudes
# Set up the bins for the histograms
scaled_vz_min = 0.6
scaled_vz_max = 1.4
Nbins = 64
scaled_vz_bins = linspace(scaled_vz_min,scaled_vz_max,Nbins+1)
vz_bins = vzpass*scaled_vz_bins
#
# Over what time interval do we integrate?
#
tmax = L/vzpass;
#
# Specify the time steps at which to report the numerical solution
#
t1 = 0.0 # initial time
t2 = tmax # final scaled time
N = 1000 # number of time steps
h = (t2-t1)/N # time step size
# The array of time values at which to store the solution
tpoints = arange(t1,t2,h)
# Specify initial conditions that don't change
x0 = 0.0 # initial x-coordinate of the charged particle [m]
dxdt0 = 0.0 # initial x-velocity of the charged particle [m/s]
y0 = 0.0 # initial y-coordinate of the charged particle [m]
dydt0 = 0.0 # initial y-velocity of the charged particle [m/s]
z0 = 0.0 # initial z-coordinate of the charged particle [m]
#
# Here are the derivatives of position and velocity
def derivs(r,t):
# derivatives of position components
xp = r[1]
yp = r[3]
zp = r[5]
dx = xp
dy = yp
dz = zp
# derivatives of velocity components
ddx = qoverm*(Ex + yp*Bz - zp*By)
ddy = qoverm*(Ey + zp*Bx - xp*Bz)
ddz = qoverm*(Ez + xp*By - yp*Bx)
return array([dx,ddx,dy,ddy,dz,ddz],float)
# Start the loop over the initial velocities
for i in range (0,Ntraj-1):
# Specify initial conditions
dzdt0 = vz[i] # initial z-velocity of the charged particle [m/s]
r0 = array([x0,dxdt0,y0,dydt0,z0,dzdt0],float)
# Calculate the numerical solution using odeint
r = odeint(derivs,r0,tpoints)
# Extract the 1D matrices of position values
position_x = r[:,0]
position_y = r[:,2]
position_z = r[:,4]
# Extract the 1D matrices of velocity values and final velocity
v_x = r[:,1]
v_y = r[:,3]
v_z = r[:,5]
vxf = v_x[N-1]
vyf = v_y[N-1]
vzf = v_z[N-1]
vf = sqrt(vxf*vxf + vyf*vyf + vzf*vzf)
# If the particle made it through the aperture, save the velocity
if absolute(position_x[N-1]) < R:
if absolute(position_y[N-1]) < sqrt(R*R - position_x[N-1]*position_x[N-1]):
transmitted_v[n_transmitted] = vf
n_transmitted = n_transmitted + 1
# Only save the non-zero values for the histogram
transmitted_v_f = transmitted_v[0:n_transmitted]
scaled_transmitted_v_f = transmitted_v_f/vzpass
# Let the user know how many particles were transmitted
print("The number of incident particles is %d"%Ntraj) #Frem: Added Brackets
print("The number of transmitted particles is %d"%n_transmitted) #Frem: Added Brackets
# start a new figure
figure()
# plot the histogram of scaled initial velocities
n, bins, patches = hist(scaled_vz, scaled_vz_bins, normed=0, facecolor='orange', alpha=0.75)
xlabel('$v_z/v_{z,pass}$ [m/s]',size = 16)
ylabel('$N$',size = 16)
title('Histogram of initial $v_z/v_{z,pass}$ values')
xlim(scaled_vz_min,scaled_vz_max)
ylim(0,0.075*Ntraj)
grid(True)
show()
# start a new figure
figure()
# plot the histogram of scaled final velocities (transmitted particles)
n, bins, patches = hist(scaled_transmitted_v_f, scaled_vz_bins, normed=0, facecolor='purple', alpha=0.75)
xlabel('$v_z/v_{z,pass}$',size = 16)
ylabel('$N$',size = 16)
title('Histogram of $v_z/v_{z,pass}$ values for transmitted particles')
xlim(scaled_vz_min,scaled_vz_max)
ylim(0,0.075*Ntraj)
grid(True)
text(0.65,2750,"R = %.2f mm"%R_mm,size=16)
show()