[BUG] transform() accepts complex input but gives incorrect results

[jmtayloruk]

If the function passed to transform() returns complex results, the Hankel transform is not computed correctly (and warnings are emitted from numpy and scipy). I think this is a regression - I am almost certain this used to work in previous versions of the hankel library.

I cannot find any explicit mention in the docs or changlog about complex functions, but this is certainly a common use-case in optical diffraction.

It looks like there is some intention to support complex functions (test for np.isrealobj(k) in the source code). Although not documented anywhere, this line implies that passing in a complex-typed k could signal that the function f() will return a complex type. That is rather odd for an undocumented behaviour (why shouldn't I have real k but complex f()?). However, even if I do that it still does not quite work (which may in part be because quad() does not accept complex functions).

There is a (rather clumsy) external workaround - see below. While a more efficient fix is surely available, one simple option would be to build this workaround in to transform().

Code To Reproduce

import numpy as np
import matplotlib.pyplot as plt
import warnings
import hankel
print(hankel.__version__)

(N,h) = (1e4, 1e-7)
k = np.arange(0, 30, 0.5)

def PupilFuncReal(r):
    r = np.array(r)
    result = np.zeros(r.shape)
    result[r<1.0] = 1.0
    return result

# Perform Hankel transform over a real function - SUCCESS
warnings.resetwarnings()
hankelTransform = hankel.HankelTransform(nu=0, N=N, h=h)
hhat = hankelTransform.transform(PupilFuncReal, k, ret_err=False) / (2*np.pi)
plt.plot(hhat); plt.show()

def PupilFuncComplex(r):
    r = np.array(r)
    result = np.zeros(r.shape, dtype='complex128')
    result[r<1.0] = 1.0j
    return result

# Perform Hankel transform over a function returning a complex result - FAIL
# Emits warnings:
# [...]anaconda/lib/python3.5/site-packages/hankel/hankel.py:170: ComplexWarning: Casting complex values to real discards the imaginary part
# [...]anaconda/lib/python3.5/site-packages/scipy/integrate/quadpack.py:455: ComplexWarning: Casting complex values to real discards the imaginary part
warnings.resetwarnings()
hhat = hankelTransform.transform(PupilFuncComplex, k, ret_err=False) / (2*np.pi)
print(hhat.dtype)
plt.plot(hhat); plt.show()

# Perform Hankel transform over a function returning a complex result,
# with k cast to complex (since hankel.py:168 transform() tests for this case)
# ALMOST SUCCEEDS. Value at k=0 is incorrect, and still emits concerning warning.
# Emits warnings:
# [...]anaconda/lib/python3.5/site-packages/scipy/integrate/quadpack.py:455: ComplexWarning: Casting complex values to real discards the imaginary part
warnings.resetwarnings()
hhat = hankelTransform.transform(PupilFuncComplex, k.astype('complex128'), ret_err=False) / (2*np.pi)
plt.plot(hhat); plt.show()

Workaround

# Perform Hankel transform over a function returning a complex result,
# by splitting into a real and an imaginary integral. SUCCESS
warnings.resetwarnings()
hhat_r = hankelTransform.transform(lambda r: PupilFuncComplex(r).real, k, ret_err=False) / (2*np.pi)
hhat_i = hankelTransform.transform(lambda r: PupilFuncComplex(r).imag, k, ret_err=False) / (2*np.pi)
hhat = hhat_r + 1j*hhat_i
print(hhat.dtype)
plt.plot(hhat.real)
plt.plot(hhat.imag)
plt.show()

[steven-murray]

Ah! Good catch. I will implement a fix ASAP.

[steven-murray]

@jmtayloruk I've made a PR for this (noted above). Would love if you could give it a run and see if there are any remaining issues!