From 5f8d785ab41be0a3b2e3717d09c9a53221b43e30 Mon Sep 17 00:00:00 2001
From: Andrew Kerr <andrewghkerr@gmail.com>
Date: Tue, 15 Jun 2021 18:14:11 -0400
Subject: [PATCH 01/10] Additional explanation and example in diffraction docs.

---
 freud/diffraction.pyx | 16 +++++++++++++++-
 1 file changed, 15 insertions(+), 1 deletion(-)

diff --git a/freud/diffraction.pyx b/freud/diffraction.pyx
index d48f354df..b014896d9 100644
--- a/freud/diffraction.pyx
+++ b/freud/diffraction.pyx
@@ -217,7 +217,21 @@ cdef class DiffractionPattern(_Compute):
                 Whether to erase the previously computed values before adding
                 the new computations; if False, will accumulate data (Default
                 value: True).
-        """
+
+            For a target view axis in Cartesian coordinates,
+            `rowan <https://rowan.readthedocs.io/en/latest/index.html>__` provides
+            ``vector_vector_rotation()``, a function that will return the necessary
+            quaternion for rotating the default view orientation to the desired axis.
+
+            Example::
+
+                >>> import rowan
+                >>> import numpy as np
+                >>> default_view_axis = np.array([0, 0, 1])
+                >>> desired_view_axis = np.array([1, 1, 1])
+                >>> quaternion = rowan.vector_vector_rotation(desired_view_axis, default_view_axis) # use for view_orientation
+
+        """ # noqa: E501
         if reset:
             self._diffraction = np.zeros((self.output_size, self.output_size))
             self._frame_counter = 0

From 805aadf886a1f910d62c5fd5f7f97b8c9a0fb393 Mon Sep 17 00:00:00 2001
From: Andrew Kerr <andrewghkerr@gmail.com>
Date: Tue, 15 Jun 2021 18:56:37 -0400
Subject: [PATCH 02/10] Docs edit.

---
 freud/diffraction.pyx | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/freud/diffraction.pyx b/freud/diffraction.pyx
index b014896d9..912bc7899 100644
--- a/freud/diffraction.pyx
+++ b/freud/diffraction.pyx
@@ -219,7 +219,7 @@ cdef class DiffractionPattern(_Compute):
                 value: True).
 
             For a target view axis in Cartesian coordinates,
-            `rowan <https://rowan.readthedocs.io/en/latest/index.html>__` provides
+            `rowan <https://rowan.readthedocs.io/en/latest/index.html>`__ provides
             ``vector_vector_rotation()``, a function that will return the necessary
             quaternion for rotating the default view orientation to the desired axis.
 

From 04a65830472346825808ac98a229b71e0908d8a5 Mon Sep 17 00:00:00 2001
From: Andrew Kerr <andrewghkerr@gmail.com>
Date: Tue, 15 Jun 2021 19:26:24 -0400
Subject: [PATCH 03/10] Docs edit.

---
 freud/diffraction.pyx | 24 ++++++++++++------------
 1 file changed, 12 insertions(+), 12 deletions(-)

diff --git a/freud/diffraction.pyx b/freud/diffraction.pyx
index 912bc7899..15f22da7f 100644
--- a/freud/diffraction.pyx
+++ b/freud/diffraction.pyx
@@ -218,18 +218,18 @@ cdef class DiffractionPattern(_Compute):
                 the new computations; if False, will accumulate data (Default
                 value: True).
 
-            For a target view axis in Cartesian coordinates,
-            `rowan <https://rowan.readthedocs.io/en/latest/index.html>`__ provides
-            ``vector_vector_rotation()``, a function that will return the necessary
-            quaternion for rotating the default view orientation to the desired axis.
-
-            Example::
-
-                >>> import rowan
-                >>> import numpy as np
-                >>> default_view_axis = np.array([0, 0, 1])
-                >>> desired_view_axis = np.array([1, 1, 1])
-                >>> quaternion = rowan.vector_vector_rotation(desired_view_axis, default_view_axis) # use for view_orientation
+        For a target view axis in Cartesian coordinates,
+        `rowan <https://rowan.readthedocs.io/en/latest/index.html>`__ provides
+        ``vector_vector_rotation()``, a function that will return the necessary
+        quaternion for rotating the default view orientation to the desired axis.
+
+        Example::
+
+            >>> import rowan
+            >>> import numpy as np
+            >>> default_view_axis = np.array([0, 0, 1])
+            >>> desired_view_axis = np.array([1, 1, 1])
+            >>> quaternion = rowan.vector_vector_rotation(desired_view_axis, default_view_axis) # use for view_orientation
 
         """ # noqa: E501
         if reset:

From 8e3adbc9e5a0443328694e6ecc0b1b940b6aba53 Mon Sep 17 00:00:00 2001
From: Andrew Kerr <andrewghkerr@gmail.com>
Date: Tue, 15 Jun 2021 19:36:38 -0400
Subject: [PATCH 04/10] Docs edit.

---
 freud/diffraction.pyx | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/freud/diffraction.pyx b/freud/diffraction.pyx
index 15f22da7f..ee813e596 100644
--- a/freud/diffraction.pyx
+++ b/freud/diffraction.pyx
@@ -229,7 +229,9 @@ cdef class DiffractionPattern(_Compute):
             >>> import numpy as np
             >>> default_view_axis = np.array([0, 0, 1])
             >>> desired_view_axis = np.array([1, 1, 1])
-            >>> quaternion = rowan.vector_vector_rotation(desired_view_axis, default_view_axis) # use for view_orientation
+            >>> quaternion = rowan.vector_vector_rotation(
+                default_view_axis, desired_view_axis
+                ) # use for view_orientation
 
         """ # noqa: E501
         if reset:

From 03761e842270d58ae538883f8680063d333267a3 Mon Sep 17 00:00:00 2001
From: Andrew Kerr <andrewghkerr@gmail.com>
Date: Wed, 16 Jun 2021 12:38:36 -0400
Subject: [PATCH 05/10] Create empty tests w/ TODOs.

