diff --git a/quaternionic/algebra.py b/quaternionic/algebra.py
index 7bb075d..c4158cc 100644
--- a/quaternionic/algebra.py
+++ b/quaternionic/algebra.py
@@ -259,6 +259,25 @@ def log(q, qout):
         qout[3] = f * q[3]
 
 
+@attach_typelist_and_signature([(float64[:], float64[:])], '(n)->()')
+def angle(q, qout):
+    """Return angle (in radians) through which input quaternion rotates a vector
+
+    If `q = np.exp(v̂ * θ/2)` for some unit vector `v̂` and an angle `θ` ∈[-2π,2π],
+    then this function returns `abs(θ)`.  This equals 2*abs(log(q)), but is more
+    efficient.
+
+    """
+    b = np.sqrt(q[1]**2 + q[2]**2 + q[3]**2)
+    if b <= _quaternion_resolution * np.abs(q[0]):
+        if q[0] < 0.0:
+            qout[0] = 2*np.pi
+        else:
+            qout[0] = 0.0
+    else:
+        qout[0] = 2*np.abs(np.arctan2(b, q[0]))
+
+
 @attach_typelist_and_signature([(float64[:], float64[:])], '(n)->(n)')
 def sqrt(q, qout):
     """Return square-root of input quaternion √q.
diff --git a/quaternionic/arrays.py b/quaternionic/arrays.py
index a1821c6..5318d89 100644
--- a/quaternionic/arrays.py
+++ b/quaternionic/arrays.py
@@ -105,10 +105,18 @@ def __array_finalize__(self, obj):
                 )
 
         def __array_function__(self, func, types, args, kwargs):
-            output = super().__array_function__(func, types, args, kwargs)
-            if func in [np.ones_like]:
-                # Want the last dimension to equal [1,0,0,0] not [1,1,1,1]
-                output.vector = 0
+            from . import algebra_ufuncs as algebra
+
+            if func == np.angle:
+                output = np.zeros(args[0].shape[:-1], dtype=dtype)
+                algebra.angle(args[0].ndarray, output)
+                if kwargs.get("deg", False):
+                    output *= 180/np.pi
+            else:
+                output = super().__array_function__(func, types, args, kwargs)
+                if func in [np.ones_like]:
+                    # Want the last dimension to equal [1,0,0,0] not [1,1,1,1]
+                    output.vector = 0
             return output
 
         def __array_ufunc__(self, ufunc, method, *args, **kwargs):
diff --git a/tests/test_algebra.py b/tests/test_algebra.py
index e56c0e2..0ef5f34 100644
--- a/tests/test_algebra.py
+++ b/tests/test_algebra.py
@@ -38,6 +38,27 @@ def test_basis_multiplication():
     assert i*j*k == -one
 
 
+def test_array_function(array):
+    np.random.seed(1234)
+
+    # Test np.angle
+    assert np.angle(array(1, 0, 0, 0)) == 0.0
+    assert np.angle(array(-1, 0, 0, 0)) == 2*np.pi
+    assert np.angle(array(0, 1, 0, 0)) == np.pi
+    assert np.angle(array(0, -1, 0, 0)) == np.pi
+    assert np.angle(array(0, 0, 1, 0, 0)) == np.pi
+    assert np.angle(array(0, 0, -1, 0, 0)) == np.pi
+    assert np.angle(array(0, 0, 0, 1)) == np.pi
+    assert np.angle(array(0, 0, 0, -1)) == np.pi
+    angle_precision = 2e-15
+    for θ in np.random.uniform(-2*np.pi, 2*np.pi, size=100):
+        for _ in range(20):
+            v̂ = array.from_vector_part(np.random.normal(size=3)).normalized
+            R = np.exp(v̂ * θ/2)
+            assert abs(abs(θ) - np.angle(R)) < angle_precision
+            assert abs(abs(np.rad2deg(θ)) - np.angle(R, deg=True)) < np.rad2deg(angle_precision)
+
+
 def test_array_ufunc(array):
     np.random.seed(1234)