diff --git a/math.jl b/math.jl
index d407979..b51a39f 100644
--- a/math.jl
+++ b/math.jl
@@ -29,7 +29,14 @@ function integral(np, lp, n, l, μω)
 end
 
 "Gives <n' l'|| p ||n l> for the HO basis"
-reduced_matrix_element(np, lp, n, l, μω) = (-1)^lp * sqrt(2*lp + 1) * sqrt(2*l + 1) * wigner3j(Float64, lp, 1, l, 0, 0, 0) * integral(np, lp, n, l, μω)
+function reduced_matrix_element(np, lp, n, l, μω)::ComplexF64
+    wig::Float64 = wigner3j(Float64, lp, 1, l, 0, 0, 0)
+    if wig == 0
+        return 0
+    else
+        return (-1)^lp * sqrt(2*lp + 1) * sqrt(2*l + 1) * wig * integral(np, lp, n, l, μω)
+    end
+end
 
 "Matrix element <n1p l1p n2p l2p| p1⋅p2 |n1 l1 n2 l2> (Ref: de-Shalit & Talmi, Eq 15.5)"
 function racahs_reduction_formula(n1p, l1p, n2p, l2p, n1, l1, n2, l2, Λ, μ1ω1, μ2ω2)