Re-implementation of MatSci Part III research project, originally done as a MuPad program, now re-done for fun using Python tools.
The project did a simulation of electrons that are emitted when an electron beam hits a p-n junction in a scanning electron microscope.
This file is basically a translation of my original MuPad source code, but using features of IPython:
It isn't complete, but most of the work is done. For reference and explanation you'll need to see my original 2002 report.
import math
import numpy as np
from scipy import integrate
import sympy
# Physical constants
q_e = -1.60e-19 # charge of electron
m_e = 9.11e-31 # mass of electron
# Simulation constants. SI units unless otherwise stated
working_distance = 4.5e-3 # Electron microscope geometry
num_events = 400
num_energy_bins = 20
num_angle_bins = 10
t_finish = 2e-9
y0 = 1e-12 # arbitrary
V_BI_base = 0.7 # Corresponds to band gap for silicon
# parameters that we will vary, but have some defaults here for sake of graphs
work_function_base = 4.52 # electron volts
V_extract_base = 20.0
excitation_point_base = 0.0
Standard results for back-scattered electrons in this configuration:
$$ N(E) \propto \frac{E}{(E+\phi)^{4}} $$$$ N(\alpha) \propto cos(\alpha) $$def make_prob_distributions(work_function):
prob_angle = lambda a: np.cos(a)
prob_energy = lambda en: en/(en + work_function) ** 4
return prob_angle, prob_energy
angle_min = - math.pi / 2
angle_range = math.pi
angle_max = angle_min + angle_range
energy_min = 0
# energy_range and energy_max to be calculated.
Normalise:
$$ n_1 = \int_0^\infty p_{energy}(e) {d}e $$$$ n_2 = \int_{-\pi/2}^{\pi/2} p_{angle}(\alpha) {d}\alpha $$def make_normalised_probability_distribution(prob_angle, prob_energy):
n1, _ = integrate.quad(prob_energy, 0, 10^9)
n2, _ = integrate.quad(prob_angle, -np.pi/2, np.pi/2)
return lambda angle, en: prob_energy(en) * prob_angle(angle) / (n1 * n2)
Check that things look the way they should:
prob_angle_tmp, prob_energy_tmp = make_prob_distributions(work_function_base)
prob_tmp = make_normalised_probability_distribution(prob_angle_tmp, prob_energy_tmp)
%matplotlib inline
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
angle_axis = np.linspace(-np.pi/2, np.pi/2, 50)
energy_axis = np.linspace(0, 20, 100)
X, Y = np.meshgrid(angle_axis, energy_axis)
Z = prob_tmp(X, Y)
fig = plt.figure(figsize=(12, 9))
ax = plt.subplot(111, projection='3d')
ax.view_init(30, 30)
ax.plot_surface(X, Y, Z, cstride=2, rstride=5, cmap=plt.cm.coolwarm)
plt.title("Probability of electron")
ax.set_xlabel("Angle")
ax.set_ylabel("Energy (eV)")
None
We want to adjust the energy range so that we get a self-similar distribution of energies whatever the work function. We scale it so that we have at most 1 electron in the box that is at $e = energyMax, \alpha=0$.
We require zero electrons in the box that would be at $e=energyMax$, $\alpha=0$. This requires:
$$ prob(0, energyMax) . numEvents . \frac{energyMax}{numEnergyBins} . \frac{angleRange}{numAngleBins} < 1/2 $$
since < 0.5 would be rounded down to 0.
Using a symbolic integration on $n_1$ above, it can be shown that the equation we need to solve is:
$$ 6 (energyMax . workFunction)^2 . numEvents . angleMax - numEnergyBins . numAngleBins .(energyMax + workFunction)^4 = 0 $$There are 4 solutions to this, 2 negative and 2 positive, and we need the largest positive one.
energyMax, angleRange, numEnergyBins, numAngleBins, workFunction, numEvents = sympy.symbols("energyMax angleRange numEnergyBins numAngleBins workFunction numEvents")
energyMaxSolutions = sympy.solve(6 * (energyMax * workFunction)**2 * numEvents * angleRange - numEnergyBins * numAngleBins * (energyMax + workFunction)**4, energyMax)
def get_energy_max(work_function):
return float(max([s.evalf(subs=dict(workFunction=work_function,
angleRange=angle_range,
numAngleBins=num_angle_bins,
numEnergyBins=num_energy_bins,
numEvents=num_events))
for s in energyMaxSolutions]))
Build probability distribution:
def make_energy_bins(work_function):
energy_max = get_energy_max(work_function)
prob_angle, prob_energy = make_prob_distributions(work_function)
prob = make_normalised_probability_distribution(prob_angle, prob_energy)
d_e = (energy_max - energy_min)/ num_energy_bins
d_a = angle_range / num_angle_bins
bins = np.empty([num_energy_bins, num_angle_bins])
for i in range(0, num_energy_bins):
for j in range(0, num_angle_bins):
# Really should do double integration over prob(),
# but a good enough numerical approximation is to
# average height at the 4 corners and multiple by base area
v = (prob(angle_min + j * d_a, energy_min + i * d_e) +
prob(angle_min + (j + 1) * d_a, energy_min + i * d_e) +
prob(angle_min + j * d_a, energy_min + (i + 1) * d_e) +
prob(angle_min + (j + 1) * d_a, energy_min + (i + 1) * d_e)) / 4
# Scale by num_events
bins[i, j] = np.round(float(v * d_a * d_e * num_events))
return bins, (angle_min, angle_max), (energy_min, energy_max)
# Plot it:
bins_tmp, _, _ = make_energy_bins(work_function_base)
def hist_3d(matrix):
fig = plt.figure(figsize=(12, 9))
ax = fig.add_subplot(111, projection='3d')
x_data, y_data = np.meshgrid(np.arange(matrix.shape[1]),
np.arange(matrix.shape[0]) )
x_data = x_data.flatten()
y_data = y_data.flatten()
z_data = matrix.flatten()
ax.bar3d(x_data, y_data, np.zeros(len(z_data)),
1, 1, z_data)
return ax
ax = hist_3d(bins_tmp)
ax.view_init(30, 30)
Build a list of angles and energies of electrons according to the distribution above. Each electron is given an energy randomly selected from within the bit it is in.
def generate_electrons(work_function):
np.random.seed(1234) # same each time
r = np.random.random_sample
bins, (angle_min, angle_max), (energy_min, energy_max) = make_energy_bins(work_function)
d_a = (angle_max - angle_min) / num_angle_bins
d_e = (energy_min - energy_max) / num_energy_bins
return [(float(angle_min + (j + r()) * d_a),
float(energy_min + (i + r()) * d_e))
for i in range(0, num_energy_bins) for j in range(0, num_angle_bins) for k in range(0, int(bins[i, j]))]
The form of the potential has been calculated theoretically:
$$ V(x, y) = \frac{V_{BI}}{\pi} arctan(\frac{x}{y}), y > 0 $$The extractor bias $V_{extract}$ and thus the extractor field $E_{extract}$ also needs to be added to this:
$$ V(x,y) = \frac{V_{BI}}{\pi} \arctan \left(\frac{x}{y}\right) - E_{extract}y $$# Symbolic manipulation:
x, y, V_BI, V_extract, workingDistance = sympy.symbols("x y V_BI V_extract workingDistance")
E_extract = -1 * (V_extract / workingDistance)
V = V_BI / sympy.pi * sympy.atan(x/y) - E_extract * y
E_x = -1 * sympy.diff(V, x)
E_y = -1 * sympy.diff(V, y)
E_z = 0
from sympy.interactive import printing
printing.init_printing(use_latex=True)
E_x, E_y
/home/luke/tmpstore/anaconda/lib/python2.7/site-packages/IPython/core/formatters.py:239: FormatterWarning: Exception in image/png formatter: \begin{pmatrix}- \frac{V_{BI}}{\pi y \left(\frac{x^{2}}{y^{2}} + 1\right)}, & \frac{V_{BI} x}{\pi y^{2} \left(\frac{x^{2}}{y^{2}} + 1\right)} - V_{extract} / workingDistance\end{pmatrix} ^ Unknown symbol: \begin (at char 0), (line:1, col:1) FormatterWarning,
# As Python expression:
printing.init_printing(pretty_print=False, use_latex=False)
E_x, E_y
(-V_BI/(pi*y*(x**2/y**2 + 1)), V_BI*x/(pi*y**2*(x**2/y**2 + 1)) - V_extract/workingDistance)
def make_field(V_BI, V_extract):
# Method 1: could use sympy.evalf:
# def E_field(x, y):
# subs = dict(workingDistance=working_distance, V_BI=V_BI, V_extract=V_extract, x=x, y=y)
# return (E_x.evalf(subs=subs),
# E_y.evalf(subs=subs))
# For speed, though, we copy-paste output from above and fix where needed:
def E_field(x, y):
return (-V_BI/(math.pi*y*(x**2/y**2 + 1)),
V_BI*x/(math.pi*y**2*(x**2/y**2 + 1)) - V_extract/working_distance)
return E_field
# Plot it:
E_field_tmp = make_field(V_BI_base, V_extract_base)
E_field_x = lambda x, y: E_field_tmp(x, y)[0]
E_field_y = lambda x, y: E_field_tmp(x, y)[1]
fig = plt.figure(figsize=(15,8))
ax = fig.add_subplot(111, projection='3d')
x = np.linspace(-0.0001, 0.0001, num=50)
y = np.linspace(0.0000001, 0.0002, num=50)
X, Y = np.meshgrid(x, y)
Z = E_field_x(X, Y).clip(-30000, 0)
plt.title("$E_x$")
ax.plot_surface(X, Y, Z, cstride=2, rstride=2, cmap=plt.cm.Reds)
ax.set_xlabel("x (m)")
ax.set_ylabel("y (m)")
None
fig = plt.figure(figsize=(15,10))
ax = fig.add_subplot(111, projection='3d')
x = np.linspace(-4e-10, 4e-10, num=50)
y = np.linspace(1e-11, 5e-10, num=50)
X, Y = np.meshgrid(x, y)
Z = E_field_y(X, Y)
Z = Z.clip(-1.5e9, 1.5e9) # clip to avoid singularity
plt.title("$E_y$")
ax.plot_surface(X, Y, Z, cstride=2, rstride=2, cmap=plt.cm.Greens)
ax.set_xlabel("x (m)")
ax.set_ylabel("y (m)")
None
Equation of motion for electron:
$$ m_e \frac{{d}^2\mathbf{r}}{{d}t^2} = q_e\mathbf{E}(\mathbf{r}) $$where
$$ \mathbf{r} = \left( \begin{array}{l}x\\y\end{array} \right) $$We can convert this second order o.d.e in 2 dependent variables to a first order o.d.e. in 4 equations:
$$ \frac{d}{dt}\left( \begin{array}{l} x\\x'\\y\\y'\\ \end{array} \right) = \left( \begin{array}{l} x'\\\frac{\displaystyle q_e}{\displaystyle m_e}E_x(x,y)\\y'\\\frac{\displaystyle q_e}{\displaystyle m_e}E_y(x,y) \end{array} \right)\\ $$def make_ode(V_BI, V_extract):
E_field = make_field(V_BI, V_extract)
def ode_func(t, vec):
x, y, vx, vy = vec
E_x, E_y = E_field(x, y)
return np.array([vx,
vy,
q_e/m_e * E_x,
q_e/m_e * E_y])
return ode_func
ev_to_joule = lambda energy_ev: abs(q_e) * energy_ev
v_x_0 = lambda angle, energy: np.sqrt(2 * np.abs(energy) / m_e) * np.sin(angle)
v_y_0 = lambda angle, energy: np.sqrt(2 * np.abs(energy) / m_e) * np.cos(angle)
def make_starting_point(electron_angle, electron_energy, excitation_point):
energy = ev_to_joule(electron_energy)
return np.array([excitation_point, # The point where electron beam is hitting the silicon, with p-n junction at origin
y0,
v_x_0(electron_angle, energy),
v_y_0(electron_angle, energy)
])
# Electron micropscope geometry:
DETECTION_HEIGHT = working_distance
DETECTION_WIDTH = 2.5e-3
t_finish_last = 2e-9
# Use a search to try to home in on the point when
# electron will be at the detector height.
def detect(starting_point, ode_func, return_trajectory=False):
global t_finish_last
t_min = 0
num_guesses = 0
t_guess = t_finish_last
t_max = None
t_last_min = None
t_give_up = 0.01
detected = False
while True:
num_guesses += 1
def f(t_final):
r = integrate.ode(ode_func)
r.set_integrator('vode')
r.set_initial_value(starting_point, 0)
return r.integrate(t_final)
x, y, vx, vy = f(t_guess)
if y > DETECTION_HEIGHT:
if y/DETECTION_HEIGHT > 1.1:
# Need to go back
t_max = t_guess
t_guess = (t_min + t_max)/2
else:
# Right height
if abs(x) < DETECTION_WIDTH:
# Correct width - caught
detected = True
#print "FOUND after {0}".format(num_guesses)
else:
detected = False
#print "NOTFOUND after {0}".format(num_guesses)
# Either way, we're done
t_finish_last = t_guess
break
else:
# Need to go higher
if t_guess > t_give_up:
#print "ABANDONED after {0}".format(num_guesses)
break
if t_guess > t_min:
t_min = t_guess
if t_max is None: # We have no idea what to use for t_max
t_guess = t_guess * 2
else:
t_guess = (t_min + t_max) /2
if return_trajectory:
return detected, np.array([f(i) for i in np.arange(0, t_guess, t_guess/20)])
return detected
V_extract_test = 20.0
electrons_test = generate_electrons(work_function_base)
e_angle, e_energy = electrons_test[0]
ode_func_test = make_ode(V_BI_base, V_extract_test)
ic = make_starting_point(e_angle, e_energy, 0)
# t_finish_last = 1e-12
# %timeit detect(ic, ode_func_test)
# Sanity check: plot some trajectories
V_extract_test = 20.0
electrons_test = generate_electrons(work_function_base)
ode_func_test = make_ode(V_BI_base, V_extract_test)
fig = plt.figure(figsize=(10, 12))
for i, (e_angle, e_energy) in enumerate(electrons_test[::4]):
ic = make_starting_point(e_angle, e_energy, 0)
detected, path = detect(ic, ode_func_test, return_trajectory=True)
xs = path[:, 0].clip(-0.004, 0.004)
ys = path[:, 1].clip(0, 0.005)
plt.plot(xs, ys, color="red" if detected else "blue")
# Data structs for memoisation against all the parameters we want.
results_d = {}
electrons_d = {}
def do_sim(excitation_point, work_function, V_extract):
if (excitation_point, work_function, V_extract) in results_d:
return results_d[excitation_point, work_function, V_extract]
if work_function in electrons_d:
electrons = electrons_d[work_function]
else:
electrons = generate_electrons(work_function)
electrons_d[work_function] = electrons
ode_func = make_ode(V_BI_base, V_extract)
c = 0.0
for e_angle, e_energy in electrons:
ic = make_starting_point(e_angle, e_energy, excitation_point)
detected = detect(ic, ode_func)
if detected:
c += 1.0
result = c/float(len(electrons))
results_d[excitation_point, work_function, V_extract] = result
return result
We simulate scanning an electron beam across the p-n junction, using $V_{extract} = 20.0 V$
# More points at centre, where the activity is
excitation_points = np.array([-0.00004, -0.00003, -0.00002, -0.000015, -0.00001, -0.000006, -0.000004,
-0.000002, -0.0000015, -0.000001, -0.0000005, -0.00000025, 0.0, 0.00000025,
0.0000005, 0.000001, 0.0000015, 0.000002, 0.000004, 0.000006, 0.00001,
0.000015, 0.00002, 0.00003, 0.00004,])
# Plot for different values of work function.
fig = plt.figure(figsize=(15,8))
work_functions = [2.0, 3.0, 4.52, 5.5]
V_extract_base = 20.0
for work_function in work_functions:
yields = 100 * np.array([do_sim(excitation_point, work_function, V_extract_base) for excitation_point in excitation_points])
plt.plot(excitation_points, yields)
plt.title("Detection profile for different values of work function")
plt.xlabel("Distance from p-n junction")
plt.ylabel("Percentage detection")
plt.legend(['{:5.2f} ev'.format(w) for w in work_functions])
None