Deprecate scipy Legendre function

spharpy/special.py

@@ -375,7 +375,7 @@ def legendre_function(n, m, z, cs_phase=True):
     z : ndarray, double
         The argument as an array
     cs_phase : bool, optional
-        Whether to use include the Condon-Shotley phase term (-1)^m or not
+        Whether to include the Condon-Shotley phase term (-1)^m or not
 
     Returns
     -------
@@ -392,16 +392,15 @@ def legendre_function(n, m, z, cs_phase=True):
     """
     z = np.atleast_1d(z)
 
-    if np.abs(m) > n:
-        legendre = np.zeros(z.shape)
-    else:
-        legendre = np.zeros(z.shape)
-        for idx, arg in zip(count(), z):
-            leg, _ = _spspecial.lpmn(m, n, arg)
-            if np.mod(m, 2) != 0 and not cs_phase:
-                legendre[idx] = -leg[-1, -1]
-            else:
-                legendre[idx] = leg[-1, -1]
+    # squeeze required because the legendre function introduced in scipy 1.15
+    # returns a 2D array, whereas the previous function `lpmn` returned a 1D
+    # array
+    legendre = np.squeeze(_spspecial.assoc_legendre_p(n, m, z))
+
+    # remove Condon-Shortley phase
+    if not cs_phase and m % 2:
+        legendre *= -1
+
     return legendre
 
 

tests/test_special.py

@@ -227,3 +227,39 @@ def test_spherical_harmonic_gradient_phi_real():
         desired = desired_all[:, acn]
 
         npt.assert_allclose(actual, desired, rtol=1e-10, atol=1e-10)
+
+
+@pytest.mark.parametrize(['condon_shortley'], [(True, ), (False, )])
+@pytest.mark.parametrize(['m'], [(-1, ), (0, ), (1, )])
+def test_legendre(condon_shortley, m):
+    """
+    Test values of the Legendre functions for first order and all degrees.
+    """
+    z = np.linspace(-1, 1, 11)
+
+    # Manually computed desired values according to
+    # Rafaely (2019), Fundamentals of Spherical Array Processing, Table 1.3
+    if m == -1:
+        desired = .5 * np.sqrt(1 - z**2)
+    if m == 0:
+        desired = z.copy()
+    if m == 1:
+        desired =  -np.sqrt(1 - z**2)
+
+    # remove Condon-Shortley phase
+    if not condon_shortley and m % 2:
+        desired *= -1
+
+    # compute and compare actual values
+    actual = special.legendre_function(1, m, z, condon_shortley)
+    npt.assert_almost_equal(actual, desired, 10)
+
+
+@pytest.mark.parametrize(['m'], [(-2, ), (2, )])
+def test_legendre_degree_out_of_range(m):
+    """Test if zero is returned if the degree m is larger than the order n"""
+    n = 1
+    z = np.linspace(-1, 1, 11)
+
+    npt.assert_equal(special.legendre_function(n, m, z),
+                     np.zeros_like(z))