diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 422001df..78a54030 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -11,6 +11,7 @@ Unreleased
 
 Added
 -----
+- ``Vector3d.from_path_ends()`` class method to get vectors between two vectors.
 
 Changed
 -------
@@ -23,6 +24,8 @@ Removed
 
 Fixed
 -----
+- Transparency of polar stereographic grid lines can now be controlled by Matplotlib's
+  ``grid.alpha``, just like the azimuth grid lines.
 
 2023-03-14 - version 0.11.1
 ===========================
diff --git a/examples/rotations/README.txt b/examples/rotations/README.txt
new file mode 100644
index 00000000..d7c891a3
--- /dev/null
+++ b/examples/rotations/README.txt
@@ -0,0 +1,2 @@
+Rotations
+=========
\ No newline at end of file
diff --git a/examples/rotations/combine_rotations.py b/examples/rotations/combine_rotations.py
new file mode 100644
index 00000000..6f20230a
--- /dev/null
+++ b/examples/rotations/combine_rotations.py
@@ -0,0 +1,83 @@
+r"""
+===================
+Combining rotations
+===================
+
+This example demonstrates how to combine two rotations :math:`g_A` and
+:math:`g_B`, i.e. from left to right like so
+
+.. math::
+    g_{AB} = g_A \cdot g_B.
+
+This order follows from the convention of passive rotations chosen in
+orix which follows :cite:`rowenhorst2015consistent`.
+
+To convince ourselves that this order is correct, we will reproduce the
+example given by Rowenhorst and co-workers in section 4.2.2 in the above
+mentioned paper. We want to rotate a vector :math:`(0, 0, z)` by two
+rotations: rotation :math:`A` by :math:`120^{\circ}` around
+:math:`[1 1 1]`, and rotation :math:`B` by :math:`180^{\circ}` around
+:math:`[1 1 0]`; rotation :math:`A` will be carried out first, followed
+by rotation :math:`B`.
+
+Note that a negative angle when *defining* a rotation in the axis-angle
+representation is necessary for consistent transformations between
+rotation representations. The rotation still rotates a vector
+intuitively.
+"""
+
+import matplotlib.pyplot as plt
+
+from orix import plot
+from orix.quaternion import Rotation
+from orix.vector import Vector3d
+
+plt.rcParams.update({"font.size": 12, "grid.alpha": 0.5})
+
+gA = Rotation.from_axes_angles([1, 1, 1], -120, degrees=True)
+gB = Rotation.from_axes_angles([1, 1, 0], -180, degrees=True)
+gAB = gA * gB
+
+# Compare with quaternions and orientation matrices from section 4.2.2
+# in Rowenhorst et al. (2015)
+g_all = Rotation.stack((gA, gB, gAB)).squeeze()
+print("gA, gB and gAB:\n* As quaternions:\n", g_all)
+print("* As orientation matrices:\n", g_all.to_matrix().squeeze().round(10))
+
+v_start = Vector3d.zvector()
+v_end = gAB * v_start
+print(
+    "Point rotated by gAB:\n",
+    v_start.data.squeeze().tolist(),
+    "->",
+    v_end.data.squeeze().round(10).tolist(),
+)
+
+# Illustrate the steps of the rotation by plotting the vector before
+# (red), during (green) and after (blue) the rotation and the rotation
+# paths (first: cyan; second: magenta)
+v_intermediate = gB * v_start
+
+v_si_path = Vector3d.from_path_ends(Vector3d.stack((v_start, v_intermediate)))
+v_sie_path = Vector3d.from_path_ends(Vector3d.stack((v_intermediate, v_end)))
+
+fig = plt.figure(layout="tight")
+ax0 = fig.add_subplot(121, projection="stereographic", hemisphere="upper")
+ax1 = fig.add_subplot(122, projection="stereographic", hemisphere="lower")
+ax0.stereographic_grid(), ax1.stereographic_grid()
+Vector3d.stack((v_start, v_intermediate, v_end)).scatter(
+    figure=fig,
+    s=50,
+    c=["r", "g", "b"],
+    axes_labels=["e1", "e2"],
+)
+ax0.plot(v_si_path, color="c"), ax1.plot(v_si_path, color="c")
+ax0.plot(v_sie_path, color="m"), ax1.plot(v_sie_path, color="m")
+gA.axis.scatter(figure=fig, c="orange")
+gB.axis.scatter(figure=fig, c="k")
+text_kw = dict(bbox=dict(alpha=0.5, fc="w", boxstyle="round,pad=0.1"), ha="right")
+ax0.text(v_start, s="Start", **text_kw)
+ax1.text(v_intermediate, s="Intermediate", **text_kw)
+ax1.text(v_end, s="End", **text_kw)
+ax1.text(gA.axis, s="Axis gA", **text_kw)
+ax0.text(gB.axis, s="Axis gB", **text_kw)
diff --git a/orix/plot/stereographic_plot.py b/orix/plot/stereographic_plot.py
index 7b344a8a..3ea28fd4 100644
--- a/orix/plot/stereographic_plot.py
+++ b/orix/plot/stereographic_plot.py
@@ -746,6 +746,7 @@ def _polar_grid(self, resolution: Optional[float] = None):
             label=label,
             edgecolors=kwargs["ec"],
             facecolors=kwargs["fc"],
