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 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. 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) 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 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])) 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) 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]))] # 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 # As Python expression: printing.init_printing(pretty_print=False, use_latex=False) E_x, E_y 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 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 # 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