Calculate the standard deviation of the gaussian distribution in defocus caused by chromatic aberration in the STEM

chromatic_aberration_simulation.py

import matplotlib.pyplot as plt
import numpy as np

import sympy as sp
sp.init_printing()

E0, dE, Cc = sp.symbols('E0, dE, C_c, ')
deltaZ_sigma = Cc*dE/E0#*2*sp.sqrt(2*sp.log(2)) 

zlp = 0.9
voltage = 3e5
chrom = 1.6e-3

f = sp.lambdify([Cc, dE, E0], deltaZ_sigma)

sigma = f(chrom, zlp, voltage) / 1e-9

mu = 0
x = np.linspace(-20,20,1000)

def gaussian(x, mu, sigma, normalise = False):
    gauss = np.exp((-(x-mu)**2)/(2*sigma**2))
    if normalise:
        gauss = gauss * 1/(sigma*np.sqrt(2*np.pi)) # normalise
    return gauss

gauss = gaussian(x, mu, sigma)

defocus = np.array([-2*sigma, -sigma, 0, sigma, 2*sigma])
weighting = gaussian(defocus, mu, sigma)
print("Standard deviation of defocus is: {}".format(sigma))

print("Chosen five steps of 1sigma apart:")
print(defocus)
print('Weighted by:')
print(weighting)

fig, ax = plt.subplots()
ax.plot(x, gauss, color='grey', zorder=-1, lw=4)
plt.fill_between(x, gauss, color='lightgrey')

plt.scatter(defocus, weighting,color='green')
plt.xlabel('Temporal distribution in defocus / nm')
plt.ylabel('Gaussian weights')
plt.ylim(0,None)

data_coords = np.array([[0, i] for i in weighting])
figure_coords = ax.transData.transform(data_coords)
axes_coords = ax.transAxes.inverted().transform(figure_coords).T[1]
for x0, ymax in zip(defocus, axes_coords):
    plt.axvline(x0, ymax = ymax, ls='--', color='green')