diff --git a/diffsims/__init__.py b/diffsims/__init__.py
index 46010466..ae976f42 100644
--- a/diffsims/__init__.py
+++ b/diffsims/__init__.py
@@ -17,10 +17,6 @@
# along with diffsims. If not, see .
import logging
-import os
-import warnings
-
-import numpy as np
from .generators.diffraction_generator import (
DiffractionGenerator,
@@ -31,6 +27,9 @@
from .sims.diffraction_simulation import DiffractionSimulation
+from .crystallography import * # What's imported is specified in the modules' init
+from .structure_factor import * # What's imported is specified in the modules' init
+
from . import release_info
__version__ = release_info.version
diff --git a/diffsims/crystallography/__init__.py b/diffsims/crystallography/__init__.py
new file mode 100644
index 00000000..46372805
--- /dev/null
+++ b/diffsims/crystallography/__init__.py
@@ -0,0 +1,35 @@
+# -*- coding: utf-8 -*-
+# Copyright 2017-2020 The diffsims developers
+#
+# This file is part of diffsims.
+#
+# diffsims is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# diffsims is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with diffsims. If not, see .
+
+"""Generation of reciprocal lattice points (crystal plane, reflector, g,
+hkl) for a crystal structure.
+"""
+
+from diffsims.crystallography.reciprocal_lattice_point import (
+ ReciprocalLatticePoint,
+ get_equivalent_hkl,
+ get_highest_hkl,
+ get_hkl,
+)
+
+__all__ = [
+ "ReciprocalLatticePoint",
+ "get_equivalent_hkl",
+ "get_highest_hkl",
+ "get_hkl",
+]
diff --git a/diffsims/crystallography/reciprocal_lattice_point.py b/diffsims/crystallography/reciprocal_lattice_point.py
new file mode 100644
index 00000000..7afd1cd3
--- /dev/null
+++ b/diffsims/crystallography/reciprocal_lattice_point.py
@@ -0,0 +1,458 @@
+# -*- coding: utf-8 -*-
+# Copyright 2017-2020 The diffsims developers
+#
+# This file is part of diffsims.
+#
+# diffsims is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# diffsims is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with diffsims. If not, see .
+
+from collections import defaultdict
+from itertools import product
+
+import numpy as np
+from orix.vector import Vector3d
+
+from diffsims.structure_factor.structure_factor import (
+ get_kinematical_structure_factor,
+ get_doyleturner_structure_factor,
+ get_refraction_corrected_wavelength,
+)
+
+
+_FLOAT_EPS = np.finfo(np.float).eps # Used to round values below 1e-16 to zero
+
+
+class ReciprocalLatticePoint:
+ """Reciprocal lattice point (or crystal plane, reflector, g, etc.)
+ with Miller indices, length of the reciprocal lattice vectors and
+ other relevant structure_factor parameters.
+ """
+
+ def __init__(self, phase, hkl):
+ """A container for Miller indices, structure factors and related
+ parameters for crystal planes (reciprocal lattice points,
+ reflectors, g, etc.).
+
+ Parameters
+ ----------
+ phase : orix.crystal_map.phase_list.Phase
+ A phase container with a crystal structure and a space and
+ point group describing the allowed symmetry operations.
+ hkl : orix.vector.Vector3d, np.ndarray, list, or tuple
+ Miller indices.
+ """
+ self._hkl = Vector3d(hkl)
+ self.phase = phase
+ self._structure_factor = [None] * self.size
+ self._theta = [None] * self.size
+
+ def __repr__(self):
+ return (
+ f"{self.__class__.__name__} {self.hkl.shape}\n"
+ f"Phase: {self.phase.name} ({self.phase.point_group.name})\n"
+ f"{np.array_str(self.hkl.data, precision=4, suppress_small=True)}"
+ )
+
+ def __getitem__(self, key):
+ new_rlp = self.__class__(self.phase, self.hkl[key])
+ if self.structure_factor[0] is None:
+ new_rlp._structure_factor = [None] * new_rlp.size
+ else:
+ new_rlp._structure_factor = self.structure_factor[key]
+ if self.theta[0] is None:
+ new_rlp._theta = [None] * new_rlp.size
+ else:
+ new_rlp._theta = self.theta[key]
+ return new_rlp
+
+ @property
+ def hkl(self):
+ """Return :class:`~orix.vector.Vector3d` of Miller indices."""
+ return Vector3d(self._hkl.data.astype(int))
+
+ @property
+ def _hkldata(self):
+ """Return :class:`np.ndarray` without 1-dimensions."""
+ return np.squeeze(self.hkl.data)
+
+ @property
+ def h(self):
+ """Return :class:`np.ndarray` of Miller index h."""
+ return self._hkldata[..., 0]
+
+ @property
+ def k(self):
+ """Return :class:`np.ndarray` of Miller index k."""
+ return self._hkldata[..., 1]
+
+ @property
+ def l(self):
+ """Return :class:`np.ndarray` of Miller index l."""
+ return self._hkldata[..., 2]
+
+ @property
+ def size(self):
+ """Return `int`."""
+ return self.hkl.size
+
+ @property
+ def shape(self):
+ """Return `tuple`."""
+ return self._hkldata.shape
+
+ @property
+ def multiplicity(self):
+ """Return either `int` or :class:`np.ndarray` of `int`."""
+ return self.symmetrise(antipodal=True, return_multiplicity=True)[1]
+
+ @property
+ def gspacing(self):
+ """Return :class:`np.ndarray` of reciprocal lattice point
+ spacings.
+ """
+ return self.phase.structure.lattice.rnorm(self._hkldata)
+
+ @property
+ def dspacing(self):
+ """Return :class:`np.ndarray` of direct lattice interplanar
+ spacings.
+ """
+ return 1 / self.gspacing
+
+ @property
+ def scattering_parameter(self):
+ """Return :class:`np.ndarray` of scattering parameters s."""
+ return 0.5 * self.gspacing
+
+ @property
+ def structure_factor(self):
+ """Return :class:`np.ndarray` of structure factors F or None."""
+ return self._structure_factor
+
+ @property
+ def allowed(self):
+ """Return whether planes diffract according to structure_factor
+ selection rules assuming kinematical scattering theory.
+ """
+ # Translational symmetry
+ centering = self.phase.space_group.short_name[0]
+
+ if centering == "P": # Primitive
+ if self.phase.space_group.crystal_system == "HEXAGONAL":
+ # TODO: See rules in e.g.
+ # https://mcl1.ncifcrf.gov/dauter_pubs/284.pdf, Table 4
+ # http://xrayweb.chem.ou.edu/notes/symmetry.html, Systematic Absences
+ raise NotImplementedError
+ else: # Any hkl
+ return np.ones(self.size, dtype=bool)
+ elif centering == "F": # Face-centred, hkl all odd/even
+ selection = np.sum(np.mod(self._hkldata, 2), axis=1)
+ return np.array([i not in [1, 2] for i in selection], dtype=bool)
+ elif centering == "I": # Body-centred, h + k + l = 2n (even)
+ return np.mod(np.sum(self._hkldata, axis=1), 2) == 0
+ elif centering == "A": # Centred on A faces only
+ return np.mod(self.k + self.l, 2) == 0
+ elif centering == "B": # Centred on B faces only
+ return np.mod(self.h + self.l, 2) == 0
+ elif centering == "C": # Centred on C faces only
+ return np.mod(self.h + self.k, 2) == 0
+ elif centering in ["R", "H"]: # Rhombohedral
+ return np.mod(-self.h + self.k + self.l, 3) == 0
+
+ @property
+ def theta(self):
+ """Return :class:`np.ndarray` of twice the Bragg angle."""
+ return self._theta
+
+ @classmethod
+ def from_min_dspacing(cls, phase, min_dspacing=0.5):
+ """Create a CrystalPlane object populated by unique Miller indices
+ with a direct space interplanar spacing greater than a lower
+ threshold.
+
+ Parameters
+ ----------
+ phase : orix.crystal_map.phase_list.Phase
+ A phase container with a crystal structure and a space and
+ point group describing the allowed symmetry operations.
+ min_dspacing : float, optional
+ Smallest interplanar spacing to consider. Default is 0.5 Å.
+ """
+ highest_hkl = get_highest_hkl(
+ lattice=phase.structure.lattice, min_dspacing=min_dspacing
+ )
+ hkl = get_hkl(highest_hkl=highest_hkl)
+ return cls(phase=phase, hkl=hkl).unique()
+
+ @classmethod
+ def from_highest_hkl(cls, phase, highest_hkl):
+ """Create a CrystalPlane object populated by unique Miller indices
+ below, but including, a set of higher indices.
+
+ Parameters
+ ----------
+ phase : orix.crystal_map.phase_list.Phase
+ A phase container with a crystal structure and a space and
+ point group describing the allowed symmetry operations.
+ highest_hkl : np.ndarray, list, or tuple of int
+ Highest Miller indices to consider (including).
+ """
+ hkl = get_hkl(highest_hkl=highest_hkl)
+ return cls(phase=phase, hkl=hkl).unique()
+
+ def unique(self, use_symmetry=True):
+ """Return planes with unique Miller indices.
+
+ Parameters
+ ----------
+ use_symmetry : bool, optional
+ Whether to use symmetry to remove the planes with indices
+ symmetrically equivalent to another set of indices.
+
+ Returns
+ -------
+ ReciprocalLatticePoint
+ """
+ if use_symmetry:
+ all_hkl = self._hkldata
+ # Remove [0, 0, 0] points
+ all_hkl = all_hkl[~np.all(np.isclose(all_hkl, 0), axis=1)]
+
+ families = defaultdict(list)
+ for this_hkl in all_hkl.tolist():
+ for that_hkl in families.keys():
+ if _is_equivalent(this_hkl, that_hkl):
+ families[tuple(that_hkl)].append(this_hkl)
+ break
+ else:
+ families[tuple(this_hkl)].append(this_hkl)
+
+ n_families = len(families)
+ unique_hkl = np.zeros((n_families, 3))
+ for i, all_hkl_in_family in enumerate(families.values()):
+ unique_hkl[i] = sorted(all_hkl_in_family)[-1]
+ else:
+ unique_hkl = self.hkl.unique()
+ # TODO: Enable inheriting classes pass on their properties in this new object
+ return self.__class__(phase=self.phase, hkl=unique_hkl)
+
+ def symmetrise(
+ self,
+ antipodal=True,
+ unique=True,
+ return_multiplicity=False,
+ ):
+ """Return planes with symmetrically equivalent Miller indices.
+
+ Parameters
+ ----------
+ antipodal : bool, optional
+ Whether to include antipodal symmetry operations. Default is
+ True.
+ unique : bool, optional
+ Whether to return only distinct indices. Default is True.
+ If True, zero-entries, which are assumed to be degenerate, are
+ removed.
+ return_multiplicity : bool, optional
+ Whether to return the multiplicity of indices. This option is
+ only available if `unique` is True. Default is False.
+
+ Returns
+ -------
+ ReciprocalLatticePoint
+ Planes with Miller indices symmetrically equivalent to the
+ original planes.
+ multiplicity : np.ndarray
+ Multiplicity of the original Miller indices. Only returned if
+ `return_multiplicity` is True.
+
+ Notes
+ -----
+ Should be the same as EMsoft's CalcFamily in their symmetry.f90
+ module, although not entirely sure. Use with care.
+ """
+ # Get symmetry operations
+ pg = self.phase.point_group
+ operations = pg if antipodal else pg[~pg.improper]
+
+ out = get_equivalent_hkl(
+ hkl=self.hkl,
+ operations=operations,
+ unique=unique,
+ return_multiplicity=return_multiplicity,
+ )
+
+ # TODO: Enable inheriting classes pass on their properties in this new object
+ # Format output and return
+ if unique and return_multiplicity:
+ multiplicity = out[1]
+ if multiplicity.size == 1:
+ multiplicity = multiplicity[0]
+ return self.__class__(phase=self.phase, hkl=out[0]), multiplicity
+ else:
+ return self.__class__(phase=self.phase, hkl=out)
+
+ def calculate_structure_factor(self, method=None, voltage=None):
+ """Populate `self.structure_factor` with the structure factor F
+ for each plane.
+
+ Parameters
+ ----------
+ method : str, optional
+ Either "kinematical" for kinematical X-ray structure factors
+ or "doyleturner" for structure factors using Doyle-Turner
+ atomic scattering factors. If None (default), kinematical
+ structure factors are calculated.
+ voltage : float, optional
+ Beam energy in V used when `method=doyleturner`.
+ """
+ if method is None:
+ method = "kinematical"
+ methods = ["kinematical", "doyleturner"]
+ if method not in methods:
+ raise ValueError(f"method={method} must be among {methods}")
+ elif method == "doyleturner" and voltage is None:
+ raise ValueError(
+ "'voltage' parameter must be set when method='doyleturner'"
+ )
+
+ # TODO: Find a better way to call different methods in the loop
+ structure_factors = np.zeros(self.size)
+ for i, (hkl, s) in enumerate(
+ zip(np.atleast_2d(self._hkldata), np.atleast_1d(self.scattering_parameter))
+ ):
+ if method == "kinematical":
+ structure_factors[i] = get_kinematical_structure_factor(
+ phase=self.phase,
+ hkl=hkl,
+ scattering_parameter=s,
+ )
+ else:
+ structure_factors[i] = get_doyleturner_structure_factor(
+ phase=self.phase,
+ hkl=hkl,
+ scattering_parameter=s,
+ voltage=voltage,
+ )
+ self._structure_factor = np.where(
+ structure_factors < _FLOAT_EPS, 0, structure_factors
+ )
+
+ def calculate_theta(self, voltage):
+ """Populate `self.theta` with the Bragg angle :math:`theta_B` for
+ each plane.
+
+ Parameters
+ ----------
+ voltage : float
+ Beam energy in V.
+ """
+ wavelength = get_refraction_corrected_wavelength(self.phase, voltage)
+ self._theta = np.arcsin(0.5 * wavelength * self.gspacing)
+
+
+def get_highest_hkl(lattice, min_dspacing=0.5):
+ """Return the highest Miller indices hkl of the plane with a direct
+ space interplanar spacing greater than but closest to a lower
+ threshold.
+
+ Parameters
+ ----------
+ lattice : diffpy.structure.Lattice
+ Crystal lattice.
+ min_dspacing : float, optional
+ Smallest interplanar spacing to consider. Default is 0.5 Å.
+
+ Returns
+ -------
+ highest_hkl : np.ndarray
+ Highest Miller indices.
+ """
+ highest_hkl = np.ones(3, dtype=int)
+ for i in range(3):
+ hkl = np.zeros(3)
+ d = min_dspacing + 1
+ while d > min_dspacing:
+ hkl[i] += 1
+ d = 1 / lattice.rnorm(hkl)
+ highest_hkl[i] = hkl[i]
+ return highest_hkl
+
+
+def get_hkl(highest_hkl):
+ """Return a list of planes from a set of highest Miller indices.
+
+ Parameters
+ ----------
+ highest_hkl : orix.vector.Vector3d, np.ndarray, list, or tuple of int
+ Highest Miller indices to consider.
+
+ Returns
+ -------
+ hkl : np.ndarray
+ An array of Miller indices.
+ """
+ index_ranges = [np.arange(-i, i + 1) for i in highest_hkl]
+ return np.asarray(list(product(*index_ranges)))
+
+
+def get_equivalent_hkl(hkl, operations, unique=False, return_multiplicity=False):
+ """Return symmetrically equivalent Miller indices.
+
+ Parameters
+ ----------
+ hkl : orix.vector.Vector3d, np.ndarray, list or tuple of int
+ Miller indices.
+ operations : orix.quaternion.symmetry.Symmetry
+ Point group describing allowed symmetry operations.
+ unique : bool, optional
+ Whether to return only unique Miller indices. Default is False.
+ return_multiplicity : bool, optional
+ Whether to return the multiplicity of the input indices. Default
+ is False.
+
+ Returns
+ -------
+ new_hkl : orix.vector.Vector3d
+ The symmetrically equivalent Miller indices.
+ multiplicity : np.ndarray
+ Number of symmetrically equivalent indices. Only returned if
+ `return_multiplicity` is True.
+ """
+ new_hkl = operations.outer(Vector3d(hkl))
+ new_hkl = new_hkl.flatten().reshape(*new_hkl.shape[::-1])
+
+ multiplicity = None
+ if unique:
+ n_families = new_hkl.shape[0]
+ multiplicity = np.zeros(n_families, dtype=int)
+ temp_hkl = new_hkl[0].unique().data
+ multiplicity[0] = temp_hkl.shape[0]
+ if n_families > 1:
+ for i, hkl in enumerate(new_hkl[1:]):
+ temp_hkl2 = hkl.unique()
+ multiplicity[i + 1] = temp_hkl2.size
+ temp_hkl = np.append(temp_hkl, temp_hkl2.data, axis=0)
+ new_hkl = Vector3d(temp_hkl[: multiplicity.sum()])
+
+ # Remove 1-dimensions
+ new_hkl = new_hkl.squeeze()
+
+ if unique and return_multiplicity:
+ return new_hkl, multiplicity
+ else:
+ return new_hkl
+
+
+def _is_equivalent(this_hkl: list, that_hkl: list) -> bool:
+ return sorted(np.abs(this_hkl)) == sorted(np.abs(that_hkl))
diff --git a/diffsims/generators/diffraction_generator.py b/diffsims/generators/diffraction_generator.py
index e1b62a73..0915e63f 100644
--- a/diffsims/generators/diffraction_generator.py
+++ b/diffsims/generators/diffraction_generator.py
@@ -36,9 +36,10 @@
get_vectorized_list_for_atomic_scattering_factors,
is_lattice_hexagonal,
get_intensities_params,
- get_scattering_params_dict
+ get_scattering_params_dict,
)
from diffsims.utils.fourier_transform import from_recip
+from diffsims.utils.shape_factor_models import linear, binary, sinc
class DiffractionGenerator(object):
@@ -72,7 +73,7 @@ def __init__(
accelerating_voltage,
max_excitation_error=None,
debye_waller_factors={},
- scattering_params="lobato"
+ scattering_params="lobato",
):
if max_excitation_error is not None:
print(
@@ -80,7 +81,7 @@ def __init__(
)
self.wavelength = get_electron_wavelength(accelerating_voltage)
self.debye_waller_factors = debye_waller_factors
- if scattering_params in ["lobato","xtables"]:
+ if scattering_params in ["lobato", "xtables"]:
self.scattering_params = scattering_params
else:
raise NotImplementedError(
@@ -94,7 +95,7 @@ def calculate_ed_data(
structure,
reciprocal_radius,
rotation=(0, 0, 0),
- shape_factor_model="linear",
+ shape_factor_model=None,
max_excitation_error=1e-2,
with_direct_beam=True,
**kwargs
@@ -113,9 +114,9 @@ def calculate_ed_data(
rotation : tuple
Euler angles, in degrees, in the rzxz convention. Default is (0,0,0)
which aligns 'z' with the electron beam
- shape_factor_model : function or str
+ shape_factor_model : function or None
a function that takes excitation_error and max_excitation_error (and potentially **kwargs) and returns an intensity
- scaling factor. The code provides "linear" and "binary" options accessed with by parsing the associated strings
+ scaling factor. If None defaults to shape_factor_models.linear
max_excitation_error : float
the exctinction distance for reflections, in reciprocal Angstroms
with_direct_beam : bool
@@ -162,14 +163,12 @@ def calculate_ed_data(
excitation_error = excitation_error[intersection]
g_hkls = spot_distances[intersection]
- if shape_factor_model == "linear":
- shape_factor = 1 - (excitation_error / max_excitation_error)
- elif shape_factor_model == "binary":
- shape_factor = 1
- else:
+ if shape_factor_model is not None:
shape_factor = shape_factor_model(
excitation_error, max_excitation_error, **kwargs
)
+ else:
+ shape_factor = linear(excitation_error, max_excitation_error)
# Calculate diffracted intensities based on a kinematical model.
intensities = get_kinematical_intensities(
diff --git a/diffsims/generators/rotation_list_generators.py b/diffsims/generators/rotation_list_generators.py
index 071774be..dd27a073 100644
--- a/diffsims/generators/rotation_list_generators.py
+++ b/diffsims/generators/rotation_list_generators.py
@@ -68,7 +68,9 @@ def get_list_from_orix(grid, rounding=2):
rotation_list = z.data.tolist()
i = 0
while i < len(rotation_list):
- rotation_list[i] = tuple(np.round(np.rad2deg(rotation_list[i]), decimals=rounding))
+ rotation_list[i] = tuple(
+ np.round(np.rad2deg(rotation_list[i]), decimals=rounding)
+ )
i += 1
return rotation_list
@@ -115,9 +117,9 @@ def get_local_grid(resolution=2, center=None, grid_width=10):
-------
rotation_list : list of tuples
"""
- if isinstance(center,tuple):
+ if isinstance(center, tuple):
z = np.deg2rad(np.asarray(center))
- center = Rotation.from_euler(z,convention="bunge",direction="crystal2lab")
+ center = Rotation.from_euler(z, convention="bunge", direction="crystal2lab")
orix_grid = get_sample_local(
resolution=resolution, center=center, grid_width=grid_width
@@ -152,14 +154,16 @@ def get_grid_around_beam_direction(beam_rotation, resolution, angular_range=(0,
>>> grid = get_grid_around_beam_direction(beam_rotation,1)
"""
z = np.deg2rad(np.asarray(beam_rotation))
- beam_rotation = Rotation.from_euler(z,convention="bunge",direction="crystal2lab")
+ beam_rotation = Rotation.from_euler(z, convention="bunge", direction="crystal2lab")
- angles = np.deg2rad(np.arange(start=angular_range[0],stop=angular_range[1],step=resolution))
- axes = np.repeat([[0,0,1]],angles.shape[0],axis=0)
- in_plane_rotation = Rotation.from_neo_euler(AxAngle.from_axes_angles(axes,angles))
+ angles = np.deg2rad(
+ np.arange(start=angular_range[0], stop=angular_range[1], step=resolution)
+ )
+ axes = np.repeat([[0, 0, 1]], angles.shape[0], axis=0)
+ in_plane_rotation = Rotation.from_neo_euler(AxAngle.from_axes_angles(axes, angles))
orix_grid = beam_rotation * in_plane_rotation
- rotation_list = get_list_from_orix(orix_grid,rounding=2)
+ rotation_list = get_list_from_orix(orix_grid, rounding=2)
return rotation_list
diff --git a/diffsims/release_info.py b/diffsims/release_info.py
index 5c35c942..bc9b1382 100644
--- a/diffsims/release_info.py
+++ b/diffsims/release_info.py
@@ -2,7 +2,13 @@
version = "0.2.3"
author = "Duncan Johnstone, Phillip Crout"
copyright = "Copyright 2017-2020, The pyXem Developers"
-credits = ["Duncan Johnstone", "Phillip Crout", "Ben Martineau", "Simon Hogas"]
+credits = [
+ "Duncan Johnstone",
+ "Phillip Crout",
+ "Ben Martineau",
+ "Simon Hogas",
+ "Håkon Wiik Ånes",
+]
license = "GPLv3"
maintainer = "Duncan Johnstone, Phillip Crout"
email = "pyxem.team@gmail.com"
diff --git a/diffsims/structure_factor/__init__.py b/diffsims/structure_factor/__init__.py
new file mode 100644
index 00000000..52befcb9
--- /dev/null
+++ b/diffsims/structure_factor/__init__.py
@@ -0,0 +1,45 @@
+# -*- coding: utf-8 -*-
+# Copyright 2017-2020 The diffsims developers
+#
+# This file is part of diffsims.
+#
+# diffsims is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# diffsims is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with diffsims. If not, see .
+
+"""Calculation of scattering factors and structure factors."""
+
+from diffsims.structure_factor.atomic_scattering_factor import (
+ get_doyleturner_atomic_scattering_factor,
+ get_kinematical_atomic_scattering_factor,
+)
+from diffsims.structure_factor.atomic_scattering_parameters import (
+ get_atomic_scattering_parameters,
+ get_element_id_from_string,
+)
+from diffsims.structure_factor.structure_factor import (
+ find_asymmetric_positions,
+ get_doyleturner_structure_factor,
+ get_kinematical_structure_factor,
+ get_refraction_corrected_wavelength,
+)
+
+__all__ = [
+ "get_doyleturner_atomic_scattering_factor",
+ "get_kinematical_atomic_scattering_factor",
+ "get_atomic_scattering_parameters",
+ "get_element_id_from_string",
+ "find_asymmetric_positions",
+ "get_doyleturner_structure_factor",
+ "get_kinematical_structure_factor",
+ "get_refraction_corrected_wavelength",
+]
diff --git a/diffsims/structure_factor/atomic_scattering_factor.py b/diffsims/structure_factor/atomic_scattering_factor.py
new file mode 100644
index 00000000..9b73cadf
--- /dev/null
+++ b/diffsims/structure_factor/atomic_scattering_factor.py
@@ -0,0 +1,105 @@
+# -*- coding: utf-8 -*-
+# Copyright 2017-2020 The diffsims developers
+#
+# This file is part of diffsims.
+#
+# diffsims is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# diffsims is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with diffsims. If not, see .
+
+import numpy as np
+
+from diffsims.structure_factor.atomic_scattering_parameters import (
+ get_atomic_scattering_parameters,
+ get_element_id_from_string,
+)
+
+
+def get_kinematical_atomic_scattering_factor(atom, scattering_parameter):
+ """Return the kinematical (X-ray) atomic scattering factor f for a
+ certain atom and scattering parameter.
+
+ Assumes structure's Debye-Waller factors are expressed in Ångströms.
+
+ This function is adapted from EMsoft.
+
+ Parameters
+ ----------
+ atom : diffpy.structure.Atom
+ Atom with element type, Debye-Waller factor and occupancy number.
+ scattering_parameter : float
+ The scattering parameter s for these Miller indices describing
+ the crystal plane in which the atom lies.
+
+ Returns
+ -------
+ f : float
+ Scattering factor for this atom on this plane.
+ """
+ # Get the atomic scattering parameters
+ element_id = get_element_id_from_string(atom.element)
+ a, b = get_atomic_scattering_parameters(element_id)
+
+ # Get the scattering parameter squared
+ s2 = scattering_parameter ** 2
+
+ # Get the atomic scattering factor
+ f = element_id - (41.78214 * s2 * np.sum(a * np.exp(-b * s2)))
+
+ # Correct for occupancy and the Debye-Waller factor
+ dw_factor = np.exp(-atom.Bisoequiv * s2)
+ f *= atom.occupancy * dw_factor
+
+ return f
+
+
+def get_doyleturner_atomic_scattering_factor(
+ atom, scattering_parameter, unit_cell_volume
+):
+ """Return the atomic scattering factor f for a certain atom and
+ scattering parameter using Doyle-Turner atomic scattering parameters
+ [Doyle1968]_.
+
+ Assumes structure's Debye-Waller factors are expressed in Ångströms.
+
+ This function is adapted from EMsoft.
+
+ Parameters
+ ----------
+ atom : diffpy.structure.Atom
+ Atom with element type, Debye-Waller factor and occupancy number.
+ scattering_parameter : float
+ The scattering parameter s for these Miller indices describing
+ the crystal plane in which the atom lies.
+ unit_cell_volume : float
+ Volume of the unit cell.
+
+ Returns
+ -------
+ f : float
+ Scattering factor for this atom on this plane.
+ """
+ # Get the atomic scattering parameters
+ element_id = get_element_id_from_string(atom.element)
+ a, b = get_atomic_scattering_parameters(element_id)
+
+ # Get the scattering parameter squared
+ s2 = scattering_parameter ** 2
+
+ # Get the atomic scattering factor
+ f = (47.87801 / unit_cell_volume) * np.sum(a * np.exp(-b * s2))
+
+ # Correct for occupancy and the Debye-Waller factor
+ dw_factor = np.exp(-atom.Bisoequiv * s2)
+ f *= atom.occupancy * dw_factor
+
+ return f
diff --git a/diffsims/structure_factor/atomic_scattering_parameters.py b/diffsims/structure_factor/atomic_scattering_parameters.py
new file mode 100644
index 00000000..019b8a8c
--- /dev/null
+++ b/diffsims/structure_factor/atomic_scattering_parameters.py
@@ -0,0 +1,208 @@
+# -*- coding: utf-8 -*-
+# Copyright 2017-2020 The diffsims developers
+#
+# This file is part of diffsims.
+#
+# diffsims is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# diffsims is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with diffsims. If not, see .
+
+import numpy as np
+
+# List of elements Z = 1-118
+# fmt: off
+ELEMENTS = [
+ "h", "he", "li", "be", "b", "c", "n", "o", "f", "ne", "na", "mg", "al",
+ "si", "p", "s", "cl", "ar", "k", "ca", "sc", "ti", "v", "cr", "mn",
+ "fe", "co", "ni", "cu", "zn", "ga", "ge", "as", "se", "br", "kr", "rb",
+ "sr", "y", "zr", "nb", "mo", "tc", "ru", "rh", "pd", "ag", "cd", "in",
+ "sn", "sb", "te", "i", "xe", "cs", "ba", "la", "ce", "pr", "nd", "pm",
+ "sm", "eu", "gd", "tb", "dy", "ho", "er", "tm", "yb", "lu", "hf", "ta",
+ "w", "re", "os", "ir", "pt", "au", "hg", "tl", "pb", "bi", "po", "at",
+ "rn", "fr", "ra", "ac", "th", "pa", "u", "np", "pu", "am", "cm", "bk",
+ "cf", "es", "fm", "md", "no", "lr", "rf", "db", "sg", "bh", "hs", "mt",
+ "ds", "rg", "cn", "nh", "fl", "mc", "lv", "ts", "og"
+]
+# fmt: on
+N_ELEMENTS = len(ELEMENTS)
+
+# fmt: off
+ATOMIC_SCATTERING_PARAMETERS_DOYLETURNER = np.array([
+ # a1 a2 a3 a4 b1 b2 b3 b4
+ 0.202, 0.244, 0.082, 0.000, 30.868, 8.544, 1.273, 0.000, # H
+ 0.091, 0.181, 0.110, 0.036, 18.183, 6.212, 1.803, 0.284, # He
+ 1.611, 1.246, 0.326, 0.099, 107.638, 30.480, 4.533, 0.495, # etc.
+ 1.250, 1.334, 0.360, 0.106, 60.804, 18.591, 3.653, 0.416,
+ 0.945, 1.312, 0.419, 0.116, 46.444, 14.178, 3.223, 0.377,
+ 0.731, 1.195, 0.456, 0.125, 36.995, 11.297, 2.814, 0.346,
+ 0.572, 1.043, 0.465, 0.131, 28.847, 9.054, 2.421, 0.317,
+ 0.455, 0.917, 0.472, 0.138, 23.780, 7.622, 2.144, 0.296,
+ 0.387, 0.811, 0.475, 0.146, 20.239, 6.609, 1.931, 0.279,
+ 0.303, 0.720, 0.475, 0.153, 17.640, 5.860, 1.762, 0.266,
+ 2.241, 1.333, 0.907, 0.286, 108.004, 24.505, 3.391, 0.435,
+ 2.268, 1.803, 0.839, 0.289, 73.670, 20.175, 3.013, 0.405,
+ 2.276, 2.428, 0.858, 0.317, 72.322, 19.773, 3.080, 0.408,
+ 2.129, 2.533, 0.835, 0.322, 57.775, 16.476, 2.880, 0.386,
+ 1.888, 2.469, 0.805, 0.320, 44.876, 13.538, 2.642, 0.361,
+ 1.659, 2.386, 0.790, 0.321, 36.650, 11.488, 2.469, 0.340,
+ 1.452, 2.292, 0.787, 0.322, 30.935, 9.980, 2.234, 0.323,
+ 1.274, 2.190, 0.793, 0.326, 26.682, 8.813, 2.219, 0.307,
+ 3.951, 2.545, 1.980, 0.482, 137.075, 22.402, 4.532, 0.434,
+ 4.470, 2.971, 1.970, 0.482, 99.523, 22.696, 4.195, 0.417,
+ 3.966, 2.917, 1.925, 0.480, 88.960, 20.606, 3.856, 0.399,
+ 3.565, 2.818, 1.893, 0.483, 81.982, 19.049, 3.590, 0.386,
+ 3.245, 2.698, 1.860, 0.486, 76.379, 17.726, 3.363, 0.374,
+ 2.307, 2.334, 1.823, 0.490, 78.405, 15.785, 3.157, 0.364,
+ 2.747, 2.456, 1.792, 0.498, 67.786, 15.674, 3.000, 0.357,
+ 2.544, 2.343, 1.759, 0.506, 64.424, 14.880, 2.854, 0.350,
+ 2.367, 2.236, 1.724, 0.515, 61.431, 14.180, 2.725, 0.344,
+ 2.210, 2.134, 1.689, 0.524, 58.727, 13.553, 2.609, 0.339,
+ 1.579, 1.820, 1.658, 0.532, 62.940, 12.453, 2.504, 0.333,
+ 1.942, 1.950, 1.619, 0.543, 54.162, 12.518, 2.416, 0.330,
+ 2.321, 2.486, 1.688, 0.599, 65.602, 15.458, 2.581, 0.351,
+ 2.447, 2.702, 1.616, 0.601, 55.893, 14.393, 2.446, 0.342,
+ 2.399, 2.790, 1.529, 0.594, 45.718, 12.817, 2.280, 0.328,
+ 2.298, 2.854, 1.456, 0.590, 38.830, 11.536, 2.146, 0.316,
+ 2.166, 2.904, 1.395, 0.589, 33.899, 10.497, 2.041, 0.307,
+ 2.034, 2.927, 1.342, 0.589, 29.999, 9.598, 1.952, 0.299,
+ 4.776, 3.859, 2.234, 0.868, 140.782, 18.991, 3.701, 0.419,
+ 5.848, 4.003, 2.342, 0.880, 104.972, 19.367, 3.737, 0.414,
+ 4.129, 3.012, 1.179, 0.000, 27.548, 5.088, 0.591, 0.0,
+ 4.105, 3.144, 1.229, 0.000, 28.492, 5.277, 0.601, 0.0,
+ 4.237, 3.105, 1.234, 0.000, 27.415, 5.074, 0.593, 0.0,
+ 3.120, 3.906, 2.361, 0.850, 72.464, 14.642, 3.237, 0.366,
+ 4.318, 3.270, 1.287, 0.000, 28.246, 5.148, 0.590, 0.0,
+ 4.358, 3.298, 1.323, 0.000, 27.881, 5.179, 0.594, 0.0,
+ 4.431, 3.343, 1.345, 0.000, 27.911, 5.153, 0.592, 0.0,
+ 4.436, 3.454, 1.383, 0.000, 28.670, 5.269, 0.595, 0.0,
+ 2.036, 3.272, 2.511, 0.837, 61.497, 11.824, 2.846, 0.327,
+ 2.574, 3.259, 2.547, 0.838, 55.675, 11.838, 2.784, 0.322,
+ 3.153, 3.557, 2.818, 0.884, 66.649, 14.449, 2.976, 0.335,
+ 3.450, 3.735, 2.118, 0.877, 59.104, 14.179, 2.855, 0.327,
+ 3.564, 3.844, 2.687, 0.864, 50.487, 13.316, 2.691, 0.316,
+ 4.785, 3.688, 1.500, 0.000, 27.999, 5.083, 0.581, 0.0,
+ 3.473, 4.060, 2.522, 0.840, 39.441, 11.816, 2.415, 0.298,
+ 3.366, 4.147, 2.443, 0.829, 35.509, 11.117, 2.294, 0.289,
+ 6.062, 5.986, 3.303, 1.096, 155.837, 19.695, 3.335, 0.379,
+ 7.821, 6.004, 3.280, 1.103, 117.657, 18.778, 3.263, 0.376,
+ 4.940, 3.968, 1.663, 0.000, 28.716, 5.245, 0.594, 0.0,
+ 5.007, 3.980, 1.678, 0.000, 28.283, 5.183, 0.589, 0.0,
+ 5.085, 4.043, 1.684, 0.000, 28.588, 5.143, 0.581, 0.0,
+ 5.151, 4.075, 1.683, 0.000, 28.304, 5.073, 0.571, 0.0,
+ 5.201, 4.094, 1.719, 0.000, 28.079, 5.081, 0.576, 0.0,
+ 5.255, 4.113, 1.743, 0.000, 28.016, 5.037, 0.577, 0.0,
+ 6.267, 4.844, 3.202, 1.200, 100.298, 16.066, 2.980, 0.367,
+ 5.225, 4.314, 1.827, 0.000, 29.158, 5.259, 0.586, 0.0,
+ 5.272, 4.347, 1.844, 0.000, 29.046, 5.226, 0.585, 0.0,
+ 5.332, 4.370, 1.863, 0.000, 28.888, 5.198, 0.581, 0.0,
+ 5.376, 4.403, 1.884, 0.000, 28.773, 5.174, 0.582, 0.0,
+ 5.436, 4.437, 1.891, 0.000, 28.655, 5.117, 0.577, 0.0,
+ 5.441, 4.510, 1.956, 0.000, 29.149, 5.264, 0.590, 0.0,
+ 5.529, 4.533, 1.945, 0.000, 28.927, 5.144, 0.578, 0.0,
+ 5.553, 4.580, 1.969, 0.000, 28.907, 5.160, 0.577, 0.0,
+ 5.588, 4.619, 1.997, 0.000, 29.001, 5.164, 0.579, 0.0,
+ 5.659, 4.630, 2.014, 0.000, 28.807, 5.114, 0.578, 0.0,
+ 5.709, 4.677, 2.019, 0.000, 28.782, 5.084, 0.572, 0.0,
+ 5.695, 4.740, 2.064, 0.000, 28.968, 5.156, 0.575, 0.0,
+ 5.750, 4.773, 2.079, 0.000, 28.933, 5.139, 0.573, 0.0,
+ 5.754, 4.851, 2.096, 0.000, 29.159, 5.152, 0.570, 0.0,
+ 5.803, 4.870, 2.127, 0.000, 29.016, 5.150, 0.572, 0.0,
+ 2.388, 4.226, 2.689, 1.255, 42.866, 9.743, 2.264, 0.307,
+ 2.682, 4.241, 2.755, 1.270, 42.822, 9.856, 2.295, 0.307,
+ 5.932, 4.972, 2.195, 0.000, 29.086, 5.126, 0.572, 0.0,
+ 3.510, 4.552, 3.154, 1.359, 52.914, 11.884, 2.571, 0.321,
+ 3.841, 4.679, 3.192, 1.363, 50.261, 11.999, 2.560, 0.318,
+ 6.070, 4.997, 2.232, 0.000, 28.075, 4.999, 0.563, 0.0,
+ 6.133, 5.031, 2.239, 0.000, 28.047, 4.957, 0.558, 0.0,
+ 4.078, 4.978, 3.096, 1.326, 38.406, 11.020, 2.355, 0.299,
+ 6.201, 5.121, 2.275, 0.000, 28.200, 4.954, 0.556, 0.0,
+ 6.215, 5.170, 2.316, 0.000, 28.382, 5.002, 0.562, 0.0,
+ 6.278, 5.195, 2.321, 0.000, 28.323, 4.949, 0.557, 0.0,
+ 6.264, 5.263, 2.367, 0.000, 28.651, 5.030, 0.563, 0.0,
+ 6.306, 5.303, 2.386, 0.000, 28.688, 5.026, 0.561, 0.0,
+ 6.767, 6.729, 4.014, 1.561, 85.951, 15.642, 2.936, 0.335,
+ 6.323, 5.414, 2.453, 0.000, 29.142, 5.096, 0.568, 0.0,
+ 6.415, 5.419, 2.449, 0.000, 28.836, 5.022, 0.561, 0.0,
+ 6.378, 5.495, 2.495, 0.000, 29.156, 5.102, 0.565, 0.0,
+ 6.460, 5.469, 2.471, 0.000, 28.396, 4.970, 0.554, 0.0,
+ 6.502, 5.478, 2.510, 0.000, 28.375, 4.975, 0.561, 0.0,
+ 6.548, 5.526, 2.520, 0.000, 28.461, 4.965, 0.557, 0.0,
+]).reshape(98, 8) # 1 / Å^2
+# fmt: on
+
+
+def get_atomic_scattering_parameters(element, unit=None):
+ """Return the eight atomic scattering parameters a_1-4, b_1-4 for
+ elements with atomic numbers Z = 1-98 from Table 12.1 in
+ [DeGraef2007]_, which are themselves from [Doyle1968]_ and
+ [Smith1962]_.
+
+ Parameters
+ ----------
+ element : int or str
+ Element to return scattering parameters for. Either one-two
+ letter string or integer atomic number.
+ unit : str, optional
+ Either "nm" or "Å"/"A". Whether to return parameters in terms
+ of Å^-2 or nm^-2. If None (default), Å^-2 is used.
+
+ Returns
+ -------
+ a : np.ndarray
+ The four atomic scattering parameters a_1-4.
+ b : np.ndarray
+ The four atomic scattering parameters b_1-4.
+
+ References
+ ----------
+ .. [DeGraef2007] M. De Graef, M. E. McHenry, "Structure of\
+ Materials," Cambridge University Press (2007).
+ .. [Doyle1968] P. A. Doyle, P. S. Turner, "Relativistic Hartree-Fock
+ X-ray and electron scattering factors," *Acta Cryst.* **24**
+ (1968), doi: https://doi.org/10.1107/S0567739468000756.
+ .. [Smith1962] G. Smith, R. Burge, "The analytical representation
+ of atomic scattering amplitudes for electrons," *Acta Cryst.*
+ **A15** (1962), doi: https://doi.org/10.1107/S0365110X62000481.
+ """
+ if isinstance(element, str):
+ element_id = get_element_id_from_string(element)
+ else:
+ element_id = int(element)
+
+ factor = 1 # Ångstrøm
+ if unit is not None:
+ if unit.lower() == "nm":
+ factor = 1e-2
+
+ a = ATOMIC_SCATTERING_PARAMETERS_DOYLETURNER[element_id - 1, :4] * factor
+ b = ATOMIC_SCATTERING_PARAMETERS_DOYLETURNER[element_id - 1, 4:] * factor
+
+ return a, b
+
+
+def get_element_id_from_string(element_str):
+ """Get periodic element ID for elements Z = 1-98 from one-two letter
+ string.
+
+ Parameters
+ ----------
+ element_str : str
+ One-two letter string.
+
+ Returns
+ -------
+ element_id : int
+ Integer ID in the periodic table of elements.
+ """
+ element2periodic = dict(zip(ELEMENTS[:N_ELEMENTS], np.arange(1, N_ELEMENTS)))
+ element_id = element2periodic[element_str.lower()]
+ return element_id
diff --git a/diffsims/structure_factor/structure_factor.py b/diffsims/structure_factor/structure_factor.py
new file mode 100644
index 00000000..946c5098
--- /dev/null
+++ b/diffsims/structure_factor/structure_factor.py
@@ -0,0 +1,232 @@
+# -*- coding: utf-8 -*-
+# Copyright 2017-2020 The diffsims developers
+#
+# This file is part of diffsims.
+#
+# diffsims is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# diffsims is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with diffsims. If not, see .
+
+from diffpy.structure.symmetryutilities import expandPosition, SymmetryConstraints
+import numpy as np
+from scipy.constants import c, e, h, physical_constants
+
+from diffsims.structure_factor.atomic_scattering_factor import (
+ get_kinematical_atomic_scattering_factor,
+ get_doyleturner_atomic_scattering_factor,
+)
+
+
+rest_mass = physical_constants["atomic unit of mass"][0]
+# rest_mass = 9.109383709e-31 # Used by EMsoft
+
+
+def find_asymmetric_positions(positions, space_group):
+ """Return the asymmetric atom positions among a set of positions
+ when considering symmetry operations defined by a space group.
+
+ Parameters
+ ----------
+ positions : list
+ A list of cartesian atom positions.
+ space_group : diffpy.structure.spacegroups.SpaceGroup
+ Space group describing the symmetry operations.
+
+ Returns
+ -------
+ np.ndarray
+ Asymmetric atom positions.
+ """
+ asymmetric_positions = SymmetryConstraints(space_group, positions).corepos
+ return [
+ np.array([np.allclose(xyz, asym_xyz) for xyz in positions])
+ for asym_xyz in asymmetric_positions
+ ][0]
+
+
+def get_kinematical_structure_factor(phase, hkl, scattering_parameter):
+ """Return the kinematical (X-ray) structure factor for a given family
+ of Miller indices.
+
+ Assumes structure's lattice parameters and Debye-Waller factors are
+ expressed in Ångströms.
+
+ This function is adapted from EMsoft.
+
+ Parameters
+ ----------
+ phase : orix.crystal_map.phase_list.Phase
+ A phase container with a crystal structure and a space and point
+ group describing the allowed symmetry operations.
+ hkl : np.ndarray or list
+ Miller indices.
+ scattering_parameter : float
+ Scattering parameter for these Miller indices.
+
+ Returns
+ -------
+ structure_factor : float
+ Structure factor F.
+ """
+ # Initialize real and imaginary parts of the structure factor
+ structure_factor = 0 + 0j
+
+ structure = phase.structure
+ space_group = phase.space_group
+
+ # Loop over asymmetric unit
+ asymmetric_positions = find_asymmetric_positions(structure.xyz, space_group)
+ for is_asymmetric, atom in zip(asymmetric_positions, structure):
+ if not is_asymmetric:
+ continue
+
+ # Get atomic scattering factor for this atom
+ f = get_kinematical_atomic_scattering_factor(atom, scattering_parameter)
+
+ # Loop over all atoms in the orbit
+ equiv_pos = expandPosition(spacegroup=space_group, xyz=atom.xyz)[0]
+ for xyz in equiv_pos:
+ arg = 2 * np.pi * np.sum(hkl * xyz)
+ structure_factor += f * (np.cos(arg) - (np.sin(arg) * 1j))
+
+ return structure_factor.real
+
+
+def get_doyleturner_structure_factor(
+ phase,
+ hkl,
+ scattering_parameter,
+ voltage,
+ return_parameters=False,
+):
+ """Return the structure factor for a given family of Miller indices
+ using Doyle-Turner atomic scattering parameters [Doyle1968]_.
+
+ Assumes structure's lattice parameters and Debye-Waller factors are
+ expressed in Ångströms.
+
+ This function is adapted from EMsoft.
+
+ Parameters
+ ----------
+ phase : orix.crystal_map.phase_list.Phase
+ A phase container with a crystal structure and a space and point
+ group describing the allowed symmetry operations.
+ hkl : np.ndarray or list
+ Miller indices.
+ scattering_parameter : float
+ Scattering parameter for these Miller indices.
+ voltage : float
+ Beam energy in V.
+ return_parameters : bool, optional
+ Whether to return a set of parameters derived from the
+ calculation as a dictionary. Default is False.
+
+ Returns
+ -------
+ structure_factor : float
+ Structure factor F.
+ params : dict
+ A dictionary with (key, item) (str, float) of parameters derived
+ from the calculation. Only returned if `return_parameters=True`.
+ """
+ structure = phase.structure
+ space_group = phase.space_group
+
+ # Initialize real and imaginary parts of the structure factor
+ structure_factor = 0 + 0j
+
+ # Get unit cell volume for the atomic scattering factor calculation
+ unit_cell_volume = structure.lattice.volume
+
+ # Loop over all atoms in the asymmetric unit
+ asymmetric_positions = find_asymmetric_positions(structure.xyz, space_group)
+ for is_asymmetric, atom in zip(asymmetric_positions, structure):
+ if not is_asymmetric:
+ continue
+
+ # Get atomic scattering factor for this atom
+ f = get_doyleturner_atomic_scattering_factor(
+ atom, scattering_parameter, unit_cell_volume
+ )
+
+ # Loop over all atoms in the orbit
+ equiv_pos = expandPosition(spacegroup=space_group, xyz=atom.xyz)[0]
+ for xyz in equiv_pos:
+ arg = 2 * np.pi * np.sum(hkl * xyz)
+ structure_factor += f * np.exp(-arg * 1j)
+
+ # Relativistic correction factor of wavelength
+ gamma_relcor = 1 + (2 * e * 0.5 * voltage / rest_mass / (c ** 2))
+ # Mean inner potential
+ v_mod = abs(structure_factor) * gamma_relcor
+ v_phase = np.arctan2(structure_factor.imag, structure_factor.real)
+ v_g = v_mod * np.exp(v_phase * 1j)
+ pre = 2 * (rest_mass * e / h ** 2) * 1e-18
+ structure_factor = (pre * v_g).real
+
+ if return_parameters:
+ params = {
+ "gamma_relcor": gamma_relcor,
+ "v_mod": v_mod,
+ "v_phase": v_phase,
+ "v_g": v_g,
+ }
+ return structure_factor, params
+ else:
+ return structure_factor
+
+
+def get_refraction_corrected_wavelength(phase, voltage):
+ """Return the refraction corrected relativistic electron wavelength
+ in Ångströms for a given crystal structure and beam energy in V.
+
+ This function is adapted from EMsoft.
+
+ Parameters
+ ----------
+ phase : orix.crystal_map.phase_list.Phase
+ A phase container with a crystal structure and a space and point
+ group describing the allowed symmetry operations.
+ voltage : float
+ Beam energy in V.
+
+ Returns
+ -------
+ wavelength : float
+ Refraction corrected relativistic electron wavelength in
+ Ångströms.
+ """
+ temp1 = 1e9 * h / np.sqrt(2 * rest_mass * e)
+ temp2 = e * 0.5 * voltage / rest_mass / (c ** 2)
+
+ # Relativistic correction factor (known as gamma). (This is used by EMsoft but not
+ # here, for now.)
+ # gamma_relcor = 1 + (2 * temp2)
+
+ # Relativistic acceleration voltage
+ psi_hat = voltage * (1 + temp2)
+
+ # Compute the electron wavelength in nm
+ hkl = np.zeros(3, dtype=int)
+ scattering_parameter = 0
+ _, params = get_doyleturner_structure_factor(
+ phase, hkl, scattering_parameter, voltage, return_parameters=True
+ )
+ v_mod = params["v_mod"]
+ psi_hat += v_mod
+ wavelength = temp1 / np.sqrt(psi_hat)
+
+ # Interaction constant sigma (this is used by EMsoft but not here, for now)
+ # sigma = 1e-18 * (2 * np.pi * rest_mass * gamma_relcor * e * wavelength) / h ** 2
+
+ return wavelength
diff --git a/diffsims/tests/conftest.py b/diffsims/tests/conftest.py
index 57a43c63..70168058 100644
--- a/diffsims/tests/conftest.py
+++ b/diffsims/tests/conftest.py
@@ -17,20 +17,18 @@
# along with diffsims. If not, see .
import pytest
-import diffpy.structure
-import numpy as np
-from transforms3d.euler import euler2mat
+from diffpy.structure import Atom, Lattice, Structure
+from orix.crystal_map import Phase
-from diffsims.libraries.vector_library import DiffractionVectorLibrary
from diffsims.generators.diffraction_generator import DiffractionGenerator
@pytest.fixture
def default_structure():
"""An atomic structure represented using diffpy """
- latt = diffpy.structure.lattice.Lattice(3, 3, 5, 90, 90, 120)
- atom = diffpy.structure.atom.Atom(atype="Ni", xyz=[0, 0, 0], lattice=latt)
- hexagonal_structure = diffpy.structure.Structure(atoms=[atom], lattice=latt)
+ latt = Lattice(3, 3, 5, 90, 90, 120)
+ atom = Atom(atype="Ni", xyz=[0, 0, 0], lattice=latt)
+ hexagonal_structure = Structure(atoms=[atom], lattice=latt)
return hexagonal_structure
@@ -38,3 +36,47 @@ def default_structure():
def default_simulator():
accelerating_voltage = 300
return DiffractionGenerator(accelerating_voltage)
+
+
+@pytest.fixture
+def nickel_phase():
+ return Phase(
+ name="nickel",
+ space_group=225,
+ structure=Structure(
+ lattice=Lattice(3.5236, 3.5236, 3.5236, 90, 90, 90),
+ atoms=[Atom(xyz=[0, 0, 0], atype="Ni", Uisoequiv=0.006332)]
+ )
+ )
+
+
+@pytest.fixture
+def ferrite_phase():
+ return Phase(
+ name="ferrite",
+ space_group=229,
+ structure=Structure(
+ lattice=Lattice(2.8665, 2.8665, 2.8665, 90, 90, 90),
+ atoms=[
+ Atom(xyz=[0, 0, 0], atype="Fe", Uisoequiv=0.006332),
+ Atom(xyz=[0.5, 0.5, 0.5], atype="Fe", Uisoequiv=0.006332),
+ ]
+ )
+ )
+
+
+@pytest.fixture
+def silicon_carbide_phase():
+ """Silicon Carbide 4H polytype (hexagonal, space group 186)."""
+ return Phase(
+ space_group=186,
+ structure=Structure(
+ lattice=Lattice(3.073, 3.073, 10.053, 90, 90, 120),
+ atoms=[
+ Atom(atype="Si", xyz=[0, 0, 0]),
+ Atom(atype="Si", xyz=[0.33, 0.667, 0.25]),
+ Atom(atype="C", xyz=[0, 0, 0.188]),
+ Atom(atype="C", xyz=[0.333, 0.667, 0.438]),
+ ],
+ )
+ )
diff --git a/diffsims/tests/test_crystallography/__init__.py b/diffsims/tests/test_crystallography/__init__.py
new file mode 100644
index 00000000..063cffcc
--- /dev/null
+++ b/diffsims/tests/test_crystallography/__init__.py
@@ -0,0 +1,17 @@
+# -*- coding: utf-8 -*-
+# Copyright 2017-2020 The diffsims developers
+#
+# This file is part of diffsims.
+#
+# diffsims is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# diffsims is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with diffsims. If not, see .
diff --git a/diffsims/tests/test_crystallography/test_reciprocal_lattice_point.py b/diffsims/tests/test_crystallography/test_reciprocal_lattice_point.py
new file mode 100644
index 00000000..f144a812
--- /dev/null
+++ b/diffsims/tests/test_crystallography/test_reciprocal_lattice_point.py
@@ -0,0 +1,286 @@
+# -*- coding: utf-8 -*-
+# Copyright 2017-2020 The diffsims developers
+#
+# This file is part of diffsims.
+#
+# diffsims is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# diffsims is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with diffsims. If not, see .
+
+from diffsims.crystallography import ReciprocalLatticePoint
+import numpy as np
+from orix.crystal_map import Phase
+from orix.vector import Vector3d
+import pytest
+
+
+class TestReciprocalLatticePoint:
+ @pytest.mark.parametrize(
+ "hkl", [[[1, 1, 1], [2, 0, 0]], np.array([[1, 1, 1], [2, 0, 0]])]
+ )
+ def test_init_rlp(self, nickel_phase, hkl):
+ rlp = ReciprocalLatticePoint(phase=nickel_phase, hkl=hkl)
+ assert rlp.phase.name == nickel_phase.name
+ assert isinstance(rlp.hkl, Vector3d)
+ assert rlp.structure_factor[0] is None
+ assert rlp.theta[0] is None
+ assert rlp.size == 2
+ assert rlp.shape == (2, 3)
+ assert rlp.hkl[0].shape == (1,)
+ assert rlp._hkldata[0].shape == (3,)
+ assert np.issubdtype(rlp.hkl.data.dtype, int)
+
+ @pytest.mark.parametrize("min_dspacing, desired_size", [(2, 9), (1, 19), (0.5, 83)])
+ def test_init_from_min_dspacing(self, ferrite_phase, min_dspacing, desired_size):
+ assert (
+ ReciprocalLatticePoint.from_min_dspacing(
+ phase=ferrite_phase, min_dspacing=min_dspacing
+ ).size
+ == desired_size
+ )
+
+ @pytest.mark.parametrize(
+ "highest_hkl, desired_highest_hkl, desired_lowest_hkl, desired_size",
+ [
+ ([3, 3, 3], [3, 3, 3], [1, 0, 0], 19),
+ ([3, 4, 0], [3, 4, 0], [0, 4, 0], 13),
+ ([4, 3, 0], [4, 3, 0], [1, 0, 0], 13),
+ ],
+ )
+ def test_init_from_highest_hkl(
+ self,
+ silicon_carbide_phase,
+ highest_hkl,
+ desired_highest_hkl,
+ desired_lowest_hkl,
+ desired_size,
+ ):
+ rlp = ReciprocalLatticePoint.from_highest_hkl(
+ phase=silicon_carbide_phase, highest_hkl=highest_hkl
+ )
+ assert np.allclose(rlp[0]._hkldata, desired_highest_hkl)
+ assert np.allclose(rlp[-1]._hkldata, desired_lowest_hkl)
+ assert rlp.size == desired_size
+
+ def test_repr(self, ferrite_phase):
+ rlp = ReciprocalLatticePoint.from_min_dspacing(ferrite_phase, min_dspacing=2)
+ assert repr(rlp) == (
+ f"ReciprocalLatticePoint (9,)\n"
+ f"Phase: ferrite (m-3m)\n"
+ "[[2 2 2]\n [2 2 1]\n [2 2 0]\n [2 1 1]\n [2 1 0]\n [2 0 0]\n [1 1 1]\n"
+ " [1 1 0]\n [1 0 0]]"
+ )
+
+ def test_get_item(self, ferrite_phase):
+ rlp = ReciprocalLatticePoint.from_min_dspacing(
+ phase=ferrite_phase, min_dspacing=1.5
+ )
+ rlp.calculate_structure_factor()
+ rlp.calculate_theta(voltage=20e3)
+
+ assert rlp[0].size == 1
+ assert rlp[:2].size == 2
+ assert np.allclose(rlp[5:7].hkl.data, rlp.hkl[5:7].data)
+
+ assert np.allclose(rlp[10:13].structure_factor, rlp.structure_factor[10:13])
+ assert np.allclose(rlp[20:23].theta, rlp.theta[20:23])
+
+ assert rlp.phase.space_group.number == rlp[0].phase.space_group.number
+ assert rlp.phase.point_group.name == rlp[10:15].phase.point_group.name
+ assert np.allclose(
+ rlp.phase.structure.lattice.abcABG(),
+ rlp[20:23].phase.structure.lattice.abcABG(),
+ )
+
+ def test_get_hkl(self, silicon_carbide_phase):
+ rlp = ReciprocalLatticePoint.from_min_dspacing(
+ silicon_carbide_phase, min_dspacing=3
+ )
+ assert np.allclose(rlp.h, [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0])
+ assert np.allclose(rlp.k, [1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0])
+ assert np.allclose(rlp.l, [4, 3, 2, 1, 0, 4, 3, 2, 0, 4, 3, 2])
+
+ def test_multiplicity(self, nickel_phase, silicon_carbide_phase):
+ assert np.allclose(
+ ReciprocalLatticePoint.from_min_dspacing(
+ phase=nickel_phase, min_dspacing=1
+ ).multiplicity,
+ # fmt: off
+ np.array([
+ 8, 24, 24, 24, 12, 24, 48, 48, 24, 24, 48, 24, 24, 24, 6, 8, 24,
+ 24, 12, 24, 48, 24, 24, 24, 6, 8, 24, 12, 24, 24, 6, 8, 12, 6
+ ])
+ # fmt: on
+ )
+ assert np.allclose(
+ ReciprocalLatticePoint.from_min_dspacing(
+ phase=silicon_carbide_phase, min_dspacing=1
+ ).multiplicity,
+ # fmt: off
+ np.array([
+ 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
+ 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 6,
+ 6, 6, 6, 6, 6, 6, 6, 6, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
+ 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 6, 6, 6, 6, 6, 6, 6, 6, 6, 12,
+ 12, 12, 12, 12, 12, 12, 12, 12, 12, 6, 6, 6, 6, 6, 6, 6, 6, 6, 1, 1,
+ 1, 1, 1, 1, 1, 1
+ ])
+ # fmt: on
+ )
+
+ def test_gspacing_dspacing_scattering_parameter(self, ferrite_phase):
+ rlp = ReciprocalLatticePoint.from_min_dspacing(
+ phase=ferrite_phase, min_dspacing=2
+ )
+ # fmt: off
+ assert np.allclose(
+ rlp.gspacing,
+ np.array([
+ 1.2084778, 1.04657248, 0.98671799, 0.85452285, 0.78006907, 0.69771498,
+ 0.604238, 0.493359, 0.34885749
+ ])
+ )
+ assert np.allclose(
+ rlp.dspacing,
+ np.array([
+ 0.82748727, 0.9555, 1.01346079, 1.17024372, 1.28193777, 1.43325,
+ 1.65497455, 2.02692159, 2.8665
+ ])
+ )
+ assert np.allclose(
+ rlp.scattering_parameter,
+ np.array([
+ 0.6042389, 0.52328624, 0.493359, 0.42726142, 0.39003453, 0.34885749,
+ 0.30211945, 0.2466795, 0.17442875
+ ])
+ )
+ # fmt: on
+
+ @pytest.mark.parametrize(
+ "space_group, hkl, centering, desired_allowed",
+ [
+ (224, [[1, 1, -1], [-2, 2, 2], [3, 1, 0]], "P", [True, True, True]),
+ (230, [[1, 1, -1], [-2, 2, 2], [3, 1, 0]], "I", [False, True, True]),
+ (225, [[1, 1, 1], [2, 2, 1], [-3, 3, 3]], "F", [True, False, True]),
+ (167, [[1, 2, 3], [2, 2, 3], [-3, 2, 4]], "H", [False, True, True]), # R
+ (68, [[1, 2, 3], [2, 2, 3], [-2, 2, 4]], "C", [False, True, True]),
+ (41, [[1, 2, 2], [1, 1, -1], [1, 1, 2]], "A", [True, True, False]),
+ ],
+ )
+ def test_allowed(self, space_group, hkl, centering, desired_allowed):
+ rlp = ReciprocalLatticePoint(phase=Phase(space_group=space_group), hkl=hkl)
+ assert rlp.phase.space_group.short_name[0] == centering
+ assert np.allclose(rlp.allowed, desired_allowed)
+
+ def test_allowed_b_centering(self):
+ """B centering will never fire since no diffpy.structure space group has
+ 'B' first in the name.
+ """
+ phase = Phase(space_group=15)
+ phase.space_group.short_name = "B"
+ rlp = ReciprocalLatticePoint(
+ phase=phase, hkl=[[1, 2, 2], [1, 1, -1], [1, 1, 2]]
+ )
+ assert np.allclose(rlp.allowed, [False, True, False])
+
+ def test_allowed_raises(self, silicon_carbide_phase):
+ with pytest.raises(NotImplementedError):
+ _ = ReciprocalLatticePoint.from_min_dspacing(
+ phase=silicon_carbide_phase, min_dspacing=1
+ ).allowed
+
+ def test_unique(self, ferrite_phase):
+ hkl = [[-1, -1, -1], [1, 1, 1], [1, 0, 0], [0, 0, 1]]
+ rlp = ReciprocalLatticePoint(phase=ferrite_phase, hkl=hkl)
+ assert isinstance(rlp.unique(), ReciprocalLatticePoint)
+ assert np.allclose(rlp.unique(use_symmetry=False)._hkldata, hkl)
+ assert np.allclose(
+ rlp.unique(use_symmetry=True)._hkldata, [[1, 1, 1], [1, 0, 0]]
+ )
+
+ def test_symmetrise(self):
+ assert np.allclose(
+ ReciprocalLatticePoint(phase=Phase(space_group=225), hkl=[1, 1, 1])
+ .symmetrise()
+ ._hkldata,
+ np.array(
+ [
+ [1, 1, 1],
+ [-1, 1, 1],
+ [-1, -1, 1],
+ [1, -1, 1],
+ [1, -1, -1],
+ [1, 1, -1],
+ [-1, 1, -1],
+ [-1, -1, -1],
+ ]
+ ),
+ )
+ rlp2, multiplicity = ReciprocalLatticePoint(
+ phase=Phase(space_group=186), hkl=[2, 2, 0]
+ ).symmetrise(return_multiplicity=True)
+ assert multiplicity == 12
+ assert np.allclose(
+ rlp2._hkldata,
+ [
+ [2, 2, 0],
+ [-2, 0, 0],
+ [0, -2, 0],
+ [-2, -2, 0],
+ [2, 0, 0],
+ [0, 2, 0],
+ [-2, 2, 0],
+ [0, -2, 0],
+ [2, 0, 0],
+ [2, -2, 0],
+ [0, 2, 0],
+ [-2, 0, 0],
+ ],
+ )
+ rlp3 = ReciprocalLatticePoint(
+ phase=Phase(space_group=186), hkl=[2, 2, 0]
+ ).symmetrise(antipodal=False)
+ assert np.allclose(
+ rlp3._hkldata,
+ [[2, 2, 0], [-2, 0, 0], [0, -2, 0], [-2, -2, 0], [2, 0, 0], [0, 2, 0]],
+ )
+
+ @pytest.mark.parametrize(
+ "method, voltage, hkl, desired_factor",
+ [
+ ("kinematical", None, [1, 1, 0], 35.783295),
+ (None, None, [1, 1, 0], 35.783295),
+ ("doyleturner", 20e3, [[2, 0, 0], [1, 1, 0]], [5.581302, 8.096651]),
+ ],
+ )
+ def test_calculate_structure_factor(
+ self, ferrite_phase, method, voltage, hkl, desired_factor
+ ):
+ rlp = ReciprocalLatticePoint(phase=ferrite_phase, hkl=hkl)
+ rlp.calculate_structure_factor(method=method, voltage=voltage)
+ assert np.allclose(rlp.structure_factor, desired_factor)
+
+ def test_calculate_structure_factor_raises(self, ferrite_phase):
+ rlp = ReciprocalLatticePoint(phase=ferrite_phase, hkl=[1, 0, 0])
+ with pytest.raises(ValueError, match="method=man must be among"):
+ rlp.calculate_structure_factor(method="man")
+ with pytest.raises(ValueError, match="'voltage' parameter must be set when"):
+ rlp.calculate_structure_factor(method="doyleturner")
+
+ @pytest.mark.parametrize(
+ "voltage, hkl, desired_theta",
+ [(20e3, [1, 1, 1], 0.00259284), (200e3, [2, 0, 0], 0.00087484)],
+ )
+ def test_calculate_theta(self, ferrite_phase, voltage, hkl, desired_theta):
+ rlp = ReciprocalLatticePoint(phase=ferrite_phase, hkl=hkl)
+ rlp.calculate_theta(voltage=voltage)
+ assert np.allclose(rlp.theta, desired_theta)
diff --git a/diffsims/tests/test_generators/test_diffraction_generator.py b/diffsims/tests/test_generators/test_diffraction_generator.py
index a9d109a0..3fdb3ec3 100644
--- a/diffsims/tests/test_generators/test_diffraction_generator.py
+++ b/diffsims/tests/test_generators/test_diffraction_generator.py
@@ -26,6 +26,7 @@
AtomicDiffractionGenerator,
)
import diffpy.structure
+from diffsims.utils.shape_factor_models import linear, binary, sinc
@pytest.fixture(params=[(300)])
@@ -124,19 +125,25 @@ def test_appropriate_intensities(self, diffraction_calculator, local_structure):
)
assert np.all(smaller)
- @pytest.mark.parametrize("string", ["linear", "binary"])
- def test_shape_factor_strings(
- self, diffraction_calculator, local_structure, string
- ):
+ @pytest.mark.parametrize("model", [None, linear, binary, sinc])
+ def test_shape_factor_strings(self, diffraction_calculator, local_structure, model):
_ = diffraction_calculator.calculate_ed_data(
- local_structure, 2, shape_factor_model=string
+ local_structure, 2, shape_factor_model=model
)
def test_shape_factor_custom(self, diffraction_calculator, local_structure):
def local_excite(excitation_error, maximum_excitation_error, t):
return (np.sin(t) * excitation_error) / maximum_excitation_error
- _ = diffraction_calculator.calculate_ed_data(local_structure, 2,shape_factor_model=local_excite, t=0.2)
+ t1 = diffraction_calculator.calculate_ed_data(
+ local_structure, 2, shape_factor_model=local_excite, t=0.2
+ )
+ t2 = diffraction_calculator.calculate_ed_data(
+ local_structure, 2, shape_factor_model=local_excite, t=0.4
+ )
+
+ # softly makes sure the two sims are different
+ assert np.sum(t1.intensities) != np.sum(t2.intensities)
def test_calculate_profile_class(self, local_structure, diffraction_calculator):
# tests the non-hexagonal (cubic) case
diff --git a/diffsims/tests/test_generators/test_rotation_list_generator.py b/diffsims/tests/test_generators/test_rotation_list_generator.py
index 062474d3..0e8a0dea 100644
--- a/diffsims/tests/test_generators/test_rotation_list_generator.py
+++ b/diffsims/tests/test_generators/test_rotation_list_generator.py
@@ -40,12 +40,16 @@ def test_get_grid(grid):
assert len(grid) > 0
assert isinstance(grid[0], tuple)
+
def test_get_grid_around_beam_direction():
- grid = get_grid_around_beam_direction((0,90,0),resolution=2,angular_range=(0,9))
+ grid = get_grid_around_beam_direction(
+ (0, 90, 0), resolution=2, angular_range=(0, 9)
+ )
assert isinstance(grid, list)
assert isinstance(grid[0], tuple)
- assert len(grid) == 5 # should have 0,2,4,6 and 8
- assert np.allclose([x[1] for x in grid],90) #taking z to y
+ assert len(grid) == 5 # should have 0,2,4,6 and 8
+ assert np.allclose([x[1] for x in grid], 90) # taking z to y
+
@pytest.mark.parametrize(
"crystal_system",
diff --git a/diffsims/tests/test_structure_factor/__init__.py b/diffsims/tests/test_structure_factor/__init__.py
new file mode 100644
index 00000000..063cffcc
--- /dev/null
+++ b/diffsims/tests/test_structure_factor/__init__.py
@@ -0,0 +1,17 @@
+# -*- coding: utf-8 -*-
+# Copyright 2017-2020 The diffsims developers
+#
+# This file is part of diffsims.
+#
+# diffsims is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# diffsims is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with diffsims. If not, see .
diff --git a/diffsims/tests/test_structure_factor/test_atomic_scattering_factor.py b/diffsims/tests/test_structure_factor/test_atomic_scattering_factor.py
new file mode 100644
index 00000000..61861c9b
--- /dev/null
+++ b/diffsims/tests/test_structure_factor/test_atomic_scattering_factor.py
@@ -0,0 +1,96 @@
+# -*- coding: utf-8 -*-
+# Copyright 2017-2020 The diffsims developers
+#
+# This file is part of diffsims.
+#
+# diffsims is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# diffsims is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with diffsims. If not, see .
+
+from diffpy.structure import Atom, Structure, Lattice
+import numpy as np
+import pytest
+
+from diffsims.structure_factor import (
+ get_kinematical_atomic_scattering_factor,
+ get_doyleturner_atomic_scattering_factor,
+)
+
+
+@pytest.mark.parametrize(
+ "element, occupancy, displacement_factor, scattering_parameter, desired_factor",
+ [
+ ("fe", 1, 0.5, 0.17442875, 6.306987),
+ ("fe", 1, 0.5, 0.90635835, 5.911219e-14),
+ ("ni", 1, 0.5, 0.14190033, 10.964046),
+ ("ni", 0.5, 0.5, 0.14190033, 5.482023),
+ ],
+)
+def test_get_kinematical_atomic_scattering_factor(
+ element, occupancy, displacement_factor, scattering_parameter, desired_factor
+):
+ atom = Atom(atype=element, occupancy=occupancy, Uisoequiv=displacement_factor)
+ assert np.allclose(
+ get_kinematical_atomic_scattering_factor(
+ atom=atom,
+ scattering_parameter=scattering_parameter,
+ ),
+ desired_factor,
+ )
+
+
+nickel = Structure(
+ lattice=Lattice(3.5236, 3.5236, 3.5236, 90, 90, 90),
+ atoms=[Atom(xyz=[0, 0, 0], atype="Ni")],
+)
+ferrite = Structure(
+ lattice=Lattice(2.8665, 2.8665, 2.8665, 90, 90, 90),
+ atoms=[Atom(xyz=[0, 0, 0], atype="Fe")],
+)
+# Silicon Carbide 4H polytype (hexagonal, space group 186)
+sic4h = Structure(
+ lattice=Lattice(3.073, 3.073, 10.053, 90, 90, 120),
+ atoms=[
+ Atom(atype="Si", xyz=[0, 0, 0]),
+ Atom(atype="Si", xyz=[0.33, 0.667, 0.25]),
+ Atom(atype="C", xyz=[0, 0, 0.188]),
+ Atom(atype="C", xyz=[0.333, 0.667, 0.438]),
+ ],
+)
+
+
+@pytest.mark.parametrize(
+ "structure, displacement_factor, scattering_parameter, desired_factor",
+ [
+ (ferrite, 0.5, 0.17442875, 2.422658),
+ (ferrite, 0.5, 0.90635835, 9.172283e-15),
+ (ferrite, 0, 0.90635835, 1.114444),
+ (nickel, 0.5, 0.14190033, 2.186904),
+ (sic4h, 0, 0.19894559, 1.512963),
+ ],
+)
+def test_get_doyleturner_atomic_scattering_factor(
+ structure,
+ displacement_factor,
+ scattering_parameter,
+ desired_factor,
+):
+ atom = structure[0]
+ atom.Uisoequiv = displacement_factor
+ assert np.allclose(
+ get_doyleturner_atomic_scattering_factor(
+ atom=atom,
+ scattering_parameter=scattering_parameter,
+ unit_cell_volume=structure.lattice.volume,
+ ),
+ desired_factor,
+ )
diff --git a/diffsims/tests/test_structure_factor/test_atomic_scattering_parameter.py b/diffsims/tests/test_structure_factor/test_atomic_scattering_parameter.py
new file mode 100644
index 00000000..b425a4b1
--- /dev/null
+++ b/diffsims/tests/test_structure_factor/test_atomic_scattering_parameter.py
@@ -0,0 +1,64 @@
+# -*- coding: utf-8 -*-
+# Copyright 2017-2020 The diffsims developers
+#
+# This file is part of diffsims.
+#
+# diffsims is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# diffsims is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with diffsims. If not, see .
+
+import numpy as np
+import pytest
+
+from diffsims.structure_factor import (
+ get_atomic_scattering_parameters,
+ get_element_id_from_string,
+)
+
+
+@pytest.mark.parametrize(
+ "element, unit, desired_parameters",
+ [
+ (10, "m", [0.303, 0.720, 0.475, 0.153, 17.640, 5.860, 1.762, 0.266]),
+ ("Fe", None, [2.544, 2.343, 1.759, 0.506, 64.424, 14.880, 2.854, 0.350]),
+ (
+ "Fe",
+ "NM",
+ [0.02544, 0.02343, 0.01759, 0.00506, 0.64424, 0.14880, 0.02854, 0.00350],
+ ),
+ ("bk", "Å", [6.502, 5.478, 2.510, 0.000, 28.375, 4.975, 0.561, 0.0]),
+ ],
+)
+def test_get_atomic_scattering_parameters(element, unit, desired_parameters):
+ a, b = get_atomic_scattering_parameters(element, unit=unit)
+ assert np.allclose(a, desired_parameters[:4])
+ assert np.allclose(b, desired_parameters[4:])
+
+
+def test_get_element_id_from_string():
+ # fmt: off
+ elements = [
+ "H", "he", "li", "be", "b", "c", "n", "o", "f", "ne", "na", "mg", "al",
+ "sI", "p", "s", "cl", "ar", "k", "ca", "sc", "ti", "v", "cr", "mn",
+ "fe", "co", "ni", "cu", "zn", "ga", "ge", "as", "se", "br", "kr", "rb",
+ "sr", "y", "zr", "nb", "mo", "tc", "ru", "rh", "pd", "ag", "cd", "in",
+ "sn", "SB", "te", "i", "xe", "cs", "ba", "la", "ce", "pr", "nd", "pm",
+ "sm", "eu", "gd", "tb", "dy", "ho", "er", "tm", "yb", "lu", "hf", "ta",
+ "w", "re", "os", "ir", "pt", "au", "hg", "tl", "pb", "bi", "po", "at",
+ "rn", "FR", "ra", "ac", "th", "pa", "u", "np", "pu", "am", "cm", "bk",
+ "cf", "es", "fm", "md", "no", "lr", "rf", "db", "sg", "bh", "hs", "mt",
+ "ds", "rg", "cn", "nh", "fl", "mc", "lv", "ts", "OG"
+ ]
+ # fmt: on
+ elements_id = np.arange(98) + 1
+ for i, element in zip(elements_id, elements):
+ assert get_element_id_from_string(element) == i
diff --git a/diffsims/tests/test_structure_factor/test_structure_factor.py b/diffsims/tests/test_structure_factor/test_structure_factor.py
new file mode 100644
index 00000000..618b32a3
--- /dev/null
+++ b/diffsims/tests/test_structure_factor/test_structure_factor.py
@@ -0,0 +1,161 @@
+# -*- coding: utf-8 -*-
+# Copyright 2017-2020 The diffsims developers
+#
+# This file is part of diffsims.
+#
+# diffsims is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# diffsims is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with diffsims. If not, see .
+
+from diffpy.structure.spacegroups import GetSpaceGroup
+from diffpy.structure import Atom, Lattice, Structure
+import numpy as np
+from orix.crystal_map import Phase
+import pytest
+
+from diffsims.structure_factor import (
+ find_asymmetric_positions,
+ get_kinematical_structure_factor,
+ get_doyleturner_structure_factor,
+ get_refraction_corrected_wavelength,
+)
+
+# Debye-Waller factor in Å^-2: Biosequiv = 8 * np.pi ** 2 * Uisoequiv
+nickel = Phase(
+ space_group=225,
+ structure=Structure(
+ lattice=Lattice(3.5236, 3.5236, 3.5236, 90, 90, 90),
+ atoms=[Atom(xyz=[0, 0, 0], atype="Ni", Uisoequiv=0.006332)],
+ ),
+)
+ferrite = Phase(
+ space_group=229,
+ structure=Structure(
+ lattice=Lattice(2.8665, 2.8665, 2.8665, 90, 90, 90),
+ atoms=[
+ Atom(xyz=[0, 0, 0], atype="Fe", Uisoequiv=0.006332),
+ Atom(xyz=[0.5, 0.5, 0.5], atype="Fe", Uisoequiv=0.006332),
+ ],
+ ),
+)
+sic4h = Phase(
+ space_group=186,
+ structure=Structure(
+ lattice=Lattice(3.073, 3.073, 10.053, 90, 90, 120),
+ atoms=[
+ Atom(atype="Si", xyz=[0, 0, 0], Uisoequiv=0.006332),
+ Atom(atype="Si", xyz=[0.33, 0.667, 0.25], Uisoequiv=0.006332),
+ Atom(atype="C", xyz=[0, 0, 0.188], Uisoequiv=0.006332),
+ Atom(atype="C", xyz=[0.333, 0.667, 0.438], Uisoequiv=0.006332),
+ ],
+ ),
+)
+
+
+@pytest.mark.parametrize(
+ "positions, space_group, desired_asymmetric",
+ [
+ ([[0, 0, 0], [0, 0.5, 0.5]], 229, [True, False]),
+ ([[0.0349, 0.3106, 0.2607], [0.0349, 0.3106, 0.2607]], 229, [True, True]),
+ ],
+)
+def test_find_asymmetric_positions(positions, space_group, desired_asymmetric):
+ assert np.allclose(
+ find_asymmetric_positions(
+ positions=positions, space_group=GetSpaceGroup(space_group)
+ ),
+ desired_asymmetric,
+ )
+
+
+@pytest.mark.parametrize(
+ "phase, hkl, scattering_parameter, desired_factor",
+ [
+ (nickel, [1, 1, 1], 0.245778, 79.665427),
+ (ferrite, [1, 1, 0], 0.2466795, 35.783295),
+ (ferrite, [1, 1, 1], 0.2466795, 0),
+ (sic4h, [0, 0, 4], 0.198945, 19.027004),
+ ],
+)
+def test_get_kinematical_structure_factor(
+ phase, hkl, scattering_parameter, desired_factor
+):
+ assert np.allclose(
+ get_kinematical_structure_factor(
+ phase=phase, hkl=hkl, scattering_parameter=scattering_parameter
+ ),
+ desired_factor,
+ )
+
+
+@pytest.mark.parametrize(
+ "phase, hkl, scattering_parameter, voltage, desired_factor",
+ [
+ (nickel, [1, 1, 1], 0.245778, 15e3, 8.606363),
+ (ferrite, [1, 1, 0], 0.2466795, 20e3, 8.096651),
+ (ferrite, [1, 1, 1], 0.2466795, 20e3, 0),
+ (sic4h, [0, 0, 4], 0.198945, 200e3, 2.744304),
+ ],
+)
+def test_get_doyleturner_structure_factor(
+ phase, hkl, scattering_parameter, voltage, desired_factor
+):
+ """Tested against EMsoft v5."""
+ assert np.allclose(
+ get_doyleturner_structure_factor(
+ phase=phase,
+ hkl=hkl,
+ scattering_parameter=scattering_parameter,
+ voltage=voltage,
+ ),
+ desired_factor,
+ )
+
+
+@pytest.mark.parametrize(
+ "hkl, scattering_parameter, desired_factor, desired_params",
+ [
+ ([1, 1, 0], 0.2466795, 8.096628, [1.03, 12.17, 1.22e-16, 12.17 + 2.45e-15j]),
+ ([2, 0, 0], 0.348857, 5.581301, [1.03, 8.39, 1.22e-16, 8.39 + 1.02e-15j]),
+ ],
+)
+def test_get_doyleturner_structure_factor_returns(
+ ferrite_phase, hkl, scattering_parameter, desired_factor, desired_params
+):
+ """Tested against EMsoft v5. Only the imaginary part of v_g is off."""
+ sf, params = get_doyleturner_structure_factor(
+ phase=ferrite_phase,
+ hkl=hkl,
+ scattering_parameter=scattering_parameter,
+ voltage=20e3,
+ return_parameters=True,
+ )
+ assert np.allclose(sf, desired_factor)
+ assert np.allclose(list(params.values()), desired_params, atol=1e-2)
+
+
+@pytest.mark.parametrize(
+ "phase, voltage, desired_wavelength",
+ [
+ (nickel, 20e3, 0.008582),
+ (nickel, 200e3, 0.002507),
+ (ferrite, 10e3, 0.012186),
+ (sic4h, 100e3, 0.003701),
+ ],
+)
+def test_get_refraction_corrected_wavelength(phase, voltage, desired_wavelength):
+ """Tested against EMsoft v5."""
+ assert np.allclose(
+ get_refraction_corrected_wavelength(phase=phase, voltage=voltage),
+ desired_wavelength,
+ atol=1e-6,
+ )
diff --git a/diffsims/utils/shape_factor_models.py b/diffsims/utils/shape_factor_models.py
new file mode 100644
index 00000000..e7af7d65
--- /dev/null
+++ b/diffsims/utils/shape_factor_models.py
@@ -0,0 +1,83 @@
+# -*- coding: utf-8 -*-
+# Copyright 2017-2020 The diffsims developers
+#
+# This file is part of diffsims.
+#
+# diffsims is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# diffsims is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with diffsims. If not, see .
+
+import numpy as np
+
+
+def binary(excitation_error, max_excitation_error):
+ """
+ Returns a unit intensity for all reflections
+
+ Parameters
+ ----------
+ excitation_error : array-like or float
+ The distance (reciprocal) from a reflection to the Ewald sphere
+
+ max_excitation_error : float
+ The distance at which a reflection becomes extinct
+
+ Returns
+ -------
+ intensities : array-like or float
+ """
+ return 1
+
+
+def linear(excitation_error, max_excitation_error):
+ """
+ Returns an intensity linearly scaled with by the excitation error
+
+ Parameters
+ ----------
+ excitation_error : array-like or float
+ The distance (reciprocal) from a reflection to the Ewald sphere
+
+ max_excitation_error : float
+ The distance at which a reflection becomes extinct
+
+ Returns
+ -------
+ intensities : array-like or float
+ """
+
+ return 1 - excitation_error / max_excitation_error
+
+
+def sinc(excitation_error, max_excitation_error, minima_number=5):
+ """
+ Returns an intensity with a sinc profile
+
+ Parameters
+ ----------
+ excitation_error : array-like or float
+ The distance (reciprocal) from a reflection to the Ewald sphere
+
+ max_excitation_error : float
+ The distance at which a reflection becomes extinct
+
+ minima_number : int
+ The minima_number'th minima lies at max_excitation_error from 0
+
+ Returns
+ -------
+ intensity : array-like or float
+ """
+
+ num = np.sin(np.pi * minima_number * excitation_error / max_excitation_error)
+ denom = excitation_error
+ return np.abs(np.divide(num, denom, out=np.zeros_like(num), where=denom != 0))