Revised polynomials.py and created test_polynomial.py

Author: mirjagranfors

deeptrack/backend/polynomials.py

@@ -1,12 +1,86 @@
-""" Expands the set of polynomials available through scipy
+"""Bessel and Riccati-Bessel polynomials.
+
+This file defines a set of functions for computing Bessel and Riccati-Bessel 
+polynomials and their derivatives. It expands the capabilities of `scipy`.
+
+Functions
+---------
+besselj(
+    l: Union[int, float],
+    x: Union[int, float, np.ndarray]
+) -> Union[float, np.ndarray]
+    Computes the Bessel polynomial of the first kind for a given order `l`
+    and input `x`.
+
+dbesselj(
+    l: Union[int, float],
+    x: Union[int, float, np.ndarray]
+) -> Union[float, np.ndarray]
+    Computes the first derivative of the Bessel polynomial of the first kind.
+
+bessely(
+    l: Union[int, float], 
+    x: Union[int, float, np.ndarray]
+) -> Union[float, np.ndarray]
+    Computes the Bessel polynomial of the second kind for a given order `l`
+    and input `x`.
+
+dbessely(
+    l: Union[int, float], 
+    x: Union[int, float, np.ndarray]
+) -> Union[float, np.ndarray]
+    Computes the first derivative of the Bessel polynomial of the second kind.
+
+ricbesj(
+    l: Union[int, float], 
+    x: Union[int, float, np.ndarray]
+) -> Union[float, np.ndarray]
+    Computes the Riccati-Bessel polynomial of the first kind.
+
+dricbesj(
+    l: Union[int, float], 
+    x: Union[int, float, np.ndarray]
+) -> Union[float, np.ndarray]
+    Computes the first derivative of the Riccati-Bessel polynomial of the
+    first kind.
+
+ricbesy(
+    l: Union[int, float], 
+    x: Union[int, float, np.ndarray]
+) -> Union[float, np.ndarray]
+    Computes the Riccati-Bessel polynomial of the second kind.
+
+dricbesy(
+    l: Union[int, float], 
+    x: Union[int, float, np.ndarray]
+) -> Union[float, np.ndarray]
+    Computes the first derivative of the Riccati-Bessel polynomial of the
+    second kind.
+
+ricbesh(
+    l: Union[int, float], 
+    x: Union[int, float, np.ndarray]
+) -> Union[float, np.ndarray]
+    Computes the Riccati-Bessel polynomial of the third kind.
+
+dricbesh(
+    l: Union[int, float], 
+    x: Union[int, float, np.ndarray]
+) -> Union[float, np.ndarray]
+    Computes the first derivative of the Riccati-Bessel polynomial of the
+    third kind.
 """
 
 import numpy as np
-from scipy.special import jv, spherical_jn, h1vp, jvp, yv
+from scipy.special import jv, h1vp, yv
+from typing import Union
 
 
-def ricbesj(l, x):
-    """The Riccati-Bessel polynomial of the first kind.
+def besselj(
+        l: Union[int, float],
+        x: Union[int, float, np.ndarray],
+    ) -> Union[float, np.ndarray]:
+    """The Bessel polynomial of the first kind.
 
     Parameters
     ----------
@@ -21,11 +95,14 @@ def ricbesj(l, x):
         The polynomial evaluated at x
     """
 
-    return np.sqrt(np.pi * x / 2) * besselj(l + 0.5, x)
+    return jv(l, x)
 
 
-def dricbesj(l, x):
-    """The first derivative of the Riccati-Bessel polynomial of the first kind.
+def dbesselj(
+        l: Union[int, float],
+        x: Union[int, float, np.ndarray],
+    ) -> Union[float, np.ndarray]:
+    """The first derivative of the Bessel polynomial of the first kind.
 
     Parameters
     ----------
