Feature - ZAP map

diffsims/generators/rotation_list_generators.py

@@ -24,11 +24,13 @@
 import warnings
 from itertools import product
 
+from transforms3d.euler import euler2axangle,axangle2euler
+
 from diffsims.utils.rotation_conversion_utils import Euler
 from diffsims.utils.fundemental_zone_utils import get_proper_point_group_string, reduce_to_fundemental_zone
 from diffsims.utils.gridding_utils import create_linearly_spaced_array_in_rzxz,rotate_axangle, \
                                           _create_advanced_linearly_spaced_array_in_rzxz, \
-                                          _get_rotation_to_beam_direction, get_beam_directions, beam_directions_to_euler_angles
+                                          get_beam_directions, beam_directions_to_euler_angles
 
 
 def _returnable_eulers_from_axangle(grid,axis_convention,round_to):
@@ -82,7 +84,7 @@ def get_grid_streographic(crystal_system,resolution,equal='angle'):
     Returns
     -------
     rotation_list : list of tuples
-        List of rotations 
+        List of rotations
     """
     beam_directions_rzxz = beam_directions_to_euler_angles(get_beam_directions(crystal_system,resolution,equal=equal))
     beam_directions_szxz = beam_directions_rzxz.to_AxAngle().to_Euler(axis_convention='szxz') # convert to high speed convention
@@ -131,33 +133,37 @@ def get_local_grid(center, max_rotation, resolution):
     return _returnable_eulers_from_axangle(raw_grid_axangle,'rzxz',round_to=2)
 
 
-def get_grid_around_beam_direction(beam_direction,resolution, angular_range=(0, 360),cubic=False):
+def get_grid_around_beam_direction(beam_rotation,resolution, angular_range=(0, 360)):
     """
     Creates a rotation list of rotations for which the rotation is about given beam direction
 
     Parameters
     ----------
-    beam_direction : [x,y,z]
-        A desired beam direction
+    beam_rotation : tuple
+        A desired beam direction as a rotation (rzxz eulers), usually found via get_rotation_from_z_to_direction
 
     resolution : float
-        The 'resolution' of the grid (degrees)
+        The resolution of the grid (degrees)
 
     angular_range : tuple
         The minimum (included) and maximum (excluded) rotation around the beam direction to be included
 
-    cubic : bool
-        This only works for cubic systems at the present, when False this raises a warning, set to
-        True to supress said warning. The default is False
-
     Returns
     -------
     rotation_list : list of tuples
