#!/usr/bin/python
# -*- coding: utf-8 -*-
"""
The crystal module handles Bragg diffraction on crystal lattices.
"""
__copyright__ = "Copyright (C) 2019 Florian Otte"
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
import re
import math
import time
import xconf.misc.units as units
import xconf.misc.constants as c
[docs]def get_lattice_constant(material_name):
"""
Convenience method for assigning lattice constant parameters of common
analyzer materials.
"""
if material_name == "Si":
return units.Length(c.lattice_constant_si, 'angstrom')
if material_name == "Ge":
return units.Length(c.lattice_constant_ge, 'angstrom')
else:
raise ValueError("Unknown crystal material requested.")
[docs]def diffraction_allowed(hkl, max_diffraction_order):
"""
Returns a list of miller indizes for allowed diffractions from *hkl* up to
the *max_diffraction_order* diffraction order.
"""
i = []
order = range(1, max_diffraction_order+1)
h0, k0, l0 = hkl
for o in order:
h,k,l = o*h0, o*k0, o*l0
if h % 2 == 1 and k % 2 == 1 and l % 2 == 1:
i.append((h,k,l))
elif h % 2 == 0 and k % 2 == 0 and l % 2 == 0 and (h+k+l) % 4 == 0:
i.append((h,k,l))
else:
continue
return i
[docs]def get_cut(hkl):
"""Divides a set of miller indices *hkl* by its greatest common factor."""
def gcf(x, y):
while(y):
x, y = y, x % y
return x
h,k,l = hkl
factor = gcf(gcf(h,k), l)
return tuple([int(i / factor) for i in hkl])
[docs]class Crystal(object):
"""
Create an analyzer crystal by specifying the cut via the Miller indices *h*,
*k* and *l*, the radius of curvature *roc* and a material string. Valid
material strings are 'Si' or 'Ge'.
"""
def __init__(self, h, k, l, roc, material, **kwargs):
self.hkl = get_cut([h,k,l])
self.roc = roc
self.material_string = material
self.lattice_constant = get_lattice_constant(material)
[docs] def get_bragg_angle(self, energy, hkl):
"""
Get the Bragg angle under which the diffraction with Miller indices
*hkl* selects the X-ray energy *energy*.
"""
lam = self.get_wavelength(energy)
d = self.get_interplanar_distance(hkl)
quot = lam/(2*d)
if quot < 1:
quot = quot.free()
angle = units.Angle(math.asin(quot), 'rad')
return(angle)
else:
return None
[docs] def get_energy(self, bragg_angle, hkl):
"""
Get the energy which is selected under a given angle *bragg_angle* for a
diffraction with Miller indices *hkl*.
"""
units.check_parameter("Bragg angle", bragg_angle)
d = self.get_interplanar_distance(hkl)
angle = bragg_angle.convert("rad").free()
wavelength = 2*d*math.sin(angle)
hc = c.planck_constant_ev*c.light_speed
energy_val = hc / wavelength.convert('m').free()
energy = units.Energy( energy_val, 'eV')
return(energy)
[docs] def get_wavelength(self, energy):
"""Get the wavelength from X-ray energy *energy*."""
units.check_parameter("Energy", energy)
eh = energy.convert('eV').free() / c.planck_constant_ev
return units.Length(c.light_speed / eh, 'm')
[docs] def get_interplanar_distance(self, hkl):
"""Get the interplanar distance for a set of Miller indices *hkl*."""
h, k, l = hkl
d_val = self.lattice_constant.free() / math.sqrt(h**2+k**2+l**2)
return units.Length(d_val, self.lattice_constant.dimension)
[docs] def get_spectral_bandwidth(self, energy, hkl, theta, polarization = 'sigma'):
"""
Compute the spectral bandwidth :math:`\\epsilon^{(s)}_H` [#f1]_ . Debye-Waller factor and anomalous
scattering are neglected, i.e. the structure factor for diamond cubic
lattice is
.. math::
F_{hkl} = f [ 1 + (-1)^{h + k} + (-1)^{k + l} + (-1)^{h + l}] * [ 1 + (-i)^{h + k + l}]
with atomic form factor :math:`f` and Miller indizes :math:`h,k,l`.
.. rubric:: References
.. [#f1] Shvyd’ko, Y. X-Ray Optics. 98, (Springer Berlin Heidelberg, 2004), p. 76.
"""
units.check_parameter("Energy", energy)
units.check_parameter("Bragg angle", theta)
theta = theta.convert('rad').free()
# polarization factor P
if polarization is 'sigma':
p = 1
elif polarization is 'pi':
p = math.cos(2*theta)
else:
raise ValueError('Polarization must be "sigma" or "pi".')
if self.material_string == "Si" or self.material_string == "Silicon":
fac_a, fac_b, fac_c = c.f0_a_si, c.f0_b_si, c.f0_c_si
elif self.material_string == "Ge" or self.material_string == "Germanium":
fac_a, fac_b, fac_c = c.f0_a_ge, c.f0_b_ge, c.f0_c_ge
else:
fmt = "Atomic form factors only tabulated for Si or Ge, not {}."
raise ValueError(fmt.format(self.material_string))
def get_f0(k):
"""
Return the atomic form factor f0 with the interpolation formula
f0[k] = c + [SUM a_i*EXP(-b_i*(k^2)) ]
i=1,5
for a given sin(theta)/lambda *k*.
"""
r = fac_c + sum([a_i*math.exp(-b_i*k**2) for a_i, b_i in zip(fac_a,fac_b)])
return r + 0j
k = math.sin(theta) / self.get_wavelength(energy).convert('angstrom').free()
f0 = get_f0(k)
# Get the structure factor (neglecting scattering corrections and
# Debye-Waller factor)
def get_F_hkl(hkl, lattice):
"""
(e.g. Shvyd’ko 2004, pp. 39-40).
Shvyd’ko, Y. X-Ray Optics. 98, (Springer Berlin Heidelberg, 2004).
"""
h,k,l = hkl
F = f0*(1 + (-1)**(h+k) + (-1)**(k+l) + (-1)**(h+l))*(1+(-1j)**(h+k+l))
return F
f_hkl = get_F_hkl(hkl, c.diamond_fcc)
d = self.get_interplanar_distance(hkl).convert('m').free()
v = self.lattice_constant.convert('m').free()**3
eps = 4*abs(p)*c.electron_radius*d**2 / (math.pi * v) * abs(f_hkl)
return eps * energy
[docs]class CylindricalCrystal(Crystal):
"""
An analyzer of cylindrical type as used in von Hamos type spectrometers.
"""
def __init__(self, width, height, *args, **kwargs):
self.width = width
self.height = height
Crystal.__init__(self, *args, **kwargs)
[docs] def get_solid_angle(self, angle):
"""
Compute the solid angle which is covered by this analyzers surface for a
given *angle*. As a general rule of thumb, an analyzer will collect
more photons if the solid angle is bigger.
"""
units.check_parameter('Bragg angle', angle)
angle = angle.convert('rad').free()
d = self.roc / math.sin(angle)
# Compute the horizontal analyzer-covered arc around the sample point
r = self.roc.convert('mm').free()
w = self.width.convert('mm').free()
h = r - 0.5*math.sqrt(4*r**2-w**2)
phi = math.atan((self.width/2 / (d.convert('mm').free()-h)).free())
hor_arc = 2* d * phi
# Compute the vertical analyzer-covered arc around the sample point
a = math.sin(angle)*self.height/2
r1 = d - math.cos(angle)*self.height/2
r2 = d + math.cos(angle)*self.height/2
vert_arc = math.atan((a/r1).free())*r1 + math.atan((a/r2).free())*r2
# Solid angle = A / r**2
A = (hor_arc * vert_arc).convert('mm').free()
sa = A / d.convert('mm').free()**2
return units.SolidAngle(sa, 'sr')