diff --git a/src/aspire/basis/basis_utils.py b/src/aspire/basis/basis_utils.py
index 9f1bbf977a..4a7fcfe98d 100644
--- a/src/aspire/basis/basis_utils.py
+++ b/src/aspire/basis/basis_utils.py
@@ -8,7 +8,7 @@
 import numpy as np
 from numpy import diff, exp, log, pi
 from numpy.polynomial.legendre import leggauss
-from scipy.special import jn, jv, lpmv
+from scipy.special import jn, jv, lpmv, sph_harm
 
 from aspire.utils import ensure
 from aspire.utils.coor_trans import grid_2d, grid_3d
@@ -126,26 +126,14 @@ def real_sph_harmonic(j, m, theta, phi):
     :return: The real spherical harmonics evaluated at the points (theta, phi).
     """
     abs_m = abs(m)
-    y = lpmv(abs_m, j, np.cos(theta))
-
-    # Beware of using just np.prod in the denominator here
-    # Unless we use float64, values in the denominator > 13! will be incorrect
-    try:
-        y = (
-            np.sqrt(
-                (2 * j + 1)
-                / (4 * pi)
-                / np.prod(range(j - abs_m + 1, j + abs_m + 1), dtype=np.float64)
-            )
-            * y
-        )
-    except RuntimeWarning:
-        logger.error("debug")
 
+    y = sph_harm(abs_m, j, phi, theta)
     if m < 0:
-        y = np.sqrt(2) * np.sin(abs_m * phi) * y
+        y = np.sqrt(2) * np.imag(y)
     elif m > 0:
-        y = np.sqrt(2) * np.cos(abs_m * phi) * y
+        y = np.sqrt(2) * np.real(y)
+    else:
+        y = np.real(y)
 
     return y
 
diff --git a/tutorials/examples/cov3d_simulation.py b/tutorials/examples/cov3d_simulation.py
index 04da70ef68..4089af3703 100644
--- a/tutorials/examples/cov3d_simulation.py
+++ b/tutorials/examples/cov3d_simulation.py
@@ -28,6 +28,7 @@
 
 # Create a simulation object with specified filters
 sim = Simulation(
+    L=img_size,
     n=num_imgs,
     C=num_vols,
     unique_filters=[RadialCTFFilter(defocus=d) for d in np.linspace(1.5e4, 2.5e4, 7)],