Phase Contras Transfer Function

This relatively simple notebook allows to plot the phase contrast transfer function. It only takes defocus, spherical aberration and 2-fold astigmatism into account.

It is written in python.

You can edit is and play with the sliders.

As for the envelopes, only the temporal cohereny is used, by using the energy spread and the chromatic aberration coefficient.

#import matplotlib.pyplot as plt
import numpy as np


from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
#from IPython.display import HTML

from bokeh.io import push_notebook, show, output_notebook
from bokeh.plotting import figure
output_notebook()

Loading BokehJS ...

# universal constants
h=6.626*10**(-34) #plank's constant J.s
me=9.109*10**(-31) #electron's mass kg
mecsquared=511 # electron's rest mass in keV
electron=1.602*10**(-19)# electron's charge

#electron's wavelength in m for an energy in keV
def lamb(Energy):
    l=h/np.sqrt(2*me*Energy*1000*electron)*(1+Energy/(2*mecsquared))**(-0.5)
    return l

#check (should be 2.5 pm)
lamb(200)

2.5080989936058438e-12

# phase contrast transfer function 
# u reciprocal length
# deltaf: defocus
# cs:spherical aberration
# lamb: wavelength
def pctf(u,deltaf,cs,lamb):
    chi=np.pi*deltaf*lamb*u**2+0.5*np.pi*cs*lamb**3*u**4
    tf=np.sin(chi)
    return tf

# phase contrast transfer function squared (to compare with diffractogram)

def pctfsq(u,deltaf,cs,lamb):
    chi=np.pi*deltaf*lamb*u**2+0.5*np.pi*cs*lamb**3*u**4
    tf=(np.sin(chi))**2
    return tf

# 2-d PCTF. Only interesting whn astigmatism is present
def pctf2d(ux,uy,deltaf,astig,cs,lamb):
    chi=np.pi*(deltaf+astig/2)*lamb*(ux**2)+np.pi*(deltaf-astig/2)*lamb*(uy**2)+0.5*np.pi*cs*lamb**3*(ux**2+uy**2)**2
    tf=np.sin(chi)
    return tf

Temporal envelope

As defined in the book by John Spence

$$E(u)=\exp({-\frac{(\pi\Delta\lambda u^2)^2}{16 \ln(2)}})$$ with $$\Delta = C_c\frac{\Delta E}{E_0}$$ $\Delta E$ energy spread

$E_0$, incident energy

# temporal envelope, as defined in Spence book
def temporal_envelope(Energy,u,cc,Espread):
    lam=lamb(Energy)
    delta=cc*Espread/(1000*Energy)
    te=np.exp(-(np.pi*delta*lam*u**2)**2/(16*np.log(2)))
    #te=np.exp(-(np.pi*delta*lam*u**2)**2/(4))
    return te

Scherzer's defocus: $$\Delta f_S=-\sqrt{4*Cs*\lambda/3}$$

def scherzer(cs,Energy):
    lam=lamb(Energy)
    sch=-(4*cs*10**(-6)*lam/3)**(0.5)
    return sch

# Check. Should be -63 nm
scherzer(1200,200)

-6.334791543349591e-08

# Interactive plot of the phase contrast transfer function


u = np.linspace(0,10,2000)
Energy0=200 #kV
deltaf0=-63 #nm
cs0=1200 #microns
cc0=1.2 #mm
Espread0=1.6

t =temporal_envelope(Energy0,10**9*u,cc0*10**(-3),Espread0)*pctf(10**9*u,deltaf0*10**(-9),cs0*10**(-6),lamb(Energy0))
e = temporal_envelope(Energy0,10**9*u,cc0*10**(-3),Espread0)

p = figure(title="Phase contrast transfer function", height=500, width=950,y_range=(-1.1,1.1)) 
p.xaxis.axis_label = 'u (1/nm)'
p.yaxis.axis_label = 'PCTF'


r1 = p.line(u,t,legend_label="PCTF", line_width=2.5)
r2 = p.line(u,e,legend_label="Envelope",line_width=2.5,color="orange")
#show(p)

def update(deltaf=-63, Energy=200, cs=1200,cc=1,Espread=1):
    r1.data_source.data['y']=temporal_envelope(Energy,10**9*u,cc*10**(-3),Espread)*\
                            pctf(10**9*u,deltaf*10**(-9),cs*10**(-6),lamb(Energy))
    r2.data_source.data['y']=temporal_envelope(Energy,10**9*u,cc*10**(-3),Espread)
    push_notebook()
    
p.legend.click_policy="hide"    
show(p, notebook_handle=True)

interact(update, \
         deltaf=widgets.IntSlider(min=-100,max=100,value=deltaf0,description="defocus (nm)"), \
         Energy=widgets.IntSlider(min=60,max=300,value=Energy0,description="Energy (kV)"), \
         cs=widgets.IntSlider(min=-100,max=1500, value=cs0,description="Cs (microns)"), \
         cc=widgets.FloatSlider(min=0,max=3, value=cc0,description="Cc(mm)"), \
         Espread=widgets.FloatSlider(min=0,max=8,value=Espread0,description="E spread (eV)"));

