From df1948d0a3e842166c600694f61b030ea3645574 Mon Sep 17 00:00:00 2001
From: Carter Francis <csfrancis@wisc.edu>
Date: Thu, 26 Oct 2023 10:47:35 -0500
Subject: [PATCH 01/11] New Feature: Add get_sample_reduced_fundamental

---
 orix/sampling/sample_generators.py | 58 ++++++++++++++++++++++++++++++
 1 file changed, 58 insertions(+)

diff --git a/orix/sampling/sample_generators.py b/orix/sampling/sample_generators.py
index aa2aae3b..e72b3feb 100644
--- a/orix/sampling/sample_generators.py
+++ b/orix/sampling/sample_generators.py
@@ -24,6 +24,11 @@
 from orix.quaternion.symmetry import get_point_group
 from orix.sampling.SO3_sampling import _three_uniform_samples_method, uniform_SO3_sample
 from orix.sampling._cubochoric_sampling import cubochoric_sampling
+from orix.quaternion import symmetry
+from orix.sampling import sample_S2
+from orix.vector import Vector3d
+
+from diffsims.rotations import ConstrainedRotation
 
 
 def get_sample_fundamental(
@@ -165,3 +170,56 @@ def _remove_larger_than_angle(rot: Rotation, max_angle: Union[int, float]) -> Ro
     mask = half_angles < half_angle
     rot_out = rot[mask]
     return rot_out
+
+
+def get_sample_reduced_fundamental(
+    resolution: float,
+    mesh: str = None,
+    point_group: Symmetry = None,
+) -> Rotation:
+    """Produces orientations to align various crystallographic directions with
+    the z-axis, with the constraint that the first Euler angle phi_1=0.
+    The crystallographic directions sample the fundamental zone, representing
+    the smallest region of symmetrically unique directions of the relevant
+    crystal system or point group.
+    Parameters
+    ----------
+    resolution
+        An angle in degrees representing the maximum angular distance to a
+        first nearest neighbor grid point.
+    mesh
+        Type of meshing of the sphere that defines how the grid is created. See
+        orix.sampling.sample_S2 for all the options. A suitable default is
+        chosen depending on the crystal system.
+        point_group
+        Symmetry operations that determines the unique directions. Defaults to
+        no symmetry, which means sampling all 3D unit vectors.
+    Returns
+    -------
+    ConstrainedRotation
+        (N, 3) array representing Euler angles for the different orientations
+    """
+    if point_group is None:
+        point_group = symmetry.C1
+
+    if mesh is None:
+        s2_auto_sampling_map = {
+            "triclinic": "icosahedral",
+            "monoclinic": "icosahedral",
+            "orthorhombic": "spherified_cube_edge",
+            "tetragonal": "spherified_cube_edge",
+            "cubic": "spherified_cube_edge",
+            "trigonal": "hexagonal",
+            "hexagonal": "hexagonal",
+        }
+        mesh = s2_auto_sampling_map[point_group.system]
+
+    s2_sample = sample_S2(resolution, method=mesh)
+    fundamental = s2_sample[s2_sample <= point_group.fundamental_sector]
+
+    phi = fundamental.polar
+    phi2 = (np.pi / 2 - fundamental.azimuth) % (2 * np.pi)
+    phi1 = np.zeros(phi2.shape[0])
+    euler_angles = np.vstack([phi1, phi, phi2]).T
+
+    return Rotation.from_euler(euler_angles, degrees=False)

From 50b4ad6899be3a438a1fc38d6a31e21b6b0ecfa2 Mon Sep 17 00:00:00 2001
From: Carter Francis <csfrancis@wisc.edu>
Date: Thu, 26 Oct 2023 11:04:31 -0500
Subject: [PATCH 02/11] Testing: Add Tests for get_sample_reduced_fundamental

---
 orix/sampling/__init__.py            |  3 +-
 orix/sampling/sample_generators.py   | 57 ++++++++++++++++++++++++++--
 orix/tests/sampling/test_sampling.py | 25 +++++++++++-
 3 files changed, 78 insertions(+), 7 deletions(-)

diff --git a/orix/sampling/__init__.py b/orix/sampling/__init__.py
index 400bbfa1..7c943b3e 100644
--- a/orix/sampling/__init__.py
+++ b/orix/sampling/__init__.py
@@ -29,10 +29,11 @@
 )
 from orix.sampling.S2_sampling import sampling_methods as sample_S2_methods
 from orix.sampling.SO3_sampling import uniform_SO3_sample
