diff --git a/defdap/__init__.py b/defdap/__init__.py index 000a951..75f37f8 100644 --- a/defdap/__init__.py +++ b/defdap/__init__.py @@ -18,7 +18,8 @@ 'BCC': 'cubic_bcc', 'HCP': 'hexagonal_withca', }, - 'ipf_triangle_convention': 'aztec' + # up or down + 'ipf_triangle_convention': 'up' } anonymous_experiment = Experiment() diff --git a/defdap/crystal.py b/defdap/crystal.py index e6e4fa7..72df6fd 100755 --- a/defdap/crystal.py +++ b/defdap/crystal.py @@ -19,6 +19,7 @@ from defdap import defaults from defdap.quat import Quat +from defdap.crystal_utils import * class Phase(object): @@ -104,78 +105,6 @@ def __init__(self, name, symmetries, vertices, faces): self.vertices = vertices self.faces = faces - # TODO: Move these to the phase class where the lattice parameters - # can be accessed - @staticmethod - def l_matrix(a, b, c, alpha, beta, gamma, convention=None): - """ Construct L matrix based on Page 22 of - Randle and Engle - Introduction to texture analysis""" - l_matrix = np.zeros((3, 3)) - - cos_alpha = np.cos(alpha) - cos_beta = np.cos(beta) - cos_gamma = np.cos(gamma) - - sin_gamma = np.sin(gamma) - - l_matrix[0, 0] = a - l_matrix[0, 1] = b * cos_gamma - l_matrix[0, 2] = c * cos_beta - - l_matrix[1, 1] = b * sin_gamma - l_matrix[1, 2] = c * (cos_alpha - cos_beta * cos_gamma) / sin_gamma - - l_matrix[2, 2] = c * np.sqrt( - 1 + 2 * cos_alpha * cos_beta * cos_gamma - - cos_alpha**2 - cos_beta**2 - cos_gamma**2 - ) / sin_gamma - - # OI/HKL convention - x // [10-10], y // a2 [-12-10] - # TSL convention - x // a1 [2-1-10], y // [01-10] - if convention is None: - convention = defaults['crystal_ortho_conv'] - - if convention.lower() in ['hkl', 'oi']: - # Swap 00 with 11 and 01 with 10 due to how OI orthonormalises - # From Brad Wynne - t1 = l_matrix[0, 0] - t2 = l_matrix[1, 0] - - l_matrix[0, 0] = l_matrix[1, 1] - l_matrix[1, 0] = l_matrix[0, 1] - - l_matrix[1, 1] = t1 - l_matrix[0, 1] = t2 - - elif convention.lower() != 'tsl': - raise ValueError( - f"Unknown convention '{convention}' for orthonormalisation of " - f"crystal structure, can be 'hkl' or 'tsl'" - ) - - # Set small components to 0 - l_matrix[np.abs(l_matrix) < 1e-10] = 0 - - return l_matrix - - @staticmethod - def q_matrix(l_matrix): - """ Construct matrix of reciprocal lattice vectors to transform - plane normals See C. T. Young and J. L. Lytton, J. Appl. Phys., - vol. 43, no. 4, pp. 1408–1417, 1972.""" - a = l_matrix[:, 0] - b = l_matrix[:, 1] - c = l_matrix[:, 2] - - volume = abs(np.dot(a, np.cross(b, c))) - a_star = np.cross(b, c) / volume - b_star = np.cross(c, a) / volume - c_star = np.cross(a, b) / volume - - q_matrix = np.stack((a_star, b_star, c_star), axis=1) - - return q_matrix - over_root2 = np.sqrt(2) / 2 sqrt3over2 = np.sqrt(3) / 2 @@ -327,14 +256,14 @@ def __init__(self, slip_plane, slip_dir, crystal_structure, c_over_a=None): slip_dir_m = convert_idc('mb', dir=slip_dir) # Transformation from crystal to orthonormal coords - l_matrix = CrystalStructure.l_matrix( + l_matrix = create_l_matrix( 1, 1, c_over_a, np.pi / 2, np.pi / 2, np.pi * 2 / 3 ) # Q matrix for transforming planes - qMatrix = CrystalStructure.q_matrix(l_matrix) + q_matrix = create_q_matrix(l_matrix) # Transform into orthonormal basis and then normalise - self.slip_plane = np.matmul(qMatrix, slip_plane_m) + self.slip_plane = np.matmul(q_matrix, slip_plane_m) self.slip_plane /= norm(self.slip_plane) self.slip_dir = np.matmul(l_matrix, slip_dir_m) self.slip_dir /= norm(self.slip_dir) @@ -367,7 +296,7 @@ def slip_plane_label(self): Slip plane label. """ - return '(' + ''.join(map(str_idx, self.plane_idc)) + ')' + return idc_to_string(self.plane_idc, '()') @property def slip_dir_label(self): @@ -379,7 +308,7 @@ def slip_dir_label(self): Slip direction label. """ - return '[' + ''.join(map(str_idx, self.dir_idc)) + ']' + return idc_to_string(self.dir_idc, '[]') def generate_family(self): """Generate the family of slip systems which this system belongs to. @@ -400,11 +329,11 @@ def generate_family(self): if self.crystal_structure.name == 'hexagonal': # Transformation from crystal to orthonormal coords - l_matrix = CrystalStructure.l_matrix( + l_matrix = create_l_matrix( 1, 1, self.c_over_a, np.pi / 2, np.pi / 2, np.pi * 2 / 3 ) # Q matrix for transforming planes - q_matrix = CrystalStructure.q_matrix(l_matrix) + q_matrix = create_q_matrix(l_matrix) # Transform into orthonormal basis plane = np.matmul(q_matrix, convert_idc('mb', plane=plane)) @@ -566,163 +495,3 @@ def print_slip_system_directory(): package_dir, _ = os.path.split(__file__) print("Slip system definition files are stored in directory:") print(f"{package_dir}/slip_systems/") - - -def convert_idc(in_type, *, dir=None, plane=None): - """ - Convert between Miller and Miller-Bravais indices. - - Parameters - ---------- - in_type : str {'m', 'mb'} - Type of indices provided. If 'm' converts from Miller to - Miller-Bravais, opposite for 'mb'. - dir : tuple of int or equivalent, optional - Direction to convert. This OR `plane` must me provided. - plane : tuple of int or equivalent, optional - Plane to convert. This OR `direction` must me provided. - - Returns - ------- - tuple of int - The converted plane or direction. - - """ - if dir is None and plane is None: - raise ValueError("One of either `direction` or `plane` must be " - "provided.") - if dir is not None and plane is not None: - raise ValueError("One of either `direction` or `plane` must be " - "provided, not both.") - - def check_len(val, length): - if len(val) != length: - raise ValueError(f"Vector must have {length} values.") - - if in_type.lower() == 'm': - if dir is None: - # plane M->MB - check_len(plane, 3) - out = np.array(plane)[[0, 1, 0, 2]] - out[2] += plane[1] - out[2] *= -1 - - else: - # direction M->MB - check_len(dir, 3) - u, v, w = dir - out = np.array([2*u-v, 2*v-u, -u-v, 3*w]) / 3 - try: - # Attempt to cast to integers - out = safe_int_cast(out) - except ValueError: - pass - - elif in_type.lower() == 'mb': - if dir is None: - # plane MB->M - check_len(plane, 4) - out = np.array(plane)[[0, 1, 3]] - - else: - # direction MB->M - check_len(dir, 4) - out = np.array(dir)[[0, 1, 3]] - out[[0, 1]] -= dir[2] - - else: - raise ValueError("`inType` must be either 'm' or 'mb'.") - - return tuple(out) - - -def pos_idc(vec): - """ - Return a consistent positive version of a set of indices. - - Parameters - ---------- - vec : tuple of int or equivalent - Indices to convert. - - Returns - ------- - tuple of int - Positive version of indices. - - """ - for idx in vec: - if idx == 0: - continue - if idx > 0: - return tuple(vec) - else: - return tuple(-np.array(vec)) - - -def reduce_idc(vec): - """ - Reduce indices to lowest integers - - Parameters - ---------- - vec : tuple of int or equivalent - Indices to reduce. - - Returns - ------- - tuple of int - The reduced indices. - - """ - return tuple((np.array(vec) / np.gcd.reduce(vec)).astype(np.int8)) - - -def safe_int_cast(vec, tol=1e-3): - """ - Cast a tuple of floats to integers, raising an error if rounding is - over a tolerance. - - Parameters - ---------- - vec : tuple of float or equivalent - Vector to cast. - tol : float - Tolerance above which an error is raised. - - Returns - ------- - tuple of int - - Raises - ------ - ValueError - If the rounding is over the tolerance for any value. - - """ - vec = np.array(vec) - vec_rounded = vec.round() - - if np.any(np.abs(vec - vec_rounded) > tol): - raise ValueError('Rounding too large', np.abs(vec - vec_rounded)) - - return tuple(vec_rounded.astype(np.int8)) - - -def str_idx(idx): - """ - String representation of an index with overbars. - - Parameters - ---------- - idx : int - - Returns - ------- - str - - """ - if not isinstance(idx, (int, np.integer)): - raise ValueError("Index must be an integer.") - - return str(idx) if idx >= 0 else str(-idx) + u'\u0305' diff --git a/defdap/crystal_utils.py b/defdap/crystal_utils.py new file mode 100644 index 0000000..18018b5 --- /dev/null +++ b/defdap/crystal_utils.py @@ -0,0 +1,412 @@ +# Copyright 2024 Mechanics of Microstructures Group +# at The University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from functools import partial +import numpy as np + +from defdap import defaults + +__all__ = [ + 'create_l_matrix', + 'create_q_matrix', + 'convert_idc', + 'equavlent_indicies', + 'project_to_orth', + 'pos_idc', + 'reduce_idc', + 'safe_int_cast', + 'idc_to_string', + 'str_idx', +] + + +def create_l_matrix(a, b, c, alpha, beta, gamma, convention=None): + """ Construct L matrix based on Page 22 of + Randle and Engle - Introduction to texture analysis""" + l_matrix = np.zeros((3, 3)) + + cos_alpha = np.cos(alpha) + cos_beta = np.cos(beta) + cos_gamma = np.cos(gamma) + + sin_gamma = np.sin(gamma) + + l_matrix[0, 0] = a + l_matrix[0, 1] = b * cos_gamma + l_matrix[0, 2] = c * cos_beta + + l_matrix[1, 1] = b * sin_gamma + l_matrix[1, 2] = c * (cos_alpha - cos_beta * cos_gamma) / sin_gamma + + l_matrix[2, 2] = c * np.sqrt( + 1 + 2 * cos_alpha * cos_beta * cos_gamma - + cos_alpha**2 - cos_beta**2 - cos_gamma**2 + ) / sin_gamma + + # OI/HKL convention - x // [10-10], y // a2 [-12-10] + # TSL convention - x // a1 [2-1-10], y // [01-10] + if convention is None: + convention = defaults['crystal_ortho_conv'] + + if convention.lower() in ['hkl', 'oi']: + # Swap 00 with 11 and 01 with 10 due to how OI orthonormalises + # From Brad Wynne + t1 = l_matrix[0, 0] + t2 = l_matrix[1, 0] + + l_matrix[0, 0] = l_matrix[1, 1] + l_matrix[1, 0] = l_matrix[0, 1] + + l_matrix[1, 1] = t1 + l_matrix[0, 1] = t2 + + elif convention.lower() != 'tsl': + raise ValueError( + f"Unknown convention '{convention}' for orthonormalisation of " + f"crystal structure, can be 'hkl' or 'tsl'" + ) + + # Set small components to 0 + l_matrix[np.abs(l_matrix) < 1e-10] = 0 + + return l_matrix + + +def create_q_matrix(l_matrix): + """ Construct matrix of reciprocal lattice vectors to transform + plane normals See C. T. Young and J. L. Lytton, J. Appl. Phys., + vol. 43, no. 4, pp. 1408–1417, 1972.""" + a = l_matrix[:, 0] + b = l_matrix[:, 1] + c = l_matrix[:, 2] + + volume = abs(np.dot(a, np.cross(b, c))) + a_star = np.cross(b, c) / volume + b_star = np.cross(c, a) / volume + c_star = np.cross(a, b) / volume + + q_matrix = np.stack((a_star, b_star, c_star), axis=1) + + return q_matrix + + +def check_len(val, length): + if len(val) != length: + raise ValueError(f"Vector must have {length} values.") + + +def convert_idc(in_type, *, dir=None, plane=None): + """ + Convert between Miller and Miller-Bravais indices. + + Parameters + ---------- + in_type : str {'m', 'mb'} + Type of indices provided. If 'm' converts from Miller to + Miller-Bravais, opposite for 'mb'. + dir : tuple of int or equivalent, optional + Direction to convert. This OR `plane` must me provided. + plane : tuple of int or equivalent, optional + Plane to convert. This OR `direction` must me provided. + + Returns + ------- + tuple of int + The converted plane or direction. + + """ + if dir is None and plane is None: + raise ValueError("One of either `direction` or `plane` must be " + "provided.") + if dir is not None and plane is not None: + raise ValueError("One of either `direction` or `plane` must be " + "provided, not both.") + + if in_type.lower() == 'm': + if dir is None: + # plane M->MB + check_len(plane, 3) + out = np.array(plane)[[0, 1, 0, 2]] + out[2] += plane[1] + out[2] *= -1 + + else: + # direction M->MB + check_len(dir, 3) + u, v, w = dir + out = np.array([2*u-v, 2*v-u, -u-v, 3*w]) / 3 + try: + # Attempt to cast to integers + out = safe_int_cast(out) + except ValueError: + pass + + elif in_type.lower() == 'mb': + if dir is None: + # plane MB->M + check_len(plane, 4) + out = np.array(plane)[[0, 1, 3]] + + else: + # direction MB->M + check_len(dir, 4) + out = np.array(dir)[[0, 1, 3]] + out[[0, 1]] -= dir[2] + + else: + raise ValueError("`inType` must be either 'm' or 'mb'.") + + return tuple(out) + + +def equavlent_indicies( + crystal_symm, + symmetries, + *, + dir=None, + plane=None, + c_over_a=None, + in_type=None +): + if dir is None and plane is None: + raise ValueError("One of either `direction` or `plane` must be " + "provided.") + if dir is not None and plane is not None: + raise ValueError("One of either `direction` or `plane` must be " + "provided, not both.") + if in_type is None: + in_type = 'mb' if crystal_symm == 'hexagonal' else 'm' + + planes = [] + dirs = [] + + if in_type == 'mb': + if dir is None: + check_len(plane, 4) + plane = convert_idc('mb', plane=plane) + if plane is None: + check_len(dir, 4) + dir = convert_idc('mb', dir=dir) + elif in_type != 'm': + raise ValueError("`inType` must be either 'm' or 'mb'.") + + if dir is None: + check_len(plane, 3) + rtn = planes + else: + check_len(dir, 3) + rtn = dirs + + if crystal_symm == 'hexagonal': + # L matrix for transforming directions + l_matrix = create_l_matrix( + 1, 1, c_over_a, np.pi / 2, np.pi / 2, np.pi * 2 / 3 + ) + # Q matrix for transforming planes + q_matrix = create_q_matrix(l_matrix) + + if dir is None: + plane = np.matmul(q_matrix, plane) + else: + dir = np.matmul(l_matrix, dir) + + for i, symm in enumerate(symmetries): + if dir is None: + plane_symm = symm.transform_vector(plane) + if plane_symm[2] < 0: + plane_symm *= -1 + if crystal_symm == 'hexagonal': + # q_matrix inverse is equal to l_matrix transposed and vice-versa + plane_symm = reduce_idc(convert_idc( + 'm', plane=safe_int_cast(np.matmul(l_matrix.T, plane_symm)) + )) + planes.append(safe_int_cast(plane_symm)) + else: + dir_symm = symm.transform_vector(dir) + if dir_symm[2] < 0: + dir_symm *= -1 + if crystal_symm == 'hexagonal': + dir_symm = reduce_idc(convert_idc( + 'm', dir=safe_int_cast(np.matmul(q_matrix.T, dir_symm)) + )) + dirs.append(safe_int_cast(dir_symm)) + + return rtn + + +def project_to_orth(c_over_a, *, dir=None, plane=None, in_type='mb'): + """ + Project from crystal aligned coordinates to an orthogonal set. + + Parameters + ---------- + in_type : str {'m', 'mb'} + Type of indices provided + dir : tuple of int or equivalent, optional + Direction to convert. This OR `plane` must me provided. + plane : tuple of int or equivalent, optional + Plane to convert. This OR `direction` must me provided. + + Returns + ------- + + + """ + if dir is None and plane is None: + raise ValueError("One of either `direction` or `plane` must be " + "provided.") + if dir is not None and plane is not None: + raise ValueError("One of either `direction` or `plane` must be " + "provided, not both.") + + if in_type == 'mb': + if dir is None: + check_len(plane, 4) + plane = convert_idc('mb', plane=plane) + if plane is None: + check_len(dir, 4) + dir = convert_idc('mb', dir=dir) + elif in_type != 'm': + raise ValueError("`inType` must be either 'm' or 'mb'.") + + # L matrix for transforming directions + l_matrix = create_l_matrix( + 1, 1, c_over_a, np.pi / 2, np.pi / 2, np.pi * 2 / 3 + ) + + if dir is None: + check_len(plane, 3) + # Q matrix for transforming planes + q_matrix = create_q_matrix(l_matrix) + return np.matmul(q_matrix, plane) + else: + check_len(dir, 3) + return np.matmul(l_matrix, dir) + + +def pos_idc(vec): + """ + Return a consistent positive version of a set of indices. + + Parameters + ---------- + vec : tuple of int or equivalent + Indices to convert. + + Returns + ------- + tuple of int + Positive version of indices. + + """ + for idx in vec: + if idx == 0: + continue + if idx > 0: + return tuple(vec) + else: + return tuple(-np.array(vec)) + + +def reduce_idc(vec): + """ + Reduce indices to lowest integers + + Parameters + ---------- + vec : tuple of int or equivalent + Indices to reduce. + + Returns + ------- + tuple of int + The reduced indices. + + """ + return tuple((np.array(vec) / np.gcd.reduce(vec)).astype(np.int8)) + + +def safe_int_cast(vec, tol=1e-3): + """ + Cast a tuple of floats to integers, raising an error if rounding is + over a tolerance. + + Parameters + ---------- + vec : tuple of float or equivalent + Vector to cast. + tol : float + Tolerance above which an error is raised. + + Returns + ------- + tuple of int + + Raises + ------ + ValueError + If the rounding is over the tolerance for any value. + + """ + vec = np.array(vec) + vec_rounded = vec.round() + + if np.any(np.abs(vec - vec_rounded) > tol): + raise ValueError('Rounding too large', np.abs(vec - vec_rounded)) + + return tuple(vec_rounded.astype(np.int8)) + + +def idc_to_string(idc, brackets=None, str_type='unicode'): + """ + String representation of a set of indicies. + + Parameters + ---------- + idc : collection of int + brackets : str + String of opening and closing brackets eg '()' + str_type : str {'unicode', 'tex'} + + Returns + ------- + str + + """ + text = ''.join(map(partial(str_idx, str_type=str_type), idc)) + if brackets is not None: + text = brackets[0] + text + brackets[1] + return text + + +def str_idx(idx, str_type='unicode'): + """ + String representation of an index with overbars. + + Parameters + ---------- + idx : int + str_type : str {'unicode', 'tex'} + + Returns + ------- + str + + """ + if not isinstance(idx, (int, np.integer)): + raise ValueError("Index must be an integer.") + + pre, post = (r'$\bar{', r'}$') if str_type == 'tex' else ('', u'\u0305') + return str(idx) if idx >= 0 else pre + str(-idx) + post \ No newline at end of file diff --git a/defdap/plotting.py b/defdap/plotting.py index 0f688be..8c7f73e 100644 --- a/defdap/plotting.py +++ b/defdap/plotting.py @@ -13,6 +13,8 @@ # See the License for the specific language governing permissions and # limitations under the License. +from functools import partial + import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt @@ -28,6 +30,7 @@ from defdap import defaults from defdap import quat +from defdap.crystal_utils import project_to_orth, equavlent_indicies, idc_to_string # TODO: add plot parameter to add to current figure @@ -993,67 +996,51 @@ def add_axis(self): """ if self.plot_type == "IPF" and self.crystal_sym == "cubic": - # line between [001] and [111] - self.add_line([0, 0, 1], [1, 1, 1], c='k', lw=2) - - # line between [001] and [101] - self.add_line([0, 0, 1], [1, 0, 1], c='k', lw=2) - - # line between [101] and [111] - self.add_line([1, 0, 1], [1, 1, 1], c='k', lw=2) - - # label poles - self.label_point([0, 0, 1], '001', - pad_y=-0.005, va='top', ha='center', fontsize=12) - self.label_point([1, 0, 1], '101', - pad_y=-0.005, va='top', ha='center', fontsize=12) - self.label_point([1, 1, 1], '111', - pad_y=0.005, va='bottom', ha='center', fontsize=12) + lines = [ + ((0, 0, 1), (0, 1, 1)), + ((0, 0, 1), (-1, 1, 1)), + ((0, 1, 1), (-1, 1, 1)), + ] + labels = [ + ((0, 0, 1), -0.005, 'top'), + ((0, 1, 1), -0.005, 'top'), + ((-1, 1, 1), 0.005, 'bottom'), + ] elif self.plot_type == "IPF" and self.crystal_sym == "hexagonal": - - triangle = defaults['ipf_triangle_convention'] - - if triangle == 'aztec': - # line between [0001] and [01-10] ([001] and [2-10]) - # converted to cubic axes - self.add_line([0, 0, 1], [np.sqrt(3), -1, 0], c='k', lw=2) - - # line between [0001] and [-12-10] ([001] and [100]) - self.add_line([0, 0, 1], [1, 0, 0], c='k', lw=2) - - # line between [-12-10] and [01-10] ([100] and [2-10]) - self.add_line([1, 0, 0], [np.sqrt(3), -1, 0], c='k', lw=2) - - # label poles - self.label_point([0, 0, 1], '0001', - pad_y=0.012, va='bottom', ha='center', fontsize=12) - self.label_point([1, 0, 0], r'$\bar{1}2\bar{1}0$', - pad_y=0.012, va='bottom', ha='center', fontsize=12) - self.label_point([np.sqrt(3), -1, 0], r'$01\bar{1}0$', - pad_y=-0.006, va='top', ha='center', fontsize=12) - - elif triangle == 'mtex': - # line between [0001] and [10-10] ([001] and [210]) - # converted to cubic axes - self.add_line([0, 0, 1], [np.sqrt(3), 1, 0], c='k', lw=2) - - # line between [0001] and [2-1-10] ([001] and [100]) - self.add_line([0, 0, 1], [1, 0, 0], c='k', lw=2) - - # line between [2-1-10] and [10-10] ([100] and [210]) - self.add_line([1, 0, 0], [np.sqrt(3), 1, 0], c='k', lw=2) - - # label poles - self.label_point([0, 0, 1], '0001', - pad_y=-0.012, va='top', ha='center', fontsize=12) - self.label_point([1, 0, 0], r'$\bar{1}2\bar{1}0$', - pad_y=-0.012, va='top', ha='center', fontsize=12) - self.label_point([np.sqrt(3), 1, 0], r'$\bar{1}100$', - pad_y=0.009, va='bottom', ha='center', fontsize=12) - + if defaults['ipf_triangle_convention'] == 'down': + lines = [ + ((0, 0, 0, 1), (-1, 2, -1, 0)), + ((0, 0, 0, 1), (0, 1, -1, 0)), + ((-1, 2, -1, 0), (0, 1, -1, 0)), + ] + labels = [ + ((0, 0, 0, 1), 0.012, 'bottom'), + ((-1, 2, -1, 0), 0.012, 'bottom'), + ((0, 1, -1, 0), -0.006, 'top'), + ] + + else: + lines = [ + ((0, 0, 0, 1), (-1, 2, -1, 0)), + ((0, 0, 0, 1), (-1, 1, 0, 0)), + ((-1, 2, -1, 0), (-1, 1, 0, 0)), + ] + labels = [ + ((0, 0, 0, 1), -0.006, 'top'), + ((-1, 2, -1, 0), -0.006, 'top'), + ((-1, 1, 0, 0), 0.012, 'bottom'), + ] else: raise NotImplementedError("Only works for cubic and hexagonal.") + + for line in lines: + self.add_line(*line, c='k', lw=2) + for label in labels: + self.label_point( + label[0], pad_y=label[1], va=label[2], + ha='center', fontsize=12 + ) self.ax.axis('equal') self.ax.axis('off') @@ -1063,9 +1050,9 @@ def add_line(self, start_point, end_point, plot_syms=False, res=100, **kwargs): Parameters ---------- - start_point : numpy.ndarray + start_point : tuple Start point in crystal coordinates (i.e. [0,0,1]). - end_point : numpy.ndarray + end_point : tuple End point in crystal coordinates, (i.e. [1,0,0]). plot_syms : bool, optional If true, plot all symmetrically equivelant points. @@ -1075,11 +1062,15 @@ def add_line(self, start_point, end_point, plot_syms=False, res=100, **kwargs): All other arguments are passed to :func:`matplotlib.pyplot.plot`. """ + if self.crystal_sym == 'hexagonal': + start_point = project_to_orth(0.8165, dir=start_point, in_type='mb') + end_point = project_to_orth(0.8165, dir=end_point, in_type='mb') + lines = [(start_point, end_point)] if plot_syms: for symm in quat.Quat.sym_eqv(self.crystal_sym)[1:]: - start_point_symm = symm.transform_vector(start_point).astype(int) - end_point_symm = symm.transform_vector(end_point).astype(int) + start_point_symm = symm.transform_vector(start_point) + end_point_symm = symm.transform_vector(end_point) if start_point_symm[2] < 0: start_point_symm *= -1 @@ -1099,14 +1090,14 @@ def add_line(self, start_point, end_point, plot_syms=False, res=100, **kwargs): xp, yp = self.projection(line_points[0], line_points[1], line_points[2]) self.ax.plot(xp, yp, **kwargs) - def label_point(self, point, label, pad_x=0, pad_y=0, **kwargs): + def label_point(self, point, label=None, plot_syms=False, pad_x=0, pad_y=0, **kwargs): """Place a label near a coordinate in the pole plot. Parameters ---------- point : tuple (x, y) coordinate to place text. - label : str + label : str, optional Text to use in label. pad_x : int, optional Pad added to x coordinate. @@ -1116,8 +1107,37 @@ def label_point(self, point, label, pad_x=0, pad_y=0, **kwargs): Other arguments are passed to :func:`matplotlib.axes.Axes.text`. """ - xp, yp = self.projection(*point) - self.ax.text(xp + pad_x, yp + pad_y, label, **kwargs) + labels = [idc_to_string(point, str_type='tex')] if label is None else [label] + + point_idc = point + if self.crystal_sym == 'hexagonal': + point = project_to_orth(0.8165, dir=point, in_type='mb') + + points = [point] + + if plot_syms: + for symm in quat.Quat.sym_eqv(self.crystal_sym)[1:]: + point_symm = symm.transform_vector(point) + if point_symm[2] < 0: + point_symm *= -1 + points.append(point_symm) + + if label is None: + labels = map( + partial(idc_to_string, str_type='tex'), + equavlent_indicies( + self.crystal_sym, + quat.Quat.sym_eqv(self.crystal_sym), + dir=point_idc, + c_over_a=0.8165 + ) + ) + else: + labels *= len(quat.Quat.sym_eqv(self.crystal_sym)) + + for point, label in zip(points, labels): + xp, yp = self.projection(*point) + self.ax.text(xp + pad_x, yp + pad_y, label, **kwargs) def add_points(self, alpha_ang, beta_ang, marker_colour=None, marker_size=None, **kwargs): """Add a point to the pole plot. @@ -1196,7 +1216,14 @@ def add_colour_bar(self, label, layer=0, **kwargs): img = self.img_layers[layer] self.colour_bar = plt.colorbar(img, ax=self.ax, label=label, **kwargs) - def add_legend(self, label='Grain area (μm$^2$)', number=6, layer=0, scaling=1, **kwargs): + def add_legend( + self, + label='Grain area (μm$^2$)', + number=6, + layer=0, + scaling=1, + **kwargs + ): """Add a marker size legend to the pole plot. Parameters @@ -1214,8 +1241,11 @@ def add_legend(self, label='Grain area (μm$^2$)', number=6, layer=0, scaling=1, """ img = self.img_layers[layer] - self.legend = plt.legend(*img.legend_elements("sizes", num=number, - func=lambda s: s / scaling), title=label, **kwargs) + self.legend = plt.legend( + *img.legend_elements("sizes", num=number, func=lambda s: s / scaling), + title=label, + **kwargs + ) @staticmethod def _validateProjection(projection_in, validate_default=False): @@ -1279,8 +1309,8 @@ def stereo_project(*args): raise Exception("3 arguments for pole directions and 2 for polar angles.") alpha_comp = np.tan(alpha / 2) - xp = alpha_comp * np.cos(beta) - yp = alpha_comp * np.sin(beta) + xp = alpha_comp * np.cos(beta - np.pi/2) + yp = alpha_comp * np.sin(beta - np.pi/2) return xp, yp @@ -1312,8 +1342,8 @@ def lambert_project(*args): raise Exception("3 arguments for pole directions and 2 for polar angles.") alpha_comp = np.sqrt(2 * (1 - np.cos(alpha))) - xp = alpha_comp * np.cos(beta) - yp = alpha_comp * np.sin(beta) + xp = alpha_comp * np.cos(beta - np.pi/2) + yp = alpha_comp * np.sin(beta - np.pi/2) return xp, yp diff --git a/defdap/quat.py b/defdap/quat.py index aec971a..1b00b8d 100755 --- a/defdap/quat.py +++ b/defdap/quat.py @@ -495,12 +495,6 @@ def plot_ipf( alpha_fund, beta_fund = Quat.calc_fund_dirs(quats, direction, sym_group) - if defaults['ipf_triangle_convention'] == 'aztec': - beta_fund -= np.pi/6 - - if defaults['ipf_triangle_convention'] == 'mtex': - beta_fund = - beta_fund + np.pi/6 - if plot is None: plot = plotting.PolePlot( "IPF", sym_group, projection=projection, @@ -955,7 +949,8 @@ def calc_fund_dirs( quats: np.ndarray, direction: np.ndarray, sym_group: str, - dtype: Optional[type] = float + dtype: Optional[type] = float, + triangle: Optional[str] = None, ) -> Tuple[np.ndarray, np.ndarray]: """ Transform the sample direction to crystal coords based on the quats @@ -1041,45 +1036,41 @@ def calc_fund_dirs( # find the poles in the fundamental triangle if sym_group == "cubic": - # first beta should be between 0 and 45 deg leaving 3 - # symmetric equivalents per orientation - trial_poles = np.logical_and(beta >= 0, beta <= np.pi / 4) + beta_range = (np.pi / 2, 3/4 * np.pi, 3) + + elif sym_group == "hexagonal": + if triangle is None: + triangle = defaults['ipf_triangle_convention'] - # if less than 3 left need to expand search slightly to - # catch edge cases - if np.any(np.sum(trial_poles, axis=0) < 3): - delta_beta = 1e-8 - trial_poles = np.logical_and(beta >= -delta_beta, - beta <= np.pi / 4 + delta_beta) + if triangle == 'up': + beta_range = (np.pi / 2, 2/3 * np.pi, 1) + elif triangle == 'down': + beta_range = (1/3 * np.pi, np.pi / 2, 1) + else: + ValueError("`triangle` must be 'up' or 'down'") + else: + raise ValueError("sym_group must be cubic or hexagonal") + + trial_poles = np.logical_and(beta >= beta_range[0], beta <= beta_range[1]) + # expand search slightly to catch edge cases if needed + if np.any(np.sum(trial_poles, axis=0) < beta_range[2]): + delta_beta = 1e-8 + trial_poles = np.logical_and( + beta >= beta_range[0] - delta_beta, + beta <= beta_range[1] + delta_beta + ) + if sym_group == "cubic": # now of symmetric equivalents left we want the one with # minimum alpha - min_alpha_idx = np.nanargmin(np.where(trial_poles==False, np.nan, alpha), axis=0) - beta_fund = beta[min_alpha_idx, np.arange(len(min_alpha_idx))] - alpha_fund = alpha[min_alpha_idx, np.arange(len(min_alpha_idx))] - - elif sym_group == "hexagonal": - - # first beta should be between -30 and 0 deg leaving 1 - # symmetric equivalent per orientation - trial_poles = np.logical_and(beta >= 0, beta <= np.pi / 6) - # if less than 1 left need to expand search slightly to - # catch edge cases - if np.any(np.sum(trial_poles, axis=0) < 1): - delta_beta = 1e-8 - trial_poles = np.logical_and(beta >= -delta_beta, - beta <= np.pi / 6 + delta_beta) - + fund_idx = np.nanargmin(np.where(trial_poles, alpha, np.nan), axis=0) + else: # non-indexed points cause more than 1 symmetric equivalent, use this # to pick one and filter non-indexed points later - first_idx = (trial_poles==True).argmax(axis=0) - beta_fund = beta[first_idx, np.arange(len(first_idx))] - alpha_fund = alpha[first_idx, np.arange(len(first_idx))] - - else: - raise Exception("symGroup must be cubic or hexagonal") + fund_idx = trial_poles.argmax(axis=0) - return alpha_fund, beta_fund + fund_idx = (fund_idx, range(len(fund_idx))) + return alpha[fund_idx], beta[fund_idx] @staticmethod def sym_eqv(symGroup: str) -> List['Quat']: