#

# 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()