-from orix.sampling.sample_generators import get_sample_fundamental, get_sample_local
+from orix.sampling.sample_generators import get_sample_fundamental, get_sample_local, get_sample_reduced_fundamental
 
 __all__ = [
     "get_sample_fundamental",
+    "get_sample_reduced_fundamental",
     "get_sample_local",
     "sample_S2",
     "sample_S2_methods",
diff --git a/orix/sampling/sample_generators.py b/orix/sampling/sample_generators.py
index e72b3feb..8f250b81 100644
--- a/orix/sampling/sample_generators.py
+++ b/orix/sampling/sample_generators.py
@@ -26,9 +26,6 @@
 from orix.sampling._cubochoric_sampling import cubochoric_sampling
 from orix.quaternion import symmetry
 from orix.sampling import sample_S2
-from orix.vector import Vector3d
-
-from diffsims.rotations import ConstrainedRotation
 
 
 def get_sample_fundamental(
@@ -177,7 +174,7 @@ def get_sample_reduced_fundamental(
     mesh: str = None,
     point_group: Symmetry = None,
 ) -> Rotation:
-    """Produces orientations to align various crystallographic directions with
+    """Produces rotations to align various crystallographic directions with
     the z-axis, with the constraint that the first Euler angle phi_1=0.
     The crystallographic directions sample the fundamental zone, representing
     the smallest region of symmetrically unique directions of the relevant
@@ -223,3 +220,55 @@ def get_sample_reduced_fundamental(
     euler_angles = np.vstack([phi1, phi, phi2]).T
 
     return Rotation.from_euler(euler_angles, degrees=False)
+
+
+def get_sample_zone_axis(
+    resolution: float,
+    mesh: str = None,
+    point_group: Symmetry = None,
+) -> Rotation:
+    """Produces the rotations to align various crystallographic directions with
+    the z-axis, with the constraint that the first Euler angle phi_1=0.
+
+
+    Parameters
+    ----------
+    resolution
+        An angle in degrees representing the maximum angular distance to a
+        first nearest neighbor grid point.
+    mesh
+        Type of meshing of the sphere that defines how the grid is created. See
+        orix.sampling.sample_S2 for all the options. A suitable default is
+        chosen depending on the crystal system.
+        point_group
+        Symmetry operations that determines the unique directions. Defaults to
+        no symmetry, which means sampling all 3D unit vectors.
+    Returns
+    -------
+    ConstrainedRotation
+        (N, 3) array representing Euler angles for the different orientations
+    """
+    if point_group is None:
+        point_group = symmetry.C1
+
+    if mesh is None:
+        s2_auto_sampling_map = {
+            "triclinic": "icosahedral",
+            "monoclinic": "icosahedral",
+            "orthorhombic": "spherified_cube_edge",
+            "tetragonal": "spherified_cube_edge",
+            "cubic": "spherified_cube_edge",
+            "trigonal": "hexagonal",
+            "hexagonal": "hexagonal",
+        }
+        mesh = s2_auto_sampling_map[point_group.system]
+
+    s2_sample = sample_S2(resolution, method=mesh)
+    fundamental = s2_sample[s2_sample <= point_group.fundamental_sector]
+
+    phi = fundamental.polar
+    phi2 = (np.pi / 2 - fundamental.azimuth) % (2 * np.pi)
+    phi1 = np.zeros(phi2.shape[0])
+    euler_angles = np.vstack([phi1, phi, phi2]).T
+
+    return Rotation.from_euler(euler_angles, degrees=False)
diff --git a/orix/tests/sampling/test_sampling.py b/orix/tests/sampling/test_sampling.py
index 828f465d..ce6bfa33 100644
--- a/orix/tests/sampling/test_sampling.py
+++ b/orix/tests/sampling/test_sampling.py
@@ -20,8 +20,13 @@
 import pytest
 
 from orix.quaternion import Rotation
-from orix.quaternion.symmetry import C2, C6, D6, get_point_group
-from orix.sampling import get_sample_fundamental, get_sample_local, uniform_SO3_sample
+from orix.quaternion.symmetry import C1, C2, C4, C6, D6, get_point_group
+from orix.sampling import (
+    get_sample_fundamental,
+    get_sample_local,
+    uniform_SO3_sample,
+    get_sample_reduced_fundamental,
+)
 from orix.sampling.SO3_sampling import _resolution_to_num_steps
 from orix.sampling._polyhedral_sampling import (
     _get_angles_between_nn_gridpoints,
@@ -175,3 +180,19 @@ def test_get_sample_fundamental_space_group(self, C6_sample):
         C2_sample = get_sample_fundamental(4, space_group=3, method="haar_euler")
         ratio = C2_sample.size / C6_sample.size
         assert np.isclose(ratio, 3, atol=0.2)
+
+    def test_get_sample_reduced_fundamental(self):
+        rotations = get_sample_reduced_fundamental(resolution=4, point_group=C1)
+        rotations2 = get_sample_reduced_fundamental(resolution=4, point_group=C2)
+        rotations4 = get_sample_reduced_fundamental(resolution=4, point_group=C4)
+        rotations6 = get_sample_reduced_fundamental(resolution=4, point_group=C4)
+
+        assert (
+            np.abs(rotations.size / rotations2.size) - 2 < 0.1
+        )  # about 2 times more rotations
+        assert (
+            np.abs(rotations2.size / rotations4.size) - 2 < 0.1
+        )  # about 2 times more rotations
+        assert (
+            np.abs(rotations.size / rotations6.size) - 6 < 0.1
+        )  # about 6 times more rotations

From 8c0ca8c5e7da8b9b6d3b08042bb372dd0502b5e2 Mon Sep 17 00:00:00 2001
From: Carter Francis <csfrancis@wisc.edu>
Date: Thu, 26 Oct 2023 22:09:56 -0500
Subject: [PATCH 03/11] New Feature: Add sampling for ZoneAxis patterns

---
 examples/rotations/sampling_rotations.py |  37 +++++++
 orix/sampling/__init__.py                |   6 +-
 orix/sampling/sample_generators.py       | 130 +++++++++++++++--------
 orix/tests/sampling/test_sampling.py     |  25 +++++
 4 files changed, 153 insertions(+), 45 deletions(-)
 create mode 100644 examples/rotations/sampling_rotations.py

diff --git a/examples/rotations/sampling_rotations.py b/examples/rotations/sampling_rotations.py
new file mode 100644
index 00000000..fe9cac50
--- /dev/null
+++ b/examples/rotations/sampling_rotations.py
@@ -0,0 +1,37 @@
+"""
+==================
+Sampling rotations
+==================
+
+This example shows how to sample some phase object in Orix. We will
+get both the zone axis and the reduced fundamental zone rotations for
+the phase of interest.
+"""
+from diffpy.structure import Atom, Lattice, Structure
+from orix.crystal_map import Phase
+from orix.sampling import (get_sample_reduced_fundamental, get_sample_zone_axis,)
+from orix.vector import Vector3d
+
+a = 5.431
+latt = Lattice(a, a, a, 90, 90, 90)
+atom_list = []
+for coords in [[0, 0, 0], [0.5, 0, 0.5], [0, 0.5, 0.5], [0.5, 0.5, 0]]:
+    x, y, z = coords[0], coords[1], coords[2]
+    atom_list.append(Atom(atype="Si", xyz=[x, y, z], lattice=latt))  # Motif part A
+    atom_list.append(
+        Atom(atype="Si", xyz=[x + 0.25, y + 0.25, z + 0.25], lattice=latt)
+    )  # Motif part B
+struct = Structure(atoms=atom_list, lattice=latt)
+p = Phase(structure=struct, space_group=227)
+reduced_fun = get_sample_reduced_fundamental(resolution=4, point_group=p.point_group)
+
+vect_rot = reduced_fun*Vector3d.zvector()  # get the vector representation of the rotations
+vect_rot.scatter(grid=True)  # plot the stereographic projection of the rotations
+
+# %%
+
+zone_axis_rot, directions = get_sample_zone_axis(phase=p, density="7", return_directions=True)  # get the zone axis rotations
+zone_vect_rot = zone_axis_rot*Vector3d.zvector()  # get the vector representation of the rotations
+zone_vect_rot.scatter(grid=True, vector_labels=directions, text_kwargs={"size":8, "rotation":0})  # plot the stereographic projection of the rotations
+
+# %%
diff --git a/orix/sampling/__init__.py b/orix/sampling/__init__.py
index 7c943b3e..9ae6d27b 100644
--- a/orix/sampling/__init__.py
+++ b/orix/sampling/__init__.py
@@ -29,12 +29,16 @@
 )
 from orix.sampling.S2_sampling import sampling_methods as sample_S2_methods
 from orix.sampling.SO3_sampling import uniform_SO3_sample
-from orix.sampling.sample_generators import get_sample_fundamental, get_sample_local, get_sample_reduced_fundamental
+from orix.sampling.sample_generators import (get_sample_fundamental,
+                                             get_sample_local,
+                                             get_sample_reduced_fundamental,
+                                             get_sample_zone_axis)
 
 __all__ = [
     "get_sample_fundamental",
     "get_sample_reduced_fundamental",
     "get_sample_local",
+    "get_sample_zone_axis",
     "sample_S2",
     "sample_S2_methods",
     "uniform_SO3_sample",
diff --git a/orix/sampling/sample_generators.py b/orix/sampling/sample_generators.py
index 8f250b81..2258fda8 100644
--- a/orix/sampling/sample_generators.py
+++ b/orix/sampling/sample_generators.py
@@ -19,6 +19,7 @@
 from typing import Optional, Union
 
 import numpy as np
+from transforms3d.euler import axangle2euler
 
 from orix.quaternion import OrientationRegion, Rotation, Symmetry
 from orix.quaternion.symmetry import get_point_group
@@ -26,6 +27,7 @@
 from orix.sampling._cubochoric_sampling import cubochoric_sampling
 from orix.quaternion import symmetry
 from orix.sampling import sample_S2
+from orix.crystal_map import Phase
 
 
 def get_sample_fundamental(
@@ -193,7 +195,7 @@ def get_sample_reduced_fundamental(
         no symmetry, which means sampling all 3D unit vectors.
     Returns
     -------
-    ConstrainedRotation
+    Rotation
         (N, 3) array representing Euler angles for the different orientations
     """
     if point_group is None:
@@ -222,53 +224,93 @@ def get_sample_reduced_fundamental(
     return Rotation.from_euler(euler_angles, degrees=False)
 
 
-def get_sample_zone_axis(
-    resolution: float,
-    mesh: str = None,
-    point_group: Symmetry = None,
-) -> Rotation:
-    """Produces the rotations to align various crystallographic directions with
-    the z-axis, with the constraint that the first Euler angle phi_1=0.
-
-
+def _get_rotation_from_z_to_direction(structure, direction):
+    """
+    Finds the rotation that takes [001] to a given zone axis.
     Parameters
     ----------
-    resolution
-        An angle in degrees representing the maximum angular distance to a
-        first nearest neighbor grid point.
-    mesh
-        Type of meshing of the sphere that defines how the grid is created. See
-        orix.sampling.sample_S2 for all the options. A suitable default is
-        chosen depending on the crystal system.
-        point_group
-        Symmetry operations that determines the unique directions. Defaults to
-        no symmetry, which means sampling all 3D unit vectors.
+    structure : diffpy.structure.structure.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
     -------
-    ConstrainedRotation
-        (N, 3) array representing Euler angles for the different orientations
+    euler_angles : tuple
+        'rzxz' in degrees.
+    See Also
+    --------
+    generate_zap_map
+    :meth:`~diffsims.generators.rotation_list_generators.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
     """
-    if point_group is None:
-        point_group = symmetry.C1
-
-    if mesh is None:
-        s2_auto_sampling_map = {
-            "triclinic": "icosahedral",
-            "monoclinic": "icosahedral",
-            "orthorhombic": "spherified_cube_edge",
-            "tetragonal": "spherified_cube_edge",
-            "cubic": "spherified_cube_edge",
-            "trigonal": "hexagonal",
-            "hexagonal": "hexagonal",
-        }
-        mesh = s2_auto_sampling_map[point_group.system]
-
-    s2_sample = sample_S2(resolution, method=mesh)
-    fundamental = s2_sample[s2_sample <= point_group.fundamental_sector]
+    # 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 euler
+
+
+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]
 
-    phi = fundamental.polar
-    phi2 = (np.pi / 2 - fundamental.azimuth) % (2 * np.pi)
-    phi1 = np.zeros(phi2.shape[0])
-    euler_angles = np.vstack([phi1, phi, phi2]).T
 
-    return Rotation.from_euler(euler_angles, degrees=False)
+def get_sample_zone_axis(
+    density: str = "3",
+    phase: Phase = None,
+    return_directions: bool = False,
+) -> Rotation:
+    """Produces rotations to align various crystallographic directions with
+    the sample zone axes.
+    """
+    system = phase.point_group.system
+    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])
+    else:
+        raise ValueError("Density must be either 3 or 7")
+
+    euler_angles = np.array([
+        _get_rotation_from_z_to_direction(phase.structure, d) for d in direction_list
+    ])
+    if return_directions:
+        return Rotation.from_euler(euler_angles), direction_list
+    else:
+        return Rotation.from_euler(euler_angles)
diff --git a/orix/tests/sampling/test_sampling.py b/orix/tests/sampling/test_sampling.py
index ce6bfa33..e79dd4d6 100644
--- a/orix/tests/sampling/test_sampling.py
+++ b/orix/tests/sampling/test_sampling.py
@@ -19,13 +19,17 @@
 import numpy as np
 import pytest
 
+from diffpy.structure import Atom, Lattice, Structure
+
 from orix.quaternion import Rotation
+from orix.crystal_map import Phase
 from orix.quaternion.symmetry import C1, C2, C4, C6, D6, get_point_group
 from orix.sampling import (
     get_sample_fundamental,
     get_sample_local,
     uniform_SO3_sample,
     get_sample_reduced_fundamental,
+    get_sample_zone_axis,
 )
 from orix.sampling.SO3_sampling import _resolution_to_num_steps
 from orix.sampling._polyhedral_sampling import (
@@ -196,3 +200,24 @@ def test_get_sample_reduced_fundamental(self):
         assert (
             np.abs(rotations.size / rotations6.size) - 6 < 0.1
         )  # about 6 times more rotations
+
+    @pytest.mark.parametrize("density", ("3", "7"))
+    @pytest.mark.parametrize("get_directions", (True, False))
+    def test_get_zone_axis(self, density, get_directions):
+        a = 5.431
+        latt = Lattice(a, a, a, 90, 90, 90)
+        atom_list = []
+        for coords in [[0, 0, 0], [0.5, 0, 0.5], [0, 0.5, 0.5], [0.5, 0.5, 0]]:
+            x, y, z = coords[0], coords[1], coords[2]
+            atom_list.append(Atom(atype="Si", xyz=[x, y, z], lattice=latt))  # Motif part A
+            atom_list.append(
+                Atom(atype="Si", xyz=[x + 0.25, y + 0.25, z + 0.25], lattice=latt)
+            )  # Motif part B
+        struct = Structure(atoms=atom_list, lattice=latt)
+        p = Phase(structure=struct, space_group=227)
+        if get_directions:
+            rot, _ = get_sample_zone_axis(phase=p, density=density, return_directions=True)
+        else:
+            rot = get_sample_zone_axis(phase=p, density=density)
+        assert isinstance(rot, Rotation)
+

From 406704bb2cd3663cbb64440492caedca5fdaa3a2 Mon Sep 17 00:00:00 2001
From: Carter Francis <csfrancis@wisc.edu>
Date: Fri, 27 Oct 2023 08:44:51 -0500
Subject: [PATCH 04/11] Refactor: Removed transform3d dependency

---
 examples/rotations/sampling_rotations.py | 21 ++++++++++++++++-----
 orix/sampling/__init__.py                | 10 ++++++----
 orix/sampling/sample_generators.py       | 17 +++++++++++------
 orix/tests/sampling/test_sampling.py     |  9 ++++++---
 4 files changed, 39 insertions(+), 18 deletions(-)

diff --git a/examples/rotations/sampling_rotations.py b/examples/rotations/sampling_rotations.py
index fe9cac50..630d8856 100644
--- a/examples/rotations/sampling_rotations.py
+++ b/examples/rotations/sampling_rotations.py
@@ -9,7 +9,10 @@
 """
 from diffpy.structure import Atom, Lattice, Structure
 from orix.crystal_map import Phase
-from orix.sampling import (get_sample_reduced_fundamental, get_sample_zone_axis,)
+from orix.sampling import (
+    get_sample_reduced_fundamental,
+    get_sample_zone_axis,
+)
 from orix.vector import Vector3d
 
 a = 5.431
@@ -25,13 +28,21 @@
 p = Phase(structure=struct, space_group=227)
 reduced_fun = get_sample_reduced_fundamental(resolution=4, point_group=p.point_group)
 
-vect_rot = reduced_fun*Vector3d.zvector()  # get the vector representation of the rotations
+vect_rot = (
+    reduced_fun * Vector3d.zvector()
+)  # get the vector representation of the rotations
 vect_rot.scatter(grid=True)  # plot the stereographic projection of the rotations
 
 # %%
 
-zone_axis_rot, directions = get_sample_zone_axis(phase=p, density="7", return_directions=True)  # get the zone axis rotations
-zone_vect_rot = zone_axis_rot*Vector3d.zvector()  # get the vector representation of the rotations
-zone_vect_rot.scatter(grid=True, vector_labels=directions, text_kwargs={"size":8, "rotation":0})  # plot the stereographic projection of the rotations
+zone_axis_rot, directions = get_sample_zone_axis(
+    phase=p, density="7", return_directions=True
+)  # get the zone axis rotations
+zone_vect_rot = (
+    zone_axis_rot * Vector3d.zvector()
+)  # get the vector representation of the rotations
+zone_vect_rot.scatter(
+    grid=True, vector_labels=directions, text_kwargs={"size": 8, "rotation": 0}
+)  # plot the stereographic projection of the rotations
 
 # %%
diff --git a/orix/sampling/__init__.py b/orix/sampling/__init__.py
index 9ae6d27b..e0a94b16 100644
--- a/orix/sampling/__init__.py
+++ b/orix/sampling/__init__.py
@@ -29,10 +29,12 @@
 )
 from orix.sampling.S2_sampling import sampling_methods as sample_S2_methods
 from orix.sampling.SO3_sampling import uniform_SO3_sample
-from orix.sampling.sample_generators import (get_sample_fundamental,
-                                             get_sample_local,
-                                             get_sample_reduced_fundamental,
-                                             get_sample_zone_axis)
+from orix.sampling.sample_generators import (
+    get_sample_fundamental,
+    get_sample_local,
+    get_sample_reduced_fundamental,
+    get_sample_zone_axis,
+)
 
 __all__ = [
     "get_sample_fundamental",
diff --git a/orix/sampling/sample_generators.py b/orix/sampling/sample_generators.py
index 2258fda8..e96bd77c 100644
--- a/orix/sampling/sample_generators.py
+++ b/orix/sampling/sample_generators.py
@@ -19,9 +19,9 @@
 from typing import Optional, Union
 
 import numpy as np
-from transforms3d.euler import axangle2euler
 
 from orix.quaternion import OrientationRegion, Rotation, Symmetry
+from orix.vector import Vector3d
 from orix.quaternion.symmetry import get_point_group
 from orix.sampling.SO3_sampling import _three_uniform_samples_method, uniform_SO3_sample
 from orix.sampling._cubochoric_sampling import cubochoric_sampling
@@ -307,10 +307,15 @@ def get_sample_zone_axis(
     else:
         raise ValueError("Density must be either 3 or 7")
 
-    euler_angles = np.array([
-        _get_rotation_from_z_to_direction(phase.structure, d) for d in direction_list
-    ])
+    # rotate the directions to the z axis
+    rots = np.stack(
+        [
+            Rotation.from_align_vectors(v, Vector3d.zvector()).data
+            for v in direction_list
+        ]
+    )
+    rotations = Rotation(rots)
     if return_directions:
-        return Rotation.from_euler(euler_angles), direction_list
+        return rotations, direction_list
     else:
-        return Rotation.from_euler(euler_angles)
+        return rotations
diff --git a/orix/tests/sampling/test_sampling.py b/orix/tests/sampling/test_sampling.py
index e79dd4d6..97fd3875 100644
--- a/orix/tests/sampling/test_sampling.py
+++ b/orix/tests/sampling/test_sampling.py
@@ -209,15 +209,18 @@ def test_get_zone_axis(self, density, get_directions):
         atom_list = []
         for coords in [[0, 0, 0], [0.5, 0, 0.5], [0, 0.5, 0.5], [0.5, 0.5, 0]]:
             x, y, z = coords[0], coords[1], coords[2]
-            atom_list.append(Atom(atype="Si", xyz=[x, y, z], lattice=latt))  # Motif part A
+            atom_list.append(
+                Atom(atype="Si", xyz=[x, y, z], lattice=latt)
+            )  # Motif part A
             atom_list.append(
                 Atom(atype="Si", xyz=[x + 0.25, y + 0.25, z + 0.25], lattice=latt)
             )  # Motif part B
         struct = Structure(atoms=atom_list, lattice=latt)
         p = Phase(structure=struct, space_group=227)
         if get_directions:
-            rot, _ = get_sample_zone_axis(phase=p, density=density, return_directions=True)
+            rot, _ = get_sample_zone_axis(
+                phase=p, density=density, return_directions=True
+            )
         else:
             rot = get_sample_zone_axis(phase=p, density=density)
         assert isinstance(rot, Rotation)
-

From 7c8a04ccd7f42f365cdfff0c0a0cacdd7ab58f25 Mon Sep 17 00:00:00 2001
From: Carter Francis <csfrancis@wisc.edu>
Date: Fri, 27 Oct 2023 08:47:53 -0500
Subject: [PATCH 05/11] Refactor: Cleaned up doc strings and unused code

---
 orix/sampling/sample_generators.py   | 48 ++++++----------------------
 orix/tests/sampling/test_sampling.py | 25 ++++++++-------
 2 files changed, 24 insertions(+), 49 deletions(-)

diff --git a/orix/sampling/sample_generators.py b/orix/sampling/sample_generators.py
index e96bd77c..0e986524 100644
--- a/orix/sampling/sample_generators.py
+++ b/orix/sampling/sample_generators.py
@@ -190,7 +190,7 @@ def get_sample_reduced_fundamental(
         Type of meshing of the sphere that defines how the grid is created. See
         orix.sampling.sample_S2 for all the options. A suitable default is
         chosen depending on the crystal system.
-        point_group
+    point_group
         Symmetry operations that determines the unique directions. Defaults to
         no symmetry, which means sampling all 3D unit vectors.
     Returns
@@ -224,43 +224,6 @@ def get_sample_reduced_fundamental(
     return Rotation.from_euler(euler_angles, degrees=False)
 
 
-def _get_rotation_from_z_to_direction(structure, direction):
-    """
-    Finds the rotation that takes [001] to a given zone axis.
-    Parameters
-    ----------
-    structure : diffpy.structure.structure.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
-    :meth:`~diffsims.generators.rotation_list_generators.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 euler
-
-
 def _corners_to_centroid_and_edge_centers(corners):
     """
     Produces the midpoints and center of a trio of corners
@@ -290,6 +253,15 @@ def get_sample_zone_axis(
 ) -> Rotation:
     """Produces rotations to align various crystallographic directions with
     the sample zone axes.
+
+    Parameters
+    ----------
+    density
+        Either '3' or '7' for the number of directions to return.
+    phase
+        The phase for which the zone axis rotations are required.
+    return_directions
+        If True, returns the directions as well as the rotations.
     """
     system = phase.point_group.system
     corners_dict = {
diff --git a/orix/tests/sampling/test_sampling.py b/orix/tests/sampling/test_sampling.py
index 97fd3875..7e7a9072 100644
--- a/orix/tests/sampling/test_sampling.py
+++ b/orix/tests/sampling/test_sampling.py
@@ -16,20 +16,19 @@
 # You should have received a copy of the GNU General Public License
 # along with orix.  If not, see <http://www.gnu.org/licenses/>.
 
+from diffpy.structure import Atom, Lattice, Structure
 import numpy as np
 import pytest
 
-from diffpy.structure import Atom, Lattice, Structure
-
-from orix.quaternion import Rotation
 from orix.crystal_map import Phase
+from orix.quaternion import Rotation
 from orix.quaternion.symmetry import C1, C2, C4, C6, D6, get_point_group
 from orix.sampling import (
     get_sample_fundamental,
     get_sample_local,
-    uniform_SO3_sample,
     get_sample_reduced_fundamental,
     get_sample_zone_axis,
+    uniform_SO3_sample,
 )
 from orix.sampling.SO3_sampling import _resolution_to_num_steps
 from orix.sampling._polyhedral_sampling import (
@@ -201,7 +200,7 @@ def test_get_sample_reduced_fundamental(self):
             np.abs(rotations.size / rotations6.size) - 6 < 0.1
         )  # about 6 times more rotations
 
-    @pytest.mark.parametrize("density", ("3", "7"))
+    @pytest.mark.parametrize("density", ("3", "7", "5"))
     @pytest.mark.parametrize("get_directions", (True, False))
     def test_get_zone_axis(self, density, get_directions):
         a = 5.431
@@ -217,10 +216,14 @@ def test_get_zone_axis(self, density, get_directions):
             )  # Motif part B
         struct = Structure(atoms=atom_list, lattice=latt)
         p = Phase(structure=struct, space_group=227)
-        if get_directions:
-            rot, _ = get_sample_zone_axis(
-                phase=p, density=density, return_directions=True
-            )
+        if density == "5":
+            with pytest.raises(ValueError):
+                get_sample_zone_axis(phase=p, density=density)
         else:
-            rot = get_sample_zone_axis(phase=p, density=density)
-        assert isinstance(rot, Rotation)
+            if get_directions:
+                rot, _ = get_sample_zone_axis(
+                    phase=p, density=density, return_directions=True
+                )
+            else:
+                rot = get_sample_zone_axis(phase=p, density=density)
+            assert isinstance(rot, Rotation)

From d0fa073ab14238274d2397ace7ba3503fc3e4268 Mon Sep 17 00:00:00 2001
From: Carter Francis <csfrancis@wisc.edu>
Date: Thu, 2 Nov 2023 14:27:52 -0400
Subject: [PATCH 06/11] BugFix: Fixed Quotes around error

---
 examples/rotations/sampling_rotations.py | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/examples/rotations/sampling_rotations.py b/examples/rotations/sampling_rotations.py
index 630d8856..9d59d6c8 100644
--- a/examples/rotations/sampling_rotations.py
+++ b/examples/rotations/sampling_rotations.py
@@ -8,11 +8,9 @@
 the phase of interest.
 """
 from diffpy.structure import Atom, Lattice, Structure
+
 from orix.crystal_map import Phase
-from orix.sampling import (
-    get_sample_reduced_fundamental,
-    get_sample_zone_axis,
-)
+from orix.sampling import get_sample_reduced_fundamental, get_sample_zone_axis
 from orix.vector import Vector3d
 
 a = 5.431

From 94d3df883714d733c74de14fc58fdab375ed68ae Mon Sep 17 00:00:00 2001
From: Carter Francis <csfrancis@wisc.edu>
Date: Mon, 6 Nov 2023 12:13:11 -0600
Subject: [PATCH 07/11] Formatting: Run isort

---
 orix/sampling/sample_generators.py | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/orix/sampling/sample_generators.py b/orix/sampling/sample_generators.py
index 0e986524..7600a4fe 100644
--- a/orix/sampling/sample_generators.py
+++ b/orix/sampling/sample_generators.py
@@ -20,14 +20,13 @@
 
 import numpy as np
 
-from orix.quaternion import OrientationRegion, Rotation, Symmetry
-from orix.vector import Vector3d
+from orix.crystal_map import Phase
+from orix.quaternion import OrientationRegion, Rotation, Symmetry, symmetry
 from orix.quaternion.symmetry import get_point_group
+from orix.sampling import sample_S2
 from orix.sampling.SO3_sampling import _three_uniform_samples_method, uniform_SO3_sample
 from orix.sampling._cubochoric_sampling import cubochoric_sampling
-from orix.quaternion import symmetry
-from orix.sampling import sample_S2
-from orix.crystal_map import Phase
+from orix.vector import Vector3d
 
 
 def get_sample_fundamental(

From 75974fe791d238e3c8c34719e29ecf490feb2869 Mon Sep 17 00:00:00 2001
From: Carter Francis <csfrancis@wisc.edu>
Date: Tue, 26 Mar 2024 11:22:04 -0500
Subject: [PATCH 08/11] Documentation: Add orcid ID

---
 .zenodo.json | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/.zenodo.json b/.zenodo.json
index f85c7be3..60b9cc92 100644
--- a/.zenodo.json
+++ b/.zenodo.json
@@ -34,6 +34,11 @@
       "name": "Anders Christian Mathisen",
       "affiliation": "Norwegian University of Science and Technology"
     },
+    {
+      "name": "Carter Francis",
+      "orcid": "0000-0003-2564-1851",
+      "affiliation": "University of Wisconsin Madison"
+    },
     {
       "name": "Simon Høgås"
     },

From 0700b5d5794cb52db6d91f8ed44aadbc930e44b0 Mon Sep 17 00:00:00 2001
From: Carter Francis <csfrancis@wisc.edu>
Date: Tue, 26 Mar 2024 11:53:09 -0500
Subject: [PATCH 09/11] Coverage: Update coverage

---
 orix/tests/sampling/test_sampling.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/orix/tests/sampling/test_sampling.py b/orix/tests/sampling/test_sampling.py
index 7e7a9072..6517ddfa 100644
--- a/orix/tests/sampling/test_sampling.py
+++ b/orix/tests/sampling/test_sampling.py
@@ -185,7 +185,7 @@ def test_get_sample_fundamental_space_group(self, C6_sample):
         assert np.isclose(ratio, 3, atol=0.2)
 
     def test_get_sample_reduced_fundamental(self):
-        rotations = get_sample_reduced_fundamental(resolution=4, point_group=C1)
+        rotations = get_sample_reduced_fundamental(resolution=4)
         rotations2 = get_sample_reduced_fundamental(resolution=4, point_group=C2)
         rotations4 = get_sample_reduced_fundamental(resolution=4, point_group=C4)
         rotations6 = get_sample_reduced_fundamental(resolution=4, point_group=C4)

From 00f6bf169851c9cc6fb6cb1ff5c97efa62cfbd6e Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Tue, 26 Mar 2024 17:26:37 +0000
Subject: [PATCH 10/11] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 examples/rotations/sampling_rotations.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/examples/rotations/sampling_rotations.py b/examples/rotations/sampling_rotations.py
index 9d59d6c8..cc3c32a6 100644
--- a/examples/rotations/sampling_rotations.py
+++ b/examples/rotations/sampling_rotations.py
@@ -7,6 +7,7 @@
 get both the zone axis and the reduced fundamental zone rotations for
 the phase of interest.
 """
+
 from diffpy.structure import Atom, Lattice, Structure
 
 from orix.crystal_map import Phase

From b45df097c6f3af85d0abaf604f39980dac6c9699 Mon Sep 17 00:00:00 2001
From: Carter Francis <csfrancis@wisc.edu>
Date: Fri, 29 Mar 2024 12:04:57 -0500
Subject: [PATCH 11/11] Documentation: Update CHANGELOG.rst

---
 CHANGELOG.rst | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 947aff41..b4b8b878 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -31,6 +31,9 @@ Added
 - The ``random()`` methods of ``Orientation`` and ``Misorientation`` now accept
   ``symmetry``. A ``random()`` method is also added to ``Vector3d`` and ``Miller``, the
   latter accepting a ``phase``.
+- ``Added orix.sampling.get_sample_zone_axis`` for getting zone axes for some point group.
+- ``Added orix.sampling.get_sample_reduced_fundamental`` for getting reduced
+  fundamental zone for some point group.
 
 Changed
 -------