Chapter 2: Diffraction
5.4. Homework 4#
Analyzing Ring Diffraction Pattern
part of
MSE672: Introduction to Transmission Electron Microscopy
by Gerd Duscher, Spring 2025
Microscopy Facilities
Institute of Advanced Materials & Manufacturing
Materials Science & Engineering
The University of Tennessee, Knoxville
Background and methods to analysis and quantification of data acquired with transmission electron microscopes.
5.4.1. Overview#
This homework follows the notebook: Analyzing Ring Diffraction Pattern
5.4.2. Load relevant python packages#
5.4.2.1. Check Installed Packages#
import sys
from pkg_resources import get_distribution, DistributionNotFound
def test_package(package_name):
"""Test if package exists and returns version or -1"""
try:
version = get_distribution(package_name).version
except (DistributionNotFound, ImportError) as err:
version = '-1'
return version
if test_package('pyTEMlib') < '0.2025.1.0':
print('installing pyTEMlib')
!{sys.executable} -m pip install --upgrade pyTEMlib
print('done')
5.4.2.2. Load the plotting and figure packages#
Import the python packages that we will use:
Beside the basic numerical (numpy) and plotting (pylab of matplotlib) libraries,
three dimensional plotting and some libraries from the book
kinematic scattering library.
%matplotlib widget
import matplotlib.pyplot as plt
import numpy as np
# 3D and selection plotting package
from mpl_toolkits.mplot3d import Axes3D # 3D plotting
from matplotlib.widgets import EllipseSelector
# additional package
import itertools
import scipy.constants as const
import os
import sys
if 'google.colab' in sys.modules:
from google.colab import output
output.enable_custom_widget_manager()
from google.colab import drive
# Import libraries from the book
import pyTEMlib
import pyTEMlib.kinematic_scattering as ks # Kinematic sCattering Library
# with Atomic form factors from Kirklands book
import pyTEMlib.file_tools as ft
# it is a good idea to show the version numbers at this point for archiving reasons.
__notebook_version__ = '2025.02.10'
print('pyTEM version: ', pyTEMlib.__version__)
print('notebook version: ', __notebook_version__)
You don't have igor2 installed. If you wish to open igor files, you will need to install it (pip install igor2) before attempting.
You don't have gwyfile installed. If you wish to open .gwy files, you will need to install it (pip install gwyfile) before attempting.
Symmetry functions of spglib enabled
Using kinematic_scattering library version {_version_ } by G.Duscher
pyTEM version: 0.2024.09.0
notebook version: 2025.02.10
5.4.3. Load Ring-Diffraction Pattern#
5.4.3.1. First we select the diffraction pattern#
In the second lab we used a sample of polycrystalline Aluminium.
If you run this notebook on your own computer you should download your images from the google drive for 2025 Lab Data, if you run it on google colab you can go to the drive directory in the dialog below.
You must log into Google with your UTK account to be able to read these data.
Go to the folder of your data in the Lab02 folder and select one
if 'google.colab' in sys.modules:
drive.mount("/content/drive")
fileWidget = ft.FileWidget(sum_frames=True)
diff_pattern = fileWidget.selected_dataset
print(f"alpha tilt {np.degrees(diff_pattern.metadata['experiment']['stage']['tilt']['alpha']):.2f}°")
print(f"beta tilt {np.degrees(diff_pattern.metadata['experiment']['stage']['tilt']['beta']):.2f}°")
view = diff_pattern.plot()
diff_pattern
alpha tilt -37.37°
beta tilt 17.90°
|
||||||||||||||||
diff_pattern.view_metadata()
1/(ks.get_wavelength(200000)/10)
experiment :
detector : BM-Ceta
acceleration_voltage : 200000.0
microscope : Titan
start_date_time : 1740155856
probe_mode : convergent
convergence_angle : 0.006196301826995508
stage :
holder :
position :
x : -0.00012349714500000015
y : -0.00020536990199999992
z : -7.549785541351843e-05
tilt :
alpha : -0.6522094449999999
beta : 0.3124780058860779
398.7345690135052
1/(ks.get_wavelength(200000)/10)
398.7345690135052
5.4.4. Finding the center#
5.4.4.1. Select the center yourself#
Select the center of the screen with the ellipse selection tool
Note: we use the logarithm to plot the diffraction pattern (look for : “np.log” in the code cell below, the number that follows is the gamma value, change it)
## Access the data of the loaded image
#diff_pattern = np.array(main_dataset.sum(axis=0))
diff_pattern = diff_pattern-diff_pattern.min()
radius = diff_pattern.shape[1]/4
center = np.array([diff_pattern.shape[0]/2, diff_pattern.shape[1]/2])
center= np.unravel_index(np.argmax(np.array(diff_pattern), axis=None), diff_pattern.shape)
plt.figure(figsize=(8, 6))
plt.imshow(np.log(3.+diff_pattern).T, origin = 'upper')
current_axis = plt.gca()
selector = EllipseSelector(current_axis,
None,
interactive=True ,
minspanx=5, minspany=5,
spancoords='pixels') # gca get current axis (plot)
center = np.array(center)
selector.extents = (center[0]-radius,center[0]+radius,center[1]-radius,center[1]+radius)
plt.show()
Get center coordinates from selection
xmin, xmax, ymin, ymax = selector.extents
x_center, y_center = selector.center
x_shift = x_center - diff_pattern.shape[0]/2
y_shift = y_center - diff_pattern.shape[1]/2
print(f'radius = {(xmax-xmin)/2:.0f} pixels')
center = (x_center, y_center )
print(f'new center = {center} [pixels]')
out_tags ={}
out_tags['center'] = center
5.4.5. Ploting Diffraction Pattern in Polar Coordinates#
5.4.5.1. The Transformation Routine#
from scipy.interpolate import interp1d
from scipy.ndimage import map_coordinates
def cartesian2polar(x, y, grid, r, t, order=3):
R,T = np.meshgrid(r, t)
new_x = R*np.cos(T)
new_y = R*np.sin(T)
ix = interp1d(x, np.arange(len(x)))
iy = interp1d(y, np.arange(len(y)))
new_ix = ix(new_x.ravel())
new_iy = iy(new_y.ravel())
return map_coordinates(grid, np.array([new_ix, new_iy]),
order=order).reshape(new_x.shape)
def warp(diff,center):
# Define original polar grid
nx = diff.shape[0]
ny = diff.shape[1]
x = np.linspace(1, nx, nx, endpoint = True)-center[0]
y = np.linspace(1, ny, ny, endpoint = True)-center[1]
z = diff
# Define new polar grid
nr = int(min([center[0], center[1], diff.shape[0]-center[0], diff.shape[1]-center[1]])-1)
print(nr)
nt = 360*3
r = np.linspace(1, nr, nr)
t = np.linspace(0., np.pi, nt, endpoint = False)
return cartesian2polar(x,y, z, r, t, order=3).T
5.4.5.2. Now we transform#
If the center is correct a ring in carthesian coordinates is a line in polar coordinates
A simple sum over all angles gives us then the diffraction profile (intensity profile of diffraction pattern)
center = np.array(center)
out_tags={'center': center}
polar_projection = warp(diff_pattern,center)
below_zero = polar_projection<0.
polar_projection[below_zero]=0.
out_tags['polar_projection'] = polar_projection
# Sum over all angles (axis 1)
profile = polar_projection.sum(axis=1)
out_tags['radial_average'] = profile
scale = ft.get_slope(diff_pattern.dim_0.values)
plt.figure()
im = plt.imshow(np.log(1000+polar_projection),extent=(0,360,polar_projection.shape[0]*scale,scale),cmap="gray")# , vmin=np.max(np.log2(1+diff_pattern))*0.5)
plt.colorbar(im)
ax = plt.gca()
ax.set_aspect("auto");
plt.xlabel('angle [degree]');
plt.ylabel('distance [1/nm]')
plt.plot(profile/profile.max()*200,np.linspace(1,len(profile),len(profile))*scale,c='r');
In the image above check:
Are the lines straight?
5.4.6. Determine Bragg Peak#
Peak finding is actually not as simple as it looks
# --- Input ------
scale = 1.
# ----------------
import scipy as sp
import scipy.signal as signal
# find_Bragg peaks in profile
peaks, g= signal.find_peaks(profile,rel_height =1.1, width=7) # np.std(second_deriv)*9)
print('Peaks are at pixels:')
print(peaks)
out_tags['ring_radii_px'] = peaks
plt.figure()
plt.imshow(np.log2(1.+polar_projection),extent=(0,360,polar_projection.shape[0]*scale,scale),cmap='gray', vmin=np.max(np.log2(1+diff_pattern))*0.5)
ax = plt.gca()
ax.set_aspect("auto");
plt.xlabel('angle [degree]');
plt.ylabel('distance [1/nm]')
plt.plot(profile/profile.max()*200,np.linspace(1,len(profile),len(profile))*scale,c='r');
for i in peaks:
if i*scale > 3.5:
plt.plot((0,360),(i*scale,i*scale), linestyle='--', c = 'steelblue')
5.4.7. Calculate Ring Pattern#
Note that you will need to change the material
see Structure Factors notebook for details.
# -------Input -----
material = 'gold'
# -------------------
# Initialize the dictionary with all the input
atoms = ks.structure_by_name(material)
#ft.h5_add_crystal_structure(main_dataset.h5_dataset.file, atoms)
#Reciprocal Lattice
# We use the linear algebra package of numpy to invert the unit_cell \"matrix\"
reciprocal_unit_cell = atoms.cell.reciprocal() # transposed of inverted unit_cell
#INPUT
hkl_max = 7# maximum allowed Miller index
acceleration_voltage = 200.0 *1000.0 #V
wave_length = ks.get_wavelength(acceleration_voltage)
h = np.linspace(-hkl_max,hkl_max,2*hkl_max+1) # all to be evaluated single Miller Index
hkl = np.array(list(itertools.product(h,h,h) )) # all to be evaluated Miller indices
g_hkl = np.dot(hkl,reciprocal_unit_cell)
# Calculate Structure Factors
structure_factors = []
base = atoms.positions # in Carthesian coordinates
for j in range(len(g_hkl)):
F = 0
for b in range(len(base)):
f = ks.feq(atoms[b].symbol,np.linalg.norm(g_hkl[j])) # Atomic form factor for element and momentum change (g vector)
F += f * np.exp(-2*np.pi*1j*(g_hkl[j]*base[b]).sum())
structure_factors.append(F)
F = structure_factors = np.array(structure_factors)
# Allowed reflections have a non zero structure factor F (with a bit of numerical error)
allowed = np.absolute(structure_factors) > 0.001
distances = np.linalg.norm(g_hkl, axis = 1)
print(f' Of the evaluated {hkl.shape[0]} Miller indices {allowed.sum()} are allowed. ')
# We select now all the
zero = distances == 0.
allowed = np.logical_and(allowed,np.logical_not(zero))
F = F[allowed]
g_hkl = g_hkl[allowed]
hkl = hkl[allowed]
distances = distances[allowed]
sorted_allowed = np.argsort(distances)
distances = distances[sorted_allowed]
hkl = hkl[sorted_allowed]
F = F[sorted_allowed]
# How many have unique distances and what is their muliplicity
unique, indices = np.unique(distances, return_index=True)
print(f' Of the {allowed.sum()} allowed Bragg reflections there are {len(unique)} families of reflections.')
intensity = np.absolute(F[indices]**2*(np.roll(indices,-1)-indices))
print('\n index \t hkl \t 1/d [1/Ang] d [pm] F multip. intensity' )
family = []
#out_tags['reflections'] = {}
reflection = 0
for j in range(len(unique)-1):
i = indices[j]
i2 = indices[j+1]
family.append(hkl[i+np.argmax(hkl[i:i2].sum(axis=1))])
index = '{'+f'{family[j][0]:.0f} {family[j][1]:.0f} {family[j][2]:.0f}'+'}'
print(f'{i:3g}\t {index} \t {distances[i]:.4f} \t {1/distances[i]*100:.0f} \t {np.absolute(F[i]):4.2f} \t {indices[j+1]-indices[j]:3g} \t {intensity[j]:.2f}')
#out_tags['reflections'+str(reflection)]={}
out_tags['reflections-'+str(reflection)+'-index'] = index
out_tags['reflections-'+str(reflection)+'-recip_distances'] = distances[i]
out_tags['reflections-'+str(reflection)+'-structure_factor'] = np.absolute(F[i])
out_tags['reflections-'+str(reflection)+'-multiplicity'] = indices[j+1]-indices[j]
out_tags['reflections-'+str(reflection)+'-intensity'] = intensity[j]
reflection +=1
5.4.8. Comparison#
Comparison between experimental profile and kinematic theory
The grain size will have an influence on the width of the diffraction rings”
# -------Input of grain size ----
first_peak_pixel = 100
first_peak_reciprocal_distance = 0.4247
pixel_size = first_peak_reciprocal_distance/first_peak_pixel
resolution = 0 # 1/nm
thickness = 100 # Ang
# -------------------------------
print(f'Pixel size is {pixel_size:.5f} 1/Ang')
from scipy import signal
width = (1/thickness + resolution) / scale
# scale = ft.get_slope(main_dataset.dim_0.values) *1.085*1.0/10
scale = pixel_size
intensity2 = intensity/intensity.max()*10
gauss = signal.windows.gaussian(len(profile), std=width)
simulated_profile = np.zeros(len(profile))
rec_dist = np.linspace(1,len(profile),len(profile))*pixel_size
plt.figure()
plt.plot(rec_dist,profile/profile.max()*150, color='blue', label='experiment');
for j in range(len(unique)-1):
if unique[j] < len(profile)*scale:
# plot lines
plt.plot([unique[j],unique[j]], [0, intensity2[j]],c='r')
# plot indices
index = '{'+f'{family[j][0]:.0f} {family[j][1]:.0f} {family[j][2]:.0f}'+'}' # pretty index string
plt.text(unique[j],-3, index, horizontalalignment='center',
verticalalignment='top', rotation = 'vertical', fontsize=8, color = 'red')
# place Gaussian with appropriate width in profile
g = np.roll(gauss,int(-len(profile)/2+unique[j]/scale))* intensity2[j]*10#rec_dist**2*10
simulated_profile = simulated_profile + g
plt.plot(np.linspace(1,len(profile),len(profile))*scale,simulated_profile/50, label='simulated');
plt.xlabel('angle (1/$\AA$)')
plt.legend()
plt.ylim(-.5,10)
5.4.9. Publication Quality Output#
Now we have all the ingredients to make a publication quality plot of the data.
from matplotlib import patches
plot_profile = profile.copy()
plot_profile[:first_peak_pixel-20] = 0
fig = plt.figure(figsize=(9, 6))
extent= np.array([-center[0], diff_pattern.shape[0]-center[0],-diff_pattern.shape[1]+center[1], center[1]])*scale
plt.imshow(np.log(3.+diff_pattern).T,cmap='gray', extent=(extent*1.0)) #, vmin=np.max(np.log2(1+diff_pattern))*0.5)
plt.xlabel(r'reciprocal distance [nm$^{-1}$]')
ax = fig.gca()
#ax.add_artist(circle1);
plt.plot(np.linspace(1,len(profile),len(profile))*scale,plot_profile/plot_profile.max(), color='y');
plt.plot((0,len(profile)*scale),(0,0),c='r')
for j in range(len(unique)-1):
i = indices[j]
if distances[i] < len(profile)*scale:
plt.plot([distances[i],distances[i]], [0, intensity2[j]/20],c='r')
arc = patches.Arc((0,0), distances[i]*2, distances[i]*2, angle=90.0, theta1=0.0, theta2=270.0, color='r', fill= False, alpha = 0.5)#, **kwargs)
ax.add_artist(arc);
plt.scatter(0,0);
for i in range(6):
index = '{'+f'{family[i][0]:.0f} {family[i][1]:.0f} {family[i][2]:.0f}'+'}' # pretty index string
plt.text(unique[i],-0.05, index, horizontalalignment='center',
verticalalignment='top', rotation = 'vertical', fontsize=8, color = 'white')
5.4.10. Homework#
Determine the pixel_size and for two different indicated camera lengths!
Submit one notebook with your diffraction pattern
Optional:
Plot the indicated camera length over the pixel size!