interactive(children=(IntSlider(value=-63, description='defocus (nm)', min=-100), IntSlider(value=200, descrip…

# Interactive plot of the PCTF squared

u = np.linspace(0,10,2000)
Energy0=200 #kV
deltaf0=-63 #nm
cs0=1200 #microns
cc0=1.2 #mm
Espread0=1.6

t =temporal_envelope(Energy0,10**9*u,cc0*10**(-3),Espread0)*pctfsq(10**9*u,deltaf0*10**(-9),cs0*10**(-6),lamb(Energy0))
e = temporal_envelope(Energy0,10**9*u,cc0*10**(-3),Espread0)

p = figure(title="Phase contrast transfer function", height=500, width=1000,y_range=(-0.1,1.1)) 
p.xaxis.axis_label = 'u (1/nm)'
p.yaxis.axis_label = 'PCTF'


r1 = p.line(u,t,legend_label="PCTF", line_width=2.5)
r2 = p.line(u,e,legend_label="Envelope",line_width=2.5,color="orange")
#show(p)

def update(deltaf=-63, Energy=200, cs=1200,cc=1,Espread=1):
    r1.data_source.data['y']=temporal_envelope(Energy,10**9*u,cc*10**(-3),Espread)*\
                            pctfsq(10**9*u,deltaf*10**(-9),cs*10**(-6),lamb(Energy))
    r2.data_source.data['y']=temporal_envelope(Energy,10**9*u,cc*10**(-3),Espread)
    push_notebook()
    
p.legend.click_policy="hide"    
show(p, notebook_handle=True)

interact(update, \
         deltaf=widgets.IntSlider(min=-100,max=100,value=-63,description="defocus (nm)"), \
         Energy=widgets.IntSlider(min=60,max=300,value=200,description="Energy (kV)"), \
         cs=widgets.IntSlider(min=-100,max=1500, value=1200,description="Cs (microns)"), \
         cc=widgets.FloatSlider(min=0,max=3, value=1.2,description="Cc(mm)"), \
         Espread=widgets.FloatSlider(min=0,max=8,value=1.6,description="E spread (eV)"));

interactive(children=(IntSlider(value=-63, description='defocus (nm)', min=-100), IntSlider(value=200, descrip…

import matplotlib.cm  
from matplotlib.colors  import rgb2hex

#mycmap = get_cmap('seismic', 256)    # PiYG
mycmap = matplotlib.colormaps.get_cmap('seismic')

mycolor=[]

for i in range(mycmap.N):
    rgb = mycmap(i)[:3] # will return rgba, we take only first 3 so we get rgb
    mycolor.append(rgb2hex(rgb))
    

#from pylab import *
from bokeh.models import LinearColorMapper, ColorBar
from bokeh.models import CustomJS, Slider
from bokeh.models import ColumnDataSource,FixedTicker
from bokeh.models import CustomJS, Slider
from bokeh.layouts import row, column

ux = np.linspace(-7, 7, 1001)
uy = np.linspace(-7, 7, 1001)
uxx, uyy = np.meshgrid(ux, uy)
uu = np.sqrt(uxx**2 + uyy**2)

Energy0=200 #kV
deltaf0=-63 #nm
cs0=1200 #microns
astig0=0
cc0=1.2 #mm
Espread0=1.6 # eV

t = pctf2d(uxx*10**9,uyy*10**9,deltaf0*10**(-9),astig0*10**(-9),cs0*10**(-6), lamb(Energy0))*\
    temporal_envelope(Energy0,10**9*uu,cc0*10**(-3),Espread0)

p = figure(x_range=(-10,10),y_range=(-10,10),height=600, width=600)
color_mapper = LinearColorMapper(palette="RdBu11", low=-1, high=1)
i1 = p.image(image=[t], x=-10, y=-10, dw=20, dh=20,  palette=mycolor)
color_bar = ColorBar(color_mapper=color_mapper, label_standoff=12, border_line_color=None, location=(0,0))

def update(deltaf=-63, astig=0, cs=1200,cc=1.2, Energy=200,Espread=1.6):
    i1.data_source.data['image'] = \
    [pctf2d(uxx*10**9,uyy*10**9,deltaf*10**(-9),astig*10**(-9),cs*10**(-6), lamb(Energy))*\
    temporal_envelope(Energy,10**9*uu,cc*10**(-3),Espread)]
    push_notebook()

show(p, notebook_handle=True)

interact(update, \
         deltaf=widgets.IntSlider(min=-100,max=100,value=-63,description="defocus (nm)"), \
         Energy=widgets.IntSlider(min=60,max=300,value=200,description="Energy (kV)"), \
         cs=widgets.IntSlider(min=-100,max=1500, value=1200,description="Cs (microns)"), \
         cc=widgets.FloatSlider(min=0,max=3, value=1.2,description="Cc(mm)"), \
         Espread=widgets.FloatSlider(min=0,max=8,value=1.6,description="E spread (eV)"),\
         astig=widgets.FloatSlider(min=0,max=100,value=0,description="Astigmatism (microns)"));


#interact(update, deltaf0=(-100,100), astig0 = (0,100), Energy0 = (60,300), cs = (-100,1500));

interactive(children=(IntSlider(value=-63, description='defocus (nm)', min=-100), FloatSlider(value=0.0, descr…

Notebook by Cécile Hébert (2020-2023). Except where otherwise noted, the content of this notebook is licensed under MIT licence.