@@ -40,13 +117,14 @@ def dricbesj(l, x):
         The polynomial evaluated at x
     """
 
-    return 0.5 * np.sqrt(np.pi / x / 2) * besselj(l + 0.5, x) + np.sqrt(
-        np.pi * x / 2
-    ) * dbesselj(l + 0.5, x)
+    return 0.5 * (besselj(l - 1, x) - besselj(l + 1, x))
 
 
-def ricbesy(l, x):
-    """The Riccati-Bessel polynomial of the second kind.
+def bessely(
+        l: Union[int, float],
+        x: Union[int, float, np.ndarray],
+    ) -> Union[float, np.ndarray]:
+    """The Bessel polynomial of the second kind.
 
     Parameters
     ----------
@@ -61,11 +139,15 @@ def ricbesy(l, x):
         The polynomial evaluated at x
     """
 
-    return -np.sqrt(np.pi * x / 2) * bessely(l + 0.5, x)
+    return yv(l, x)
 
 
-def dricbesy(l, x):
-    """The first derivative of the Riccati-Bessel polynomial of the second kind.
+def dbessely(
+        l: Union[int, float],
+        x: Union[int, float, np.ndarray],
+    ) -> Union[float, np.ndarray]:
+    
+    """The first derivative of the Bessel polynomial of the second kind.
 
     Parameters
     ----------
@@ -79,14 +161,15 @@ def dricbesy(l, x):
     float, ndarray
         The polynomial evaluated at x
     """
-
-    return -0.5 * np.sqrt(np.pi / 2 / x) * yv(l + 0.5, x) - np.sqrt(
-        np.pi * x / 2
-    ) * dbessely(l + 0.5, x)
+    
+    return 0.5 * (bessely(l - 1, x) - bessely(l + 1, x))
 
 
-def ricbesh(l, x):
-    """The Riccati-Bessel polynomial of the third kind.
+def ricbesj(
+        l: Union[int, float],
+        x: Union[int, float, np.ndarray],
+    ) -> Union[float, np.ndarray]:
+    """The Riccati-Bessel polynomial of the first kind.
 
     Parameters
     ----------
@@ -100,11 +183,16 @@ def ricbesh(l, x):
     float, ndarray
         The polynomial evaluated at x
     """
-    return np.sqrt(np.pi * x / 2) * h1vp(l + 0.5, x, False)
+
+    return np.sqrt(np.pi * x / 2) * besselj(l + 0.5, x)
 
 
-def dricbesh(l, x):
-    """The first derivative of the Riccati-Bessel polynomial of the third kind.
+def dricbesj(
+        l: Union[int, float],
+        x: Union[int, float, np.ndarray],
+    ) -> Union[float, np.ndarray]:
+    """The first derivative of the Riccati-Bessel polynomial of the first
+    kind.
 
     Parameters
     ----------
@@ -118,14 +206,17 @@ def dricbesh(l, x):
     float, ndarray
         The polynomial evaluated at x
     """
-    xi = 0.5 * np.sqrt(np.pi / 2 / x) * h1vp(l + 0.5, x, False) + np.sqrt(
+
+    return 0.5 * np.sqrt(np.pi / x / 2) * besselj(l + 0.5, x) + np.sqrt(
         np.pi * x / 2
-    ) * h1vp(l + 0.5, x, True)
-    return xi
+    ) * dbesselj(l + 0.5, x)
 
 
-def besselj(l, x):
-    """The Bessel polynomial of the first kind.
+def ricbesy(
+        l: Union[int, float],
+        x: Union[int, float, np.ndarray],
+    ) -> Union[float, np.ndarray]:
+    """The Riccati-Bessel polynomial of the second kind.
 
     Parameters
     ----------
@@ -140,11 +231,15 @@ def besselj(l, x):
         The polynomial evaluated at x
     """
 
-    return jv(l, x)
+    return -np.sqrt(np.pi * x / 2) * bessely(l + 0.5, x)
 
 
-def dbesselj(l, x):
-    """The first derivative of the Bessel polynomial of the first kind.
+def dricbesy(
+        l: Union[int, float],
+        x: Union[int, float, np.ndarray],
+    ) -> Union[float, np.ndarray]:
+    """The first derivative of the Riccati-Bessel polynomial of the second
+    kind.
 
     Parameters
     ----------
@@ -158,11 +253,17 @@ def dbesselj(l, x):
     float, ndarray
         The polynomial evaluated at x
     """
-    return 0.5 * (besselj(l - 1, x) - besselj(l + 1, x))
+
+    return -0.5 * np.sqrt(np.pi / 2 / x) * yv(l + 0.5, x) - np.sqrt(
+        np.pi * x / 2
+    ) * dbessely(l + 0.5, x)
 
 
-def bessely(l, x):
-    """The Bessel polynomial of the second kind.
+def ricbesh(
+        l: Union[int, float],
+        x: Union[int, float, np.ndarray],
+    ) -> Union[float, np.ndarray]:
+    """The Riccati-Bessel polynomial of the third kind.
 
     Parameters
     ----------
@@ -177,11 +278,14 @@ def bessely(l, x):
         The polynomial evaluated at x
     """
 
-    return yv(l, x)
+    return np.sqrt(np.pi * x / 2) * h1vp(l + 0.5, x, False)
 
 
-def dbessely(l, x):
-    """The first derivative of the Bessel polynomial of the second kind.
+def dricbesh(
+        l: Union[int, float],
+        x: Union[int, float, np.ndarray],
+    ) -> Union[float, np.ndarray]:
+    """The first derivative of the Riccati-Bessel polynomial of the third kind.
 
     Parameters
     ----------
@@ -195,4 +299,8 @@ def dbessely(l, x):
     float, ndarray
         The polynomial evaluated at x
     """
-    return 0.5 * (bessely(l - 1, x) - bessely(l + 1, x))
+
+    xi = 0.5 * np.sqrt(np.pi / 2 / x) * h1vp(l + 0.5, x, False) + np.sqrt(
+        np.pi * x / 2
+    ) * h1vp(l + 0.5, x, True)
+    return xi

deeptrack/test/backend/__init__.py

@@ -1 +1,2 @@
-from .test_core import *
+from .test_core import *
+from .test_polynomials import *

deeptrack/test/backend/test_polynomials.py

@@ -0,0 +1,72 @@
+import unittest
+
+from deeptrack.backend import polynomials
+
+
+class TestPolynomials(unittest.TestCase):
+
+    def test_Besselj(self):
+        Besselj = polynomials.besselj
+
+        self.assertEqual(Besselj(0, 0), 1)
+        self.assertEqual(Besselj(1, 0), 0)
+
+        self.assertEqual(Besselj(1, 5), -Besselj(-1, 5))
+        self.assertTrue(Besselj(1, 5) < 0)
+
+
+    def test_dBesselj(self):
+        dBesselj = polynomials.dbesselj
+
+        self.assertEqual(dBesselj(1,0), 0.5)
+
+
+    def test_Bessely(self):
+        Bessely = polynomials.bessely
+
+        self.assertEqual(Bessely(1, 3), -Bessely(-1,3))
+        self.assertTrue(Bessely(1, 3) > 0)
+
+
+    def test_dBessely(self):
+        dBessely = polynomials.dbessely
+
+        self.assertEqual(dBessely(1, 3), -dBessely(-1,3))
+        self.assertTrue(dBessely(1, 3) > 0)
+
+
+    def test_RicBesj(self):
+        Ricbesj = polynomials.ricbesj
+
+        self.assertEqual(Ricbesj(4, 0), 0)
+        self.assertTrue(abs(Ricbesj(1, 8)- 0.2691) < 1e-3)
+
+
+    def test_dRicBesj(self):
+        dRicBesj = polynomials.dricbesj
+
+        self.assertTrue(abs(dRicBesj(1, 1) - 0.5403) < 1e-3)
+
+
+    def test_RicBesy(self):
+        RicBesy = polynomials.ricbesy
+
+        self.assertTrue(abs(RicBesy(2, 3) - 0.8011) < 1e-3)
+
+
+    def test_dRicBesy(self):
+        dRicBesy = polynomials.dricbesy
+
+        self.assertTrue(abs(dRicBesy(1, 1) + 0.8414) < 1e-3)
+
+        
+    def test_RicBesh(self):
+        RicBesh = polynomials.ricbesh
+
+        self.assertTrue(abs(RicBesh(3, 2) - (0.1214 - 2.968j)) < 1e-3)
+
+
+    def test_dRicBesh(self):
+        dRicBesh = polynomials.dricbesh
+
+        self.assertTrue(abs(dRicBesh(2, 6) - (-0.9321-0.2206j)) < 1e-3)