diff --git a/cmeutils/__init__.py b/cmeutils/__init__.py
index 2fd5818..b841d21 100644
--- a/cmeutils/__init__.py
+++ b/cmeutils/__init__.py
@@ -1,4 +1,4 @@
-from . import gsd_utils
+from . import gsd_utils, structure
 from .__version__ import __version__
 
-__all__ = ["__version__", "gsd_utils"]
+__all__ = ["__version__", "gsd_utils", "structure"]
diff --git a/cmeutils/structure.py b/cmeutils/structure.py
index 55b3e2b..1d4f3ac 100644
--- a/cmeutils/structure.py
+++ b/cmeutils/structure.py
@@ -1,8 +1,53 @@
-from cmeutils import gsd_utils
 import freud
 import gsd
 import gsd.hoomd
 import numpy as np
+from rowan import vector_vector_rotation
+
+from cmeutils import gsd_utils
+
+
+def get_quaternions(n_views = 20):
+    """Get the quaternions for the specified number of views.
+
+    The first (n_view - 3) views will be the views even distributed on a sphere,
+    while the last three views will be the face-on, edge-on, and corner-on
+    views, respectively.
+
+    These quaternions are useful as input to `view_orientation` kwarg in
+    `freud.diffraction.Diffractometer.compute`.
+
+    Parameters
+    ----------
+    n_views : int, default 20
+        The number of views to compute.
+
+    Returns
+    -------
+    list of numpy.ndarray
+        Quaternions as (4,) arrays.
+    """
+    if n_views <=3 or not isinstance(n_views, int):
+        raise ValueError("Please set n_views to an integer greater than 3.")
+    # Calculate points for even distribution on a sphere
+    ga = np.pi * (3 - 5**0.5)
+    theta = ga * np.arange(n_views-3)
+    z = np.linspace(1 - 1/(n_views-3), 1/(n_views-3), n_views-3)
+    radius = np.sqrt(1 - z * z)
+    points = np.zeros((n_views, 3))
+    points[:-3,0] = radius * np.cos(theta)
+    points[:-3,1] = radius * np.sin(theta)
+    points[:-3,2] = z
+
+    # face on
+    points[-3] = np.array([0, 0, 1])
+    # edge on
+    points[-2] = np.array([0, 1, 1])
+    # corner on
+    points[-1] = np.array([1, 1, 1])
+
+    unit_z = np.array([0, 0, 1])
+    return [vector_vector_rotation(i, unit_z) for i in points]
 
 
 def gsd_rdf(
diff --git a/cmeutils/tests/test_structure.py b/cmeutils/tests/test_structure.py
index 87f1b41..1a42d8a 100644
--- a/cmeutils/tests/test_structure.py
+++ b/cmeutils/tests/test_structure.py
@@ -1,5 +1,7 @@
 import pytest
 
+import numpy as np
+
 from cmeutils.tests.base_test import BaseTest
 
 
@@ -10,7 +12,7 @@ def test_gsd_rdf(self, gsdfile_bond):
         rdf_ex, norm = gsd_rdf(gsdfile_bond, "A", "B")
         rdf_noex, norm2 = gsd_rdf(gsdfile_bond, "A", "B", exclude_bonded=False)
         assert norm2 == 1
-        assert rdf_noex != rdf_ex
+        assert not np.array_equal(rdf_noex, rdf_ex)
 
     def test_gsd_rdf_samename(self, gsdfile_bond):
         from cmeutils.structure import gsd_rdf
@@ -18,4 +20,17 @@ def test_gsd_rdf_samename(self, gsdfile_bond):
         rdf_ex, norm = gsd_rdf(gsdfile_bond, "A", "A")
         rdf_noex, norm2 = gsd_rdf(gsdfile_bond, "A", "A", exclude_bonded=False)
         assert norm2 == 1
-        assert rdf_noex != rdf_ex
+        assert not np.array_equal(rdf_noex, rdf_ex)
+
+    def test_get_quaternions(self):
+        from cmeutils.structure import get_quaternions
+
+        with pytest.raises(ValueError):
+            get_quaternions(0)
+
+        with pytest.raises(ValueError):
+            get_quaternions(5.3)
+
+        qs = get_quaternions()
+        assert len(qs) == 20
+        assert len(qs[0]) == 4
diff --git a/environment.yml b/environment.yml
index d6f2f03..baae63c 100644
--- a/environment.yml
+++ b/environment.yml
@@ -7,5 +7,7 @@ dependencies:
   - numpy
   - pip
   - python=3.7
+  - matplotlib
   - pytest
   - pytest-cov
+  - rowan