+            alpha=kwargs["alpha"],
         )
         has_collection, index = self._has_collection(label, self.collections)
         if has_collection:
diff --git a/orix/tests/plot/test_stereographic_plot.py b/orix/tests/plot/test_stereographic_plot.py
index fbb9d7b5..636965b7 100644
--- a/orix/tests/plot/test_stereographic_plot.py
+++ b/orix/tests/plot/test_stereographic_plot.py
@@ -102,10 +102,15 @@ def test_grids(self):
         assert ax._azimuth_resolution == azimuth_res
         assert ax._polar_resolution == polar_res
 
-        ax.stereographic_grid(azimuth_resolution=30, polar_resolution=45)
+        alpha = 0.5
+        with plt.rc_context({"grid.alpha": alpha}):
+            ax.stereographic_grid(azimuth_resolution=30, polar_resolution=45)
         assert ax._azimuth_resolution == 30
         assert ax._polar_resolution == 45
 
+        assert len(ax.collections) == 2
+        assert all([coll.get_alpha() for coll in ax.collections])
+
         plt.close("all")
 
     def test_set_labels(self):
diff --git a/orix/tests/test_vector3d.py b/orix/tests/test_vector3d.py
index 5ce453b3..90fc1290 100644
--- a/orix/tests/test_vector3d.py
+++ b/orix/tests/test_vector3d.py
@@ -642,6 +642,24 @@ def test_get_circle(self):
         assert np.allclose(c.mean().data, [0, 0, 0], atol=1e-2)
         assert np.allclose(v.cross(c[0, 0]).data, [1, 0, 0])
 
+    def test_from_path_ends(self):
+        vx = Vector3d.xvector()
+        vy = Vector3d.yvector()
+
+        v_xy = Vector3d.from_path_ends(Vector3d.stack((vx, vy)))
+        assert v_xy.size == 27
+        assert np.allclose(v_xy.polar, np.pi / 2)
+        assert np.allclose(v_xy[-1].data, vy.data)
+
+        v_xyz = Vector3d.from_path_ends(
+            [[1, 0, 0], [0, 1, 0], [0, 0, 1]], steps=150, close=True
+        )
+        assert v_xyz.size == 115
+        assert np.allclose(v_xyz[-1].data, vx.data)
+
+        with pytest.raises(ValueError, match="No vectors are perpendicular"):
+            _ = Vector3d.from_path_ends(Vector3d.stack((vx, -vx)))
+
 
 class TestPlotting:
     v = Vector3d(
diff --git a/orix/vector/vector3d.py b/orix/vector/vector3d.py
index c2a03ab4..5fcda85b 100644
--- a/orix/vector/vector3d.py
+++ b/orix/vector/vector3d.py
@@ -476,6 +476,76 @@ def zvector(cls) -> Vector3d:
         """Return a unit vector in the z-direction."""
         return cls((0, 0, 1))
 
+    @classmethod
+    def from_path_ends(
+        cls,
+        vectors: Union[list, tuple, Vector3d],
+        close: bool = False,
+        steps: int = 100,
+    ) -> Vector3d:
+        r"""Return vectors along the shortest path on the sphere between
+        two or more consectutive vectors.
+
+        Parameters
+        ----------
+        vectors
+            Two or more vectors to get paths between.
+        close
+            Whether to append the first to the end of ``vectors`` in
+            order to close the paths when more than two vectors are
+            passed. Default is False.
+        steps
+            Number of vectors in the great circle about the normal
+            vector between each two vectors *before* restricting the
+            circle to the path between the two. Default is 100. More
+            steps give a smoother path on the sphere.
+
+        Returns
+        -------
+        paths
+            Vectors along the shortest path(s) between given vectors.
+
+        Notes
+        -----
+        The vectors along the shortest path on the sphere between two
+        vectors :math:`v_1` and :math:`v_2` are found by first getting
+        the vectors :math:`v_i` along the great circle about the vector
+        normal to these two vectors, and then only keeping the part of
+        the circle between the two vectors. Vectors within this part
+        satisfy these two conditions
+
+        .. math::
+            (v_1 \times v_i) \cdot (v_1 \times v_2) \geq 0,
+            (v_2 \times v_i) \cdot (v_2 \times v_1) \geq 0.
+        """
+        v = Vector3d(vectors).flatten()
+
+        if close:
+            v = Vector3d(np.concatenate((v.data, v[0].data)))
+
+        paths_list = []
+
+        n = v.size - 1
+        for i in range(n):
+            v1, v2 = v[i : i + 2]
+            v_normal = v1.cross(v2)
+            v_circle = v_normal.get_circle(steps=steps)
+
+            cond1 = v1.cross(v_circle).dot(v1.cross(v2)) >= 0
+            cond2 = v2.cross(v_circle).dot(v2.cross(v1)) >= 0
+
+            v_path = v_circle[cond1 & cond2]
+
+            to_concatenate = (v1.data, v_path.data)
+            if i == n - 1:
+                to_concatenate += (v2.data,)
+            paths_list.append(np.concatenate(to_concatenate, axis=0))
+
+        paths_data = np.concatenate(paths_list, axis=0)
+        paths = Vector3d(paths_data)
+
+        return paths
+
     def angle_with(self, other: Vector3d, degrees: bool = False) -> np.ndarray:
         """Return the angles between these vectors in other vectors.