+
+    Example
+    -------
+    >>> from diffsims.generators.zap_map_generator import get_rotation_from_z_to_direction
+    >>> beam_rotation = get_rotation_from_z_to_direction(structure,[1,1,1])
+    >>> grid = get_grid_around_beam_direction(beam_rotation,1)
     """
 
-    if not cubic:
-        warnings.warn("This code only works for cubic systems at present")
-    rotation_alpha, rotation_beta = _get_rotation_to_beam_direction(beam_direction)
+    beam_rotation = np.deg2rad(beam_rotation)
+    axangle = euler2axangle(beam_rotation[0],beam_rotation[1],beam_rotation[2],'rzxz')
+    euler_szxz = axangle2euler(axangle[0],axangle[1],'szxz') # convert to szxz
+    rotation_alpha, rotation_beta = np.rad2deg(euler_szxz[0]),np.rad2deg(euler_szxz[1])
+
     # see _create_advanced_linearly_spaced_array_in_rzxz for details
     steps_gamma = int(np.ceil((angular_range[1] - angular_range[0])/resolution))
     alpha = np.asarray([rotation_alpha])

diffsims/generators/zap_map_generator.py

@@ -0,0 +1,178 @@
+# -*- coding: utf-8 -*-
+# Copyright 2017-2019 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 <http://www.gnu.org/licenses/>.
+
+import numpy as np
+from transforms3d.euler import axangle2euler
+
+def get_rotation_from_z_to_direction(structure,direction):
+    """
+    Finds the rotation that takes [001] to a given zone axis.
+
+    Parameters
+    ----------
+    structure : diffpy.structure
+        The structure for which a rotation needs to be found
+
+    direction : array like
+        [UVW] direction that the 'z' axis should end up point down.
+
+    Returns
+    -------
+    euler_angles : tuple
+        'rzxz' in degrees
+
+    See Also
+    --------
+    generate_zap_map
+    get_grid_around_beam_direction
+
+    Notes
+    -----
+    This implementation works with an axis arrangement that has +x as left to right,
+    +y as bottom to top and +z as out of the plane of a page. Rotations are counter clockwise
+    as you look from the tip of the axis towards the origin
+    """
+
+    # Case where we don't need a rotation, As axis is [0,0,z] or [0,0,0]
+    if np.dot(direction,[0,0,1]) == np.linalg.norm(direction):
+        return (0,0,0)
+
+    # Normalize our directions
+    cartesian_direction = structure.lattice.cartesian(direction)
+    cartesian_direction = cartesian_direction / np.linalg.norm(cartesian_direction)
+
+    #Find the rotation using cartesian vector geometry
+    rotation_axis = np.cross([0,0,1],cartesian_direction)
+    rotation_angle = np.arccos(np.dot([0,0,1],cartesian_direction))
+    euler = axangle2euler(rotation_axis,rotation_angle,axes='rzxz')
+    return np.rad2deg(euler)
+
+def generate_directional_simulations(structure,simulator,direction_list,reciprocal_radius=1,**kwargs):
+    """
+    Produces simualtion of a structure aligned with certain axes
+
+    Parameters
+    ----------
+    structure : diffpy.structure
+        The structure from which simulations need to be produced
+
+    simulator : DiffractionGenerator
+        The diffraction generator object used to produce the simulations
+
+    direction_list : list of lists
+        A list of [UVW] indicies, eg) [[1,0,0],[1,1,0]]
+
+    reciprocal_radius : float
+        Default to 1
+
+    Returns
+    -------
+    direction_dictionary : dict
+        Keys are zone axes, values are simulations
+    """
+
+    direction_dictionary = {}
+    for direction in direction_list:
+        if np.allclose(direction,0):
+            break
+        rotation_rzxz = get_rotation_from_z_to_direction(structure,direction)
+        simulation = simulator.calculate_ed_data(structure,reciprocal_radius,rotation=rotation_rzxz,**kwargs)
+        direction_dictionary[direction] = simulation
+
+    return direction_dictionary
+
+def corners_to_centroid_and_edge_centers(corners):
+    """
+    Produces the midpoints and center of a trio of corners
+
+    Parameters
+    ----------
+    corners : list of lists
+        Three corners of a streographic triangle
+
+    Returns
+    -------
+    list_of_corners : list
+        Length 7, elements ca,cb,cc,mean,cab,cbc,cac where naming is such that
+        ca is the first corner of the input, and cab is the midpoint between
+        corner a and corner b.
+    """
+    ca,cb,cc = corners[0],corners[1],corners[2]
+    mean = tuple(np.add(np.add(ca,cb),cc))
+    cab  = tuple(np.add(ca,cb))
+    cbc  = tuple(np.add(cb,cc))
+    cac  = tuple(np.add(ca,cc))
+    return [ca,cb,cc,mean,cab,cbc,cac]
+
+def generate_zap_map(structure,simulator,system='cubic',reciprocal_radius=1,density='7',**kwargs):
+    """
+    Produces a number of zone axis patterns for a structure
+
+    Parameters
+    ----------
+    structure : diffpy.structure
+        The structure to be simulated
+
+    simulator : DiffractionGenerator
+        The simulator used to generate the simulations
+
+    system : str
+        'cubic','hexagonal', 'trigonal', 'tetragonal','orthorhombic','monoclinic' - Defaults to 'cubic'
+
+    reciprocal_radius : float
+        The range of reciprocal lattice spots to be included. Default to 1
+
+    density : str
+        '3' for the corners or '7' (corners + midpoints + centroids). Defaults to 7
+
+    **kwargs :
+        keyword arguments to be passed to simulator.calculate_ed_data()
+
+    Returns
+    -------
+    zap_dictionary : dict
+        Keys are zone axes, values are simulations
+
+    Example
+    -------
+    Plot all of the patterns that you have generated
+
+    >>> zap_map = generate_zap_map(structure,simulator,'hexagonal',density='3')
+    >>> for k in zap_map.keys():
+    >>>     pattern = zap_map[k]
+    >>>     pattern.calibration = 4e-3
+    >>>     plt.imshow(pattern.get_diffraction_pattern(),vmax=0.02)
+    >>>     plt.figure()
+
+    """
+
+    corners_dict = {'cubic': [(0, 0, 1), (1, 0, 1), (1, 1, 1)],
+    'hexagonal': [(0, 0, 1), (2, 1, 0), (1, 1, 0)],
+    'orthorhombic': [(0, 0, 1), (1, 0, 0), (0, 1, 0)],
+    'tetragonal': [(0, 0, 1), (1, 0, 0), (1, 1, 0)],
+    'trigonal': [(0, 0, 1), (-1, -2, 0), (1, -1, 0)],
+    'monoclinic': [(0, 0, 1), (0, 1, 0), (0, -1, 0)]}
+
+    if density == '3':
+        direction_list = corners_dict[system]
+    elif density == '7':
+        direction_list  = corners_to_centroid_and_edge_centers(corners_dict[system])
+
+    zap_dictionary = generate_directional_simulations(structure,simulator,direction_list,**kwargs)
+
+    return zap_dictionary

diffsims/tests/conftest.py

@@ -22,16 +22,24 @@
 from transforms3d.euler import euler2mat
 
 from diffsims.libraries.vector_library import DiffractionVectorLibrary
+from diffsims.generators.diffraction_generator import DiffractionGenerator
 
 @pytest.fixture
 def default_structure():
-    """An atomic structure represetned using diffpy
+    """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)
     return hexagonal_structure
 
+@pytest.fixture
+def default_simulator():
+    accelerating_voltage = 300
+    max_excitation_error = 1e-2
+    return DiffractionGenerator(accelerating_voltage,max_excitation_error)
+
+
 @pytest.fixture()
 def random_eulers():
     """ Using [0,360] [0,180] and [0,360] as ranges """

diffsims/tests/test_generators/test_zap_map_generator.py

@@ -0,0 +1,80 @@
+# -*- coding: utf-8 -*-
+# Copyright 2017-2019 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 <http://www.gnu.org/licenses/>.
+
+import pytest
+
+import diffpy.structure
+import numpy as np
+from diffsims.generators.zap_map_generator import get_rotation_from_z_to_direction, generate_zap_map
+
+def test_zero_rotation_cases(default_structure):
+    r_test = get_rotation_from_z_to_direction(default_structure,[0,0,2])
+    r_test_zero = get_rotation_from_z_to_direction(default_structure,[0,0,0])
+    assert r_test == (0,0,0)
+    assert r_test_zero == (0,0,0)
+
+
+class TestOrthonormals:
+
+    @pytest.fixture(params=[(3,3,3),(3,3,4),(3,4,5)])
+    def sample_system(self,request):
+        """Orthonormal structures"""
+        a,b,c = request.param[0],request.param[1],request.param[2]
+        latt = diffpy.structure.lattice.Lattice(a,b,c,90,90,90)
+        atom = diffpy.structure.atom.Atom(atype='Ni',xyz=[0,0,0],lattice=latt)
+        return diffpy.structure.Structure(atoms=[atom],lattice=latt)
+
+    def test_rotation_to__static_x_axis(self,sample_system):
+        r_to_x = get_rotation_from_z_to_direction(sample_system,[1,0,0])
+        assert np.allclose(r_to_x,(90,90,-90))
+
+    def test_rotation_to_static_y_axis(self,sample_system):
+        r_to_y = get_rotation_from_z_to_direction(sample_system,[0,1,0])
+        assert np.allclose(r_to_y,(180,90,-180))
+
+    def test_rotations_to_static_yz(self,sample_system):
+        """ We rotate from z towards y, and compare the results to geometry"""
+        r_to_yz = get_rotation_from_z_to_direction(sample_system,[0,1,1])
+        tan_angle = np.tan(np.deg2rad(r_to_yz[1]))
+        tan_lattice = sample_system.lattice.b / sample_system.lattice.c
+        assert np.allclose(tan_angle,tan_lattice,atol=1e-5)
+
+
+@pytest.mark.parametrize('system',['cubic','hexagonal','trigonal','orthorhombic','tetragonal','monoclinic'])
+def test_zap_map_all_systems(default_structure,default_simulator,system):
+    z_dict = generate_zap_map(default_structure,default_simulator,system=system)
+    assert (0,0,1) in z_dict.keys()
+    assert (0,0,0) not in z_dict.keys()
+
+@pytest.mark.parametrize('density',['3','7'])
+def test_zap_map_density_changes(default_structure,default_simulator,density):
+    """ Checks density arguments are passed correctly """
+    z_dict = generate_zap_map(default_structure,default_simulator,density=density)
+    if density == '3':
+        assert str(len(z_dict.keys())) == '3'
+    elif density == '7':
+        assert len(z_dict.keys()) > 5 #monoclinic case gives 6 rather than 7
+
+def test_zap_map_kwargs(default_structure,default_simulator):
+    z_dict_no_beam = generate_zap_map(default_structure,default_simulator,with_direct_beam=False)
+    z_dict_yes_beam = generate_zap_map(default_structure,default_simulator,with_direct_beam=True)
+    for k in z_dict_no_beam.keys():
+        #both dictionary's have the same keys
+        assert k in z_dict_yes_beam.keys()
+        # no beam has one fewer spots than yes beam
+        assert z_dict_no_beam[k].intensities.shape[0] == z_dict_yes_beam[k].intensities.shape[0] - 1