This commit includes a rough / hacky change to diffraction.pyx that makes the long z-dimesion test pass by selecting  from only the box vectors in the x and y directions. This will likely be changed in the process of making the other new tests pass.
---
 freud/diffraction.pyx                        |  6 ++-
 tests/test_diffraction_DiffractionPattern.py | 51 ++++++++++++++++++++
 2 files changed, 55 insertions(+), 2 deletions(-)

diff --git a/freud/diffraction.pyx b/freud/diffraction.pyx
index ee813e596..16f884a7d 100644
--- a/freud/diffraction.pyx
+++ b/freud/diffraction.pyx
@@ -169,7 +169,7 @@ cdef class DiffractionPattern(_Compute):
             roll_shift -= 0.5 / zoom
 
         box_matrix = box.to_matrix()
-        ss = np.max(box_matrix) * inv_shear
+        ss = np.max(box_matrix[:, 0:2]) * inv_shear
 
         shift_matrix = np.array(
             [[1, 0, -roll],
@@ -293,7 +293,9 @@ cdef class DiffractionPattern(_Compute):
 
         # Cache the view orientation and box matrix scale factor for
         # lazy evaluation of k-values and k-vectors
-        self._box_matrix_scale_factor = np.max(system.box.to_matrix())
+        self._box_matrix_scale_factor = np.max(
+            rowan.rotate(view_orientation, system.box.to_matrix())[:, 0:2]
+        )
         self._view_orientation = view_orientation
         self._k_scale_factor = 2 * np.pi * self.output_size / (self._box_matrix_scale_factor * zoom)
         self._k_values_cached = False
diff --git a/tests/test_diffraction_DiffractionPattern.py b/tests/test_diffraction_DiffractionPattern.py
index 28720cd0a..7668b99ce 100644
--- a/tests/test_diffraction_DiffractionPattern.py
+++ b/tests/test_diffraction_DiffractionPattern.py
@@ -282,3 +282,54 @@ def test_cubic_system_parameterized(self):
                                 ideal_peaks[peak] = True
 
                     assert all(ideal_peaks.values())
+
+    def test_system_long_and_short_in_view_direction(self):
+        # Same as above tests, but with unit cell lengthed
+        # and shortened in the direction parallel to the view axis.
+        # This should produce a diffraction pattern identical to the
+        # simple cubic case.
+
+        dp = freud.diffraction.DiffractionPattern(grid_size = 512)
+
+        for Lz in [0.5, 2, 3.456]:
+            Lx, Ly = 1, 1
+            cell = freud.box.Box(Lx=Lx, Ly=Ly, Lz=Lz)
+            box, points = freud.data.UnitCell(cell).generate_system(
+                num_replicas=16, scale=1, sigma_noise=0.1
+            )
+
+            dp.compute((box, points), zoom=2, view_orientation=None)
+
+            threshold = 0.2
+            xs, ys = np.nonzero(dp.diffraction > threshold)
+            xy = np.dstack((xs, ys))[0]
+
+            ideal_peaks = {i: False for i in [-2, -1, 0, 1, 2]}
+
+            lattice_vector = np.array([Lx, Ly, Lz])
+            for peak in ideal_peaks:
+                for x, y in xy:
+                    k_vector = dp.k_vectors[x, y]
+                    dot_prod = np.dot(k_vector, lattice_vector)
+
+                    if np.isclose(dot_prod, peak * 2 * np.pi, atol=1e-2):
+                        ideal_peaks[peak] = True
+
+            assert all(ideal_peaks.values())
+
+    def test_unique_box_vector_lengths(self):
+        # TODO: k-values, k-vectors and peaks should be
+        # correct for various boxes with unique side lengths.
+        # Testing w/ variations of other parameters is also possible.
+        pass
+
+    def test_rotated_system(self):
+        # TODO: Same as above, but for boxes that are rotated
+        # such that no face normal vector aligns with the view axis.
+        pass
+
+    def test_multiple_faces_project_equally(self):
+        # TODO: Special case of above, test configurations where two
+        # or three box face normal vectors project equally onto
+        # the view axis.
+        pass

From 12e016f35e208191d61f45a157757448376b94b1 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Wed, 16 Jun 2021 16:38:57 +0000
Subject: [PATCH 06/10] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 tests/test_diffraction_DiffractionPattern.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_diffraction_DiffractionPattern.py b/tests/test_diffraction_DiffractionPattern.py
index 7668b99ce..1187b9021 100644
--- a/tests/test_diffraction_DiffractionPattern.py
+++ b/tests/test_diffraction_DiffractionPattern.py
@@ -289,7 +289,7 @@ def test_system_long_and_short_in_view_direction(self):
         # This should produce a diffraction pattern identical to the
         # simple cubic case.
 
-        dp = freud.diffraction.DiffractionPattern(grid_size = 512)
+        dp = freud.diffraction.DiffractionPattern(grid_size=512)
 
         for Lz in [0.5, 2, 3.456]:
             Lx, Ly = 1, 1

From 7ca562f91b947a5877793a90d0c389539ec84d71 Mon Sep 17 00:00:00 2001
From: Bradley Dice <bdice@bradleydice.com>
Date: Thu, 24 Jun 2021 07:08:03 -0500
Subject: [PATCH 07/10] Add intersphinx reference to vector_vector_rotation.

---
 doc/source/conf.py    | 1 +
 freud/diffraction.pyx | 9 +++++----
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/doc/source/conf.py b/doc/source/conf.py
index 40642d75f..42e9a540e 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -320,6 +320,7 @@
     "MDAnalysis": ("https://docs.mdanalysis.org/stable", None),
     "ovito": ("https://www.ovito.org/docs/current/python/", None),
     "pytest": ("https://docs.pytest.org/en/stable", None),
+    "rowan": ("https://rowan.readthedocs.io/en/stable/", None),
 }
 
 autodoc_default_options = {
diff --git a/freud/diffraction.pyx b/freud/diffraction.pyx
index 16f884a7d..1a4b454ac 100644
--- a/freud/diffraction.pyx
+++ b/freud/diffraction.pyx
@@ -218,10 +218,11 @@ cdef class DiffractionPattern(_Compute):
                 the new computations; if False, will accumulate data (Default
                 value: True).
 
-        For a target view axis in Cartesian coordinates,
-        `rowan <https://rowan.readthedocs.io/en/latest/index.html>`__ provides
-        ``vector_vector_rotation()``, a function that will return the necessary
-        quaternion for rotating the default view orientation to the desired axis.
+        For a target view axis in Cartesian coordinates, `rowan
+        <https://rowan.readthedocs.io/en/latest/index.html>`__ provides
+        :py:func:`rowan.vector_vector_rotation`, a function that will return
+        the necessary quaternion for rotating the default view orientation to
+        the desired axis.
 
         Example::
 

From f7947be38ac20e6083da4a593001845478fff404 Mon Sep 17 00:00:00 2001
From: Bradley Dice <bdice@bradleydice.com>
Date: Thu, 24 Jun 2021 07:09:36 -0500
Subject: [PATCH 08/10] Update code snippet.

---
 freud/diffraction.pyx | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/freud/diffraction.pyx b/freud/diffraction.pyx
index 1a4b454ac..fa1f3ecd2 100644
--- a/freud/diffraction.pyx
+++ b/freud/diffraction.pyx
@@ -230,11 +230,11 @@ cdef class DiffractionPattern(_Compute):
             >>> import numpy as np
             >>> default_view_axis = np.array([0, 0, 1])
             >>> desired_view_axis = np.array([1, 1, 1])
-            >>> quaternion = rowan.vector_vector_rotation(
-                default_view_axis, desired_view_axis
-                ) # use for view_orientation
+            >>> view_orientation = rowan.vector_vector_rotation(
+            ...     default_view_axis, desired_view_axis
+            ... )
 
-        """ # noqa: E501
+        """
         if reset:
             self._diffraction = np.zeros((self.output_size, self.output_size))
             self._frame_counter = 0

From a5fddf91942bc9791de9bf4a1970628cf7aae9d5 Mon Sep 17 00:00:00 2001
From: Andrew Kerr <andrewghkerr@gmail.com>
Date: Mon, 30 Aug 2021 11:50:17 -0400
Subject: [PATCH 09/10] updated tests to cover more systems and better align
 with Tim's implementation

---
 tests/test_diffraction_DiffractionPattern.py | 447 +++++++++++++++----
 1 file changed, 351 insertions(+), 96 deletions(-)

diff --git a/tests/test_diffraction_DiffractionPattern.py b/tests/test_diffraction_DiffractionPattern.py
index 7668b99ce..1ec7307d7 100644
--- a/tests/test_diffraction_DiffractionPattern.py
+++ b/tests/test_diffraction_DiffractionPattern.py
@@ -3,12 +3,13 @@
 import numpy.testing as npt
 import pytest
 import rowan
+import itertools
+import time
 
 import freud
 
 matplotlib.use("agg")
 
-
 class TestDiffractionPattern:
     def test_compute(self):
         dp = freud.diffraction.DiffractionPattern()
@@ -156,68 +157,170 @@ def test_repr(self):
         dp = freud.diffraction.DiffractionPattern(grid_size=123, output_size=234)
         assert str(dp) == str(eval(repr(dp)))
 
-    def test_k_values_and_k_vectors(self):
+    # based on the implimentation in the current Freud release
+    @pytest.mark.parametrize("size, zoom", [
+        (size, zoom) for size in (2, 5, 10) for zoom in (2, 3.4)
+        ])
+    def test_k_values_and_k_vectors(self, size, zoom):
         dp = freud.diffraction.DiffractionPattern()
 
-        for size in [2, 5, 10]:
-            box, positions = freud.data.make_random_system(size, 1)
-            zoom = 4
-            view_orientation = np.asarray([1, 0, 0, 0])
-            dp.compute((box, positions), view_orientation=view_orientation, zoom=zoom)
+        box, positions = freud.data.make_random_system(size, 1)
 
-            output_size = dp.output_size
-            npt.assert_allclose(dp.k_values[output_size // 2], 0)
-            center_index = (output_size // 2, output_size // 2)
-            npt.assert_allclose(dp.k_vectors[center_index], [0, 0, 0])
+        view_orientation = np.asarray([1, 0, 0, 0])
+        dp.compute((box, positions), view_orientation=view_orientation, zoom=zoom)
 
-            # Tests for the first and last k value bins. Uses the math for
-            # np.fft.fftfreq and the _k_scale_factor to write an expression for
-            # the first and last bins' k values and k vectors.
+        output_size = dp.output_size
+        npt.assert_allclose(dp.k_values[output_size // 2], 0)
+        center_index = (output_size // 2, output_size // 2)
+        npt.assert_allclose(dp.k_vectors[center_index], [0, 0, 0])
 
-            scale_factor = 2 * np.pi * output_size / (np.max(box.to_matrix()) * zoom)
+        # Tests for the first and last k value bins. Uses the math for
+        # np.fft.fftfreq and the _k_scale_factor to write an expression for
+        # the first and last bins' k values and k vectors.
 
-            if output_size % 2 == 0:
-                first_k_value = -0.5 * scale_factor
-                last_k_value = (0.5 - (1 / output_size)) * scale_factor
-            else:
-                first_k_value = (-0.5 + 1 / (2 * output_size)) * scale_factor
-                last_k_value = (0.5 - 1 / (2 * output_size)) * scale_factor
+        # TODO:
+        # (1) update np.max(box.to_matrix()) to a new representation
+        # from Tim's diffraction notebook
+        scale_factor = 2 * np.pi * output_size / (np.max(box.to_matrix()) * zoom)
 
-            npt.assert_allclose(dp.k_values[0], first_k_value)
-            npt.assert_allclose(dp.k_values[-1], last_k_value)
+        if output_size % 2 == 0:
+            first_k_value = -0.5 * scale_factor
+            last_k_value = (0.5 - (1 / output_size)) * scale_factor
+        else:
+            first_k_value = (-0.5 + 1 / (2 * output_size)) * scale_factor
+            last_k_value = (0.5 - 1 / (2 * output_size)) * scale_factor
 
-            first_k_vector = rowan.rotate(
-                view_orientation, [first_k_value, first_k_value, 0]
-            )
-            last_k_vector = rowan.rotate(
-                view_orientation, [last_k_value, last_k_value, 0]
-            )
-            top_right_k_vector = rowan.rotate(
-                view_orientation, [first_k_value, last_k_value, 0]
-            )
-            bottom_left_k_vector = rowan.rotate(
-                view_orientation, [last_k_value, first_k_value, 0]
-            )
+        npt.assert_allclose(len(dp.k_values), output_size)
 
-            npt.assert_allclose(dp.k_vectors[0, 0], first_k_vector)
-            npt.assert_allclose(dp.k_vectors[-1, -1], last_k_vector)
-            npt.assert_allclose(dp.k_vectors[0, -1], top_right_k_vector)
-            npt.assert_allclose(dp.k_vectors[-1, 0], bottom_left_k_vector)
+        npt.assert_allclose(dp.k_values[0], first_k_value)
+        npt.assert_allclose(dp.k_values[-1], last_k_value)
 
-            center = output_size // 2
-            top_center_k_vector = rowan.rotate(view_orientation, [0, first_k_value, 0])
-            bottom_center_k_vector = rowan.rotate(
-                view_orientation, [0, last_k_value, 0]
-            )
-            left_center_k_vector = rowan.rotate(view_orientation, [first_k_value, 0, 0])
-            right_center_k_vector = rowan.rotate(view_orientation, [last_k_value, 0, 0])
+        first_k_vector = rowan.rotate(
+            view_orientation, [first_k_value, first_k_value, 0]
+        )
+        last_k_vector = rowan.rotate(
+            view_orientation, [last_k_value, last_k_value, 0]
+        )
+        top_right_k_vector = rowan.rotate(
+            view_orientation, [first_k_value, last_k_value, 0]
+        )
+        bottom_left_k_vector = rowan.rotate(
+            view_orientation, [last_k_value, first_k_value, 0]
+        )
 
-            npt.assert_allclose(dp.k_vectors[center, 0], top_center_k_vector)
-            npt.assert_allclose(dp.k_vectors[center, -1], bottom_center_k_vector)
-            npt.assert_allclose(dp.k_vectors[0, center], left_center_k_vector)
-            npt.assert_allclose(dp.k_vectors[-1, center], right_center_k_vector)
+        npt.assert_allclose(dp.k_vectors[0, 0], first_k_vector)
+        npt.assert_allclose(dp.k_vectors[-1, -1], last_k_vector)
+        npt.assert_allclose(dp.k_vectors[0, -1], top_right_k_vector)
+        npt.assert_allclose(dp.k_vectors[-1, 0], bottom_left_k_vector)
 
-    def test_cubic_system(self):
+        center = output_size // 2
+        top_center_k_vector = rowan.rotate(view_orientation, [0, first_k_value, 0])
+        bottom_center_k_vector = rowan.rotate(
+            view_orientation, [0, last_k_value, 0]
+        )
+        left_center_k_vector = rowan.rotate(view_orientation, [first_k_value, 0, 0])
+        right_center_k_vector = rowan.rotate(view_orientation, [last_k_value, 0, 0])
+
+        npt.assert_allclose(dp.k_vectors[center, 0], top_center_k_vector)
+        npt.assert_allclose(dp.k_vectors[center, -1], bottom_center_k_vector)
+        npt.assert_allclose(dp.k_vectors[0, center], left_center_k_vector)
+        npt.assert_allclose(dp.k_vectors[-1, center], right_center_k_vector)
+
+    # based on Tim's new diffraction code
+    @pytest.mark.parametrize("Lx, Ly, Lz", [(16, 10, 12)])
+    def test_k_values_and_vectors_non_cubic(self, Lx, Ly, Lz):
+
+        dp = freud.diffraction.DiffractionPattern()
+
+        size = 30
+        _, positions = freud.data.make_random_system(size, 1)
+        box = freud.box.Box(Lx=Lx, Ly=Ly, Lz=Lz)
+        positions = box.wrap(positions)
+
+        view_orientation = np.asarray([1, 0, 0, 0])
+        zoom = 2
+        dp.compute((box, positions), view_orientation=view_orientation, zoom=zoom)
+
+        grid_size = dp.grid_size
+        output_size = dp.output_size
+
+        # new method for selecting the optimal axes for the diffraction plot,
+        # regardless of system shape or rotation. Particles are coerced into
+        # an orthorhombic new_box, whose edges perpendicular to view direction
+        # are used to determine the k-bins in x- and y-directions.
+        v1, v2, v3 = box.to_matrix().T
+        v1 = rowan.rotate(view_orientation, v1)
+        v2 = rowan.rotate(view_orientation, v2)
+        v3 = rowan.rotate(view_orientation, v3)
+        v1[2] = 0
+        v2[2] = 0
+        v3[2] = 0
+
+        pos_z = np.array([0, 0, 1])
+        f3 = np.linalg.norm(np.dot(np.cross(v1, v2), pos_z))
+        f2 = np.linalg.norm(np.dot(np.cross(v1, v3), pos_z))
+        f1 = np.linalg.norm(np.dot(np.cross(v2, v3), pos_z))
+        max_idx = np.argmax([f1, f2, f3])
+
+        if max_idx == 0:
+            v1_proj = v2
+            v2_proj = v3
+        elif max_idx == 1:
+            v1_proj = v3
+            v2_proj = v1
+        elif max_idx == 2:
+            v1_proj = v1
+            v2_proj = v2
+        v1_proj_original = np.copy(v1_proj)
+        v2_proj_original = np.copy(v2_proj)
+
+        rot_angle = -np.arctan2(v1_proj_original[1], v1_proj_original[0])
+        system_rotation = rowan.from_axis_angle([0,0,1], rot_angle)
+        v1_proj = rowan.rotate(system_rotation, v1_proj_original)
+        v2_proj = rowan.rotate(system_rotation, v2_proj_original)
+        # make sure system is right-handed, which means the y-component of
+        # v2_proj is positive; if not, place v2_proj along the x-axis
+        # note that if this is the case, we switch which vectors we call v1 and v2
+        # so that what we call v1 is the one that lies along +x
+        if v2_proj[1] < 0:
+            rot_angle = np.arctan2(v2_proj_original[1], v2_proj_original[0])
+            system_rotation = rowan.from_axis_angle([0,0,1], rot_angle)
+            v1_proj = rowan.rotate(system_rotation, v2_proj_original)
+            v2_proj = rowan.rotate(system_rotation, v1_proj_original)
+        assert v2_proj[1] > 0
+
+        new_box = freud.Box(Lx=v1_proj[0], Ly=v2_proj[1])
+
+        # TODO: verify this
+        kx_scale_factor = 2 * np.pi * grid_size / (new_box.Lx * zoom)
+        ky_scale_factor = 2 * np.pi * grid_size / (new_box.Ly * zoom)
+
+        if output_size % 2 == 0:
+            first_kx_value = -0.5 * kx_scale_factor
+            first_ky_value = -0.5 * ky_scale_factor
+            last_kx_value = (0.5 - (1 / output_size)) * kx_scale_factor
+            last_ky_value = (0.5 - (1 / output_size)) * ky_scale_factor
+        else:
+            first_kx_value = (-0.5 + 1 / (2 * output_size)) * kx_scale_factor
+            first_ky_value = (-0.5 + 1 / (2 * output_size)) * ky_scale_factor
+            last_kx_value = (0.5 - 1 / (2 * output_size)) * kx_scale_factor
+            last_ky_value = (0.5 - 1 / (2 * output_size)) * ky_scale_factor
+
+        # TODO: from Tim's notebook:
+        # kxs = np.fft.fftshift(np.fft.fftfreq(n=output_size, d=new_box.Lx/grid_size))
+        # kys = np.fft.fftshift(np.fft.fftfreq(n=output_size, d=new_box.Ly/grid_size))
+        # k_vectors = np.asarray(np.meshgrid(kxs, kys, [0])).T[0] / zoom * 2 * np.pi
+        # i.e. scale factor = 2 * np.pi * grid_size / (new_box.{axis} * zoom)
+
+        # tests to verify kx and ky bins go here. These are difficult to write
+        # until the dp object has separate arrays for kx and ky bins. See
+        # temp_k_test.py
+
+
+    # TODO: add tests for each of:
+    # diamond(?), and
+    # non-orthorhombic (sheared bcc?) lattices-
+    def test_sc_system(self):
         length = 1
         box, positions = freud.data.UnitCell.sc().generate_system(
             num_replicas=16, scale=length, sigma_noise=0.1 * length
@@ -228,11 +331,10 @@ def test_cubic_system(self):
         dp.compute((box, positions), zoom=4.123)
 
         # Locate brightest areas of diffraction pattern
-        # (intensity > threshold), and check that the ideal
-        # diffraction peak locations, given by k * R = 2*pi*N
-        # for some lattice vector R and integer N, are contained
-        # within these regions.
-        # This test only checks N in range [-2, 2].
+        # (intensity > threshold), and check that the ideal diffraction peak
+        # locations, given by k * R = 2*pi*N for some lattice vector R and
+        # integer N, are contained within these regions. This test only checks
+        # N in range [-2, 2].
         threshold = 0.2
         xs, ys = np.nonzero(dp.diffraction > threshold)
         xy = np.dstack((xs, ys))[0]
@@ -250,46 +352,78 @@ def test_cubic_system(self):
 
         assert all(ideal_peaks.values())
 
-    def test_cubic_system_parameterized(self):
+    # TODO: currently fails for Lx, Ly, Lz = 16, 14, 15, and other
+    # combinations of lengths that differ by 1. Running the same
+    # failing tuple multiple times produces nearly idential errors
+    # (same number of mismatched elements, nearly same maximum
+    # rel. and abs. difference).
+    @pytest.mark.parametrize("Lx, Ly, Lz",
+    [(16, 16, 16), (16, 8, 4), (16, 14, 15)
+    ])
+    def test_ideal_gas(self, Lx, Ly, Lz):
+
+        dp = freud.diffraction.DiffractionPattern()
+        box = freud.box.Box(Lx=Lx, Ly=Ly, Lz=Lz)
+        zoom = 1
+
+        for i in range(1, 20):
+            milliseconds = int(round(time.time() * 1000))
+
+            _, points = freud.data.make_random_system(
+                16.0, 16**5, seed=(i*milliseconds)%(2**32-1)
+            )
+            wrapped_points = box.wrap(points)
+            dp.compute((box, wrapped_points), zoom=zoom, reset=False)
+
+        cutoff = 0.05
+        avg = np.mean(dp.diffraction[dp.diffraction < cutoff])
+
+        npt.assert_allclose(dp.diffraction[dp.diffraction < cutoff], avg, atol=1e-4)
+
+    @pytest.mark.parametrize("grid_size, output_size, zoom", [
+        (grid_size, output_size, zoom)
+        for grid_size in (256, 1024)
+        for output_size in (255,256, 1023, 1024)
+        for zoom in (1, 2.5, 4.123)
+    ])
+    def test_cubic_system_parameterized(self, grid_size, output_size, zoom):
         length = 1
         box, positions = freud.data.UnitCell.sc().generate_system(
             num_replicas=16, scale=length, sigma_noise=0.1 * length
         )
-        # Same as above test but with different grid_size,
-        # output_size, and zoom values.
-        for grid_size in (256, 1024):
-            for output_size in (255, 256, 1023, 1024):
-                for zoom in (1, 2.5, 4.123):
-                    dp = freud.diffraction.DiffractionPattern(
-                        grid_size=grid_size,
-                        output_size=output_size,
-                    )
-                    dp.compute((box, positions), zoom=zoom)
 
-                    threshold = 0.2
-                    xs, ys = np.nonzero(dp.diffraction > threshold)
-                    xy = np.dstack((xs, ys))[0]
+        dp = freud.diffraction.DiffractionPattern(
+            grid_size=grid_size,
+            output_size=output_size,
+        )
+        dp.compute((box, positions), zoom=zoom)
 
-                    ideal_peaks = {i: False for i in [-2, -1, 0, 1, 2]}
+        threshold = 0.2
+        xs, ys = np.nonzero(dp.diffraction > threshold)
+        xy = np.dstack((xs, ys))[0]
 
-                    lattice_vector = np.array([length, length, length])
-                    for peak in ideal_peaks:
-                        for x, y in xy:
-                            k_vector = dp.k_vectors[x, y]
-                            dot_prod = np.dot(k_vector, lattice_vector)
+        ideal_peaks = {i: False for i in [-2, -1, 0, 1, 2]}
 
-                            if np.isclose(dot_prod, peak * 2 * np.pi, atol=1e-2):
-                                ideal_peaks[peak] = True
+        lattice_vector = np.array([length, length, length])
+        for peak in ideal_peaks:
+            for x, y in xy:
+                k_vector = dp.k_vectors[x, y]
+                dot_prod = np.dot(k_vector, lattice_vector)
 
-                    assert all(ideal_peaks.values())
+                if np.isclose(dot_prod, peak * 2 * np.pi, atol=1e-2):
+                    ideal_peaks[peak] = True
+
+        assert all(ideal_peaks.values())
 
     def test_system_long_and_short_in_view_direction(self):
-        # Same as above tests, but with unit cell lengthed
-        # and shortened in the direction parallel to the view axis.
+        # Here, the unit cell is lengthed and shortened in
+        # the direction parallel to the view axis.
         # This should produce a diffraction pattern identical to the
         # simple cubic case.
 
-        dp = freud.diffraction.DiffractionPattern(grid_size = 512)
+        # TODO: compare w/ Tim's notebook / diffraction implementation
+
+        dp = freud.diffraction.DiffractionPattern(grid_size=512)
 
         for Lz in [0.5, 2, 3.456]:
             Lx, Ly = 1, 1
@@ -317,19 +451,140 @@ def test_system_long_and_short_in_view_direction(self):
 
             assert all(ideal_peaks.values())
 
-    def test_unique_box_vector_lengths(self):
-        # TODO: k-values, k-vectors and peaks should be
-        # correct for various boxes with unique side lengths.
-        # Testing w/ variations of other parameters is also possible.
-        pass
+    @pytest.mark.parametrize("Lx, Ly, Lz", [permutation for permutation in itertools.permutations((0.6, 1.1, 2))])
+    def test_unique_box_vector_lengths(self, Lx, Ly, Lz):
+        dp = freud.diffraction.DiffractionPattern(grid_size=512)
+
+        cell = freud.box.Box(Lx, Ly, Lz)
+        box, points = freud.data.UnitCell(cell).generate_system(
+            num_replicas=16, scale=1, sigma_noise=0.1
+        )
+
+        dp.compute((box, points), zoom=2, view_orientation=None)
+
+        threshold = 0.2
+        xs, ys = np.nonzero(dp.diffraction > threshold)
+        xy = np.dstack((xs, ys))[0]
+
+        ideal_peaks = {i: False for i in [-2, -1, 0, 1, 2]}
+
+        lattice_vector = np.array([Lx, Ly, Lz])
+        for peak in ideal_peaks:
+            for x, y in xy:
+                k_vector = dp.k_vectors[x, y]
+                dot_prod = np.dot(k_vector, lattice_vector)
+
+                if np.isclose(dot_prod, peak * 2 * np.pi, atol=1e-2):
+                    ideal_peaks[peak] = True
+
+        assert all(ideal_peaks.values())
+
+    # TODO: verify the chosen fcc, bcc, hex lattice vectors
+    def test_fcc_system(self):
+        length = 1
+        box, positions = freud.data.UnitCell.fcc().generate_system(
+            num_replicas=16, scale=length, sigma_noise=0.1 * length
+        )
+        # Pick a non-integer value for zoom, to ensure that peaks besides k=0
+        # are not perfectly aligned on pixels.
+        dp = freud.diffraction.DiffractionPattern(grid_size = 512)
+        dp.compute((box, positions), zoom = 4.123)
+
+        threshold = 0.2
+        xs, ys = np.nonzero(dp.diffraction > threshold)
+        xy = np.dstack((xs, ys))[0]
+
+        ideal_peaks = {i: False for i in [-2, -1, 0, 1, 2]}
+
+        lattice_vector = np.array([length / 2, length / 2, 0])
+        for peak in ideal_peaks:
+            for x, y in xy:
+                k_vector = dp.k_vectors[x, y]
+                dot_prod = np.dot(k_vector, lattice_vector)
+
+                if np.isclose(dot_prod, peak * 2 * np.pi, atol=1e-2):
+                    ideal_peaks[peak] = True
+
+        assert all(ideal_peaks.values())
+
+    def test_bcc_system(self):
+        length = 1
+        box, positions = freud.data.UnitCell.bcc().generate_system(
+            num_replicas=16, scale=length, sigma_noise=0.1 * length
+        )
+        # Pick a non-integer value for zoom, to ensure that peaks besides k=0
+        # are not perfectly aligned on pixels.
+        dp = freud.diffraction.DiffractionPattern(grid_size = 512)
+        dp.compute((box, positions), zoom = 4.123)
+
+        threshold = 0.2
+        xs, ys = np.nonzero(dp.diffraction > threshold)
+        xy = np.dstack((xs, ys))[0]
+
+        ideal_peaks = {i: False for i in [-2, -1, 0, 1, 2]}
+
+        lattice_vector = np.array([length / 2, length / 2, length / 2])
+        for peak in ideal_peaks:
+            for x, y in xy:
+                k_vector = dp.k_vectors[x, y]
+                dot_prod = np.dot(k_vector, lattice_vector)
+
+                if np.isclose(dot_prod, peak * 2 * np.pi, atol=1e-2):
+                    ideal_peaks[peak] = True
+
+        assert all(ideal_peaks.values())
+
+    def test_hexagonal_system(self):
+        length = 1
+        box, positions = freud.data.UnitCell.hex().generate_system(
+            num_replicas=16, scale=length, sigma_noise=0.1 * length
+        )
+        # Pick a non-integer value for zoom, to ensure that peaks besides k=0
+        # are not perfectly aligned on pixels.
+        dp = freud.diffraction.DiffractionPattern(grid_size = 512)
+        dp.compute((box, positions), zoom = 4.123)
+
+        threshold = 0.2
+        xs, ys = np.nonzero(dp.diffraction > threshold)
+        xy = np.dstack((xs, ys))[0]
+
+        ideal_peaks = {i: False for i in [-2, -1, 0, 1, 2]}
+
+        lattice_vector = np.array([length * np.cos(np.pi / 3), length * np.sin(np.pi / 3), 0])
+        for peak in ideal_peaks:
+            for x, y in xy:
+                k_vector = dp.k_vectors[x, y]
+                dot_prod = np.dot(k_vector, lattice_vector)
+
+                if np.isclose(dot_prod, peak * 2 * np.pi, atol=1e-2):
+                    ideal_peaks[peak] = True
+
+        assert all(ideal_peaks.values())
+
 
     def test_rotated_system(self):
-        # TODO: Same as above, but for boxes that are rotated
-        # such that no face normal vector aligns with the view axis.
-        pass
-
-    def test_multiple_faces_project_equally(self):
-        # TODO: Special case of above, test configurations where two
-        # or three box face normal vectors project equally onto
-        # the view axis.
-        pass
+        dp = freud.diffraction.DiffractionPattern(grid_size=512)
+        box, points = freud.data.UnitCell.sc().generate_system(
+            num_replicas=16, scale=1, sigma_noise=0.1
+        )
+
+        view_orientation = rowan.vector_vector_rotation([0, 0, 1], [0.5, 1, 2])
+
+        dp.compute((box, points), zoom=2, view_orientation=view_orientation)
+
+        threshold = 0.2
+        xs, ys = np.nonzero(dp.diffraction > threshold)
+        xy = np.dstack((xs, ys))[0]
+
+        ideal_peaks = {i: False for i in [-2, -1, 0, 1, 2]}
+
+        lattice_vector = rowan.rotate(view_orientation, [1, 1, 1])
+        for peak in ideal_peaks:
+            for x, y in xy:
+                k_vector = dp.k_vectors[x, y]
+                dot_prod = np.dot(k_vector, lattice_vector)
+
+                if np.isclose(dot_prod, peak * 2 * np.pi, atol=1e-2):
+                    ideal_peaks[peak] = True
+
+        assert all(ideal_peaks.values())

From 732f63f7bc0df02d2db04d0f49a523d079cfb5bf Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Mon, 30 Aug 2021 15:58:17 +0000
Subject: [PATCH 10/10] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 tests/test_diffraction_DiffractionPattern.py | 68 ++++++++++----------
 1 file changed, 35 insertions(+), 33 deletions(-)

diff --git a/tests/test_diffraction_DiffractionPattern.py b/tests/test_diffraction_DiffractionPattern.py
index 1fe4ea7da..131983929 100644
--- a/tests/test_diffraction_DiffractionPattern.py
+++ b/tests/test_diffraction_DiffractionPattern.py
@@ -1,15 +1,17 @@
+import itertools
+import time
+
 import matplotlib
 import numpy as np
 import numpy.testing as npt
 import pytest
 import rowan
-import itertools
-import time
 
 import freud
 
 matplotlib.use("agg")
 
+
 class TestDiffractionPattern:
     def test_compute(self):
         dp = freud.diffraction.DiffractionPattern()
@@ -158,9 +160,9 @@ def test_repr(self):
         assert str(dp) == str(eval(repr(dp)))
 
     # based on the implimentation in the current Freud release
-    @pytest.mark.parametrize("size, zoom", [
-        (size, zoom) for size in (2, 5, 10) for zoom in (2, 3.4)
-        ])
+    @pytest.mark.parametrize(
+        "size, zoom", [(size, zoom) for size in (2, 5, 10) for zoom in (2, 3.4)]
+    )
     def test_k_values_and_k_vectors(self, size, zoom):
         dp = freud.diffraction.DiffractionPattern()
 
@@ -198,9 +200,7 @@ def test_k_values_and_k_vectors(self, size, zoom):
         first_k_vector = rowan.rotate(
             view_orientation, [first_k_value, first_k_value, 0]
         )
-        last_k_vector = rowan.rotate(
-            view_orientation, [last_k_value, last_k_value, 0]
-        )
+        last_k_vector = rowan.rotate(view_orientation, [last_k_value, last_k_value, 0])
         top_right_k_vector = rowan.rotate(
             view_orientation, [first_k_value, last_k_value, 0]
         )
@@ -215,9 +215,7 @@ def test_k_values_and_k_vectors(self, size, zoom):
 
         center = output_size // 2
         top_center_k_vector = rowan.rotate(view_orientation, [0, first_k_value, 0])
-        bottom_center_k_vector = rowan.rotate(
-            view_orientation, [0, last_k_value, 0]
-        )
+        bottom_center_k_vector = rowan.rotate(view_orientation, [0, last_k_value, 0])
         left_center_k_vector = rowan.rotate(view_orientation, [first_k_value, 0, 0])
         right_center_k_vector = rowan.rotate(view_orientation, [last_k_value, 0, 0])
 
@@ -275,7 +273,7 @@ def test_k_values_and_vectors_non_cubic(self, Lx, Ly, Lz):
         v2_proj_original = np.copy(v2_proj)
 
         rot_angle = -np.arctan2(v1_proj_original[1], v1_proj_original[0])
-        system_rotation = rowan.from_axis_angle([0,0,1], rot_angle)
+        system_rotation = rowan.from_axis_angle([0, 0, 1], rot_angle)
         v1_proj = rowan.rotate(system_rotation, v1_proj_original)
         v2_proj = rowan.rotate(system_rotation, v2_proj_original)
         # make sure system is right-handed, which means the y-component of
@@ -284,7 +282,7 @@ def test_k_values_and_vectors_non_cubic(self, Lx, Ly, Lz):
         # so that what we call v1 is the one that lies along +x
         if v2_proj[1] < 0:
             rot_angle = np.arctan2(v2_proj_original[1], v2_proj_original[0])
-            system_rotation = rowan.from_axis_angle([0,0,1], rot_angle)
+            system_rotation = rowan.from_axis_angle([0, 0, 1], rot_angle)
             v1_proj = rowan.rotate(system_rotation, v2_proj_original)
             v2_proj = rowan.rotate(system_rotation, v1_proj_original)
         assert v2_proj[1] > 0
@@ -316,7 +314,6 @@ def test_k_values_and_vectors_non_cubic(self, Lx, Ly, Lz):
         # until the dp object has separate arrays for kx and ky bins. See
         # temp_k_test.py
 
-
     # TODO: add tests for each of:
     # diamond(?), and
     # non-orthorhombic (sheared bcc?) lattices-
@@ -357,9 +354,7 @@ def test_sc_system(self):
     # failing tuple multiple times produces nearly idential errors
     # (same number of mismatched elements, nearly same maximum
     # rel. and abs. difference).
-    @pytest.mark.parametrize("Lx, Ly, Lz",
-    [(16, 16, 16), (16, 8, 4), (16, 14, 15)
-    ])
+    @pytest.mark.parametrize("Lx, Ly, Lz", [(16, 16, 16), (16, 8, 4), (16, 14, 15)])
     def test_ideal_gas(self, Lx, Ly, Lz):
 
         dp = freud.diffraction.DiffractionPattern()
@@ -370,7 +365,7 @@ def test_ideal_gas(self, Lx, Ly, Lz):
             milliseconds = int(round(time.time() * 1000))
 
             _, points = freud.data.make_random_system(
-                16.0, 16**5, seed=(i*milliseconds)%(2**32-1)
+                16.0, 16 ** 5, seed=(i * milliseconds) % (2 ** 32 - 1)
             )
             wrapped_points = box.wrap(points)
             dp.compute((box, wrapped_points), zoom=zoom, reset=False)
@@ -380,12 +375,15 @@ def test_ideal_gas(self, Lx, Ly, Lz):
 
         npt.assert_allclose(dp.diffraction[dp.diffraction < cutoff], avg, atol=1e-4)
 
-    @pytest.mark.parametrize("grid_size, output_size, zoom", [
-        (grid_size, output_size, zoom)
-        for grid_size in (256, 1024)
-        for output_size in (255,256, 1023, 1024)
-        for zoom in (1, 2.5, 4.123)
-    ])
+    @pytest.mark.parametrize(
+        "grid_size, output_size, zoom",
+        [
+            (grid_size, output_size, zoom)
+            for grid_size in (256, 1024)
+            for output_size in (255, 256, 1023, 1024)
+            for zoom in (1, 2.5, 4.123)
+        ],
+    )
     def test_cubic_system_parameterized(self, grid_size, output_size, zoom):
         length = 1
         box, positions = freud.data.UnitCell.sc().generate_system(
@@ -449,7 +447,10 @@ def test_system_long_and_short_in_view_direction(self):
 
             assert all(ideal_peaks.values())
 
-    @pytest.mark.parametrize("Lx, Ly, Lz", [permutation for permutation in itertools.permutations((0.6, 1.1, 2))])
+    @pytest.mark.parametrize(
+        "Lx, Ly, Lz",
+        [permutation for permutation in itertools.permutations((0.6, 1.1, 2))],
+    )
     def test_unique_box_vector_lengths(self, Lx, Ly, Lz):
         dp = freud.diffraction.DiffractionPattern(grid_size=512)
 
@@ -485,8 +486,8 @@ def test_fcc_system(self):
         )
         # Pick a non-integer value for zoom, to ensure that peaks besides k=0
         # are not perfectly aligned on pixels.
-        dp = freud.diffraction.DiffractionPattern(grid_size = 512)
-        dp.compute((box, positions), zoom = 4.123)
+        dp = freud.diffraction.DiffractionPattern(grid_size=512)
+        dp.compute((box, positions), zoom=4.123)
 
         threshold = 0.2
         xs, ys = np.nonzero(dp.diffraction > threshold)
@@ -512,8 +513,8 @@ def test_bcc_system(self):
         )
         # Pick a non-integer value for zoom, to ensure that peaks besides k=0
         # are not perfectly aligned on pixels.
-        dp = freud.diffraction.DiffractionPattern(grid_size = 512)
-        dp.compute((box, positions), zoom = 4.123)
+        dp = freud.diffraction.DiffractionPattern(grid_size=512)
+        dp.compute((box, positions), zoom=4.123)
 
         threshold = 0.2
         xs, ys = np.nonzero(dp.diffraction > threshold)
@@ -539,8 +540,8 @@ def test_hexagonal_system(self):
         )
         # Pick a non-integer value for zoom, to ensure that peaks besides k=0
         # are not perfectly aligned on pixels.
-        dp = freud.diffraction.DiffractionPattern(grid_size = 512)
-        dp.compute((box, positions), zoom = 4.123)
+        dp = freud.diffraction.DiffractionPattern(grid_size=512)
+        dp.compute((box, positions), zoom=4.123)
 
         threshold = 0.2
         xs, ys = np.nonzero(dp.diffraction > threshold)
@@ -548,7 +549,9 @@ def test_hexagonal_system(self):
 
         ideal_peaks = {i: False for i in [-2, -1, 0, 1, 2]}
 
-        lattice_vector = np.array([length * np.cos(np.pi / 3), length * np.sin(np.pi / 3), 0])
+        lattice_vector = np.array(
+            [length * np.cos(np.pi / 3), length * np.sin(np.pi / 3), 0]
+        )
         for peak in ideal_peaks:
             for x, y in xy:
                 k_vector = dp.k_vectors[x, y]
@@ -559,7 +562,6 @@ def test_hexagonal_system(self):
 
         assert all(ideal_peaks.values())
 
-
     def test_rotated_system(self):
         dp = freud.diffraction.DiffractionPattern(grid_size=512)
         box, points = freud.data.UnitCell.sc().generate_system(