Merge c40925b into 166b04e

.gitignore

@@ -10,7 +10,7 @@ docs/demos/.ipynb_checkpoints/*
 .tox/*
 \.coverage
 .eggs/
-
+hankel/_version.py
 htmlcov/
 
 *.aux

CHANGELOG.rst

@@ -1,6 +1,12 @@
 Changelog
 =========
 
+v1.1.0
+------
+**Enhancements**
+
+- Add ability to perform hankel transforms on complex functions.
+
 v1.0.0
 ------
 **Enhancements**

hankel/__init__.py

@@ -27,7 +27,7 @@
 
 try:
     from hankel._version import __version__
-except ModuleNotFoundError:  # pragma: nocover
+except ImportError:  # pragma: nocover
     # package is not installed
     pass
 

hankel/hankel.py

@@ -107,7 +107,7 @@ def _norm(self, inverse=False):
         return 1
 
     def _get_series(self, f, k=1):
-        with np.errstate(divide="ignore"):  # numpy safely divids by 0
+        with np.errstate(divide="ignore"):  # numpy safely divides by 0
             args = np.divide.outer(self.x, k).T  # x = r*k
         return self._series_fac * f(args) * safe_power(self.x, self._r_power)
 
@@ -164,11 +164,16 @@ def transform(self, f, k=1, ret_err=True, ret_cumsum=False, inverse=False):
         # The basic transform has a norm of 1.
         # But when doing FT's, this depends on the dimensionality.
         norm = self._norm(inverse)
+
         # calculate the result for non zero k (int k -> real ret)
-        ret = np.empty_like(k, float) if np.isrealobj(k) else np.empty_like(k)
         summation = self._get_series(f, k_tmp)
+        # return value should be same dtype as summation (allows complex)
+        ret = np.empty(k.shape, dtype=summation.dtype)
+
         ret[kn0] = np.array(norm * np.sum(summation, axis=-1) / knorm)
 
+        is_cmplx = not np.isrealobj(summation)
+
         # care about k=0
         ret_0 = 0
         err_0 = 0
@@ -182,11 +187,23 @@ def transform(self, f, k=1, ret_err=True, ret_cumsum=False, inverse=False):
                 int_fac = j_lim(self.nu) * norm
 
                 def integrand(r):
-                    return f(r) * safe_power(r, lim_r_pow)
+                    return f(r).real * safe_power(r, lim_r_pow)
 
                 int_res = quad(integrand, 0, np.inf)
-                ret_0 = int_res[0] * int_fac
-                err_0 = int_res[1] * int_fac
+
+                if is_cmplx:
+                    # For a complex function, need to do the complex part of the integral.
+                    def integrand(r):
+                        return f(r).imag * safe_power(r, lim_r_pow)
+
+                    int_res_cmplx = quad(integrand, 0, np.inf)
+                    ret_0 = (int_res[0] + 1j * int_res_cmplx[0]) * int_fac
+                    err_0 = (int_res[1] + 1j * int_res_cmplx[1]) * int_fac
+
+                else:
+                    ret_0 = int_res[0] * int_fac
+                    err_0 = int_res[1] * int_fac
+
             elif self.nu < nu_th:
                 ret_0 = np.nan
             ret[k_0] = ret_0

tests/test_transforms.py

@@ -137,3 +137,24 @@ def test_k_zero(nu, alt):
         assert np.isclose(ans, 0, rtol=1e-3)
     else:
         assert np.isclose(ans, 1, rtol=1e-3)
+
+
+@pytest.mark.parametrize(
+    "s, nu, k, N, h", [[-2, 1, 0.01, 300, 10 ** -3.2], [-2, 1, 1, 300, 10 ** -3.2],],
+)
+def test_complex(s, nu, k, N, h):
+    """
+    Test f(r) = 1/r^2 j, nu=0
+    """
+
+    ht = HankelTransform(nu=nu, N=N, h=h)
+    ans = ht.transform(lambda x: x ** s * 1j, k, False, False)
+    if nu - s <= 0 and (nu - s) % 2 == 0:
+        raise Exception("Can't have a negative integer for gamma")
+
+    anl = (
+        2 ** (s + 1) * gamma(0.5 * (2 + nu + s)) / k ** (s + 2) / gamma(0.5 * (nu - s))
+    ) * 1j
+
+    print("Numerical Result: ", ans, " (required %s)" % anl)
+    assert np.isclose(ans, anl, rtol=1e-3)