hankel: A Python library for performing simple and accurate Hankel transformations

TL;DR

The paper introduces 'hankel', a Python library that simplifies and accelerates the computation of Hankel transforms, which are vital in physical sciences for radial Fourier analysis in multiple dimensions.

Contribution

It provides a pure Python implementation of Hankel transforms, emphasizing simplicity, efficiency, and usability for scientific computations.

Findings

Efficient computation of Hankel transforms in Python.

Applicable to radial Fourier transforms in physics.

Extensive documentation and examples provided.

Abstract

This paper presents \textsc{hankel}, a pure-python code for solving Hankel-type integrals and transforms. Such transforms are common in the physical sciences, especially appearing as the radial solution to angularly symmetric Fourier Transforms in arbitrary dimensions. The code harnesses the advantages of solving such transforms via the one-dimensional Hankel transform -- an increase in conceptual simplicity and efficiency -- and implements them in the user-friendly and flexible Python language. We discuss several limitations of the adopted method, and point to the code's extensive documentation for further examples.

Full text

hankel: A Python library for performing simple and accurate Hankel

transformations

Steven G. Murray

International Centre for Radio Astronomy Research (ICRAR), Curtin University, Bentley, WA 6102, Australia

ARC Centre of Excellence for All-Sky Astrophysics in 3 Dimensions (ASTRO 3D)

School of Earth and Space Exploration, Arizona State University, Tempe, AZ, 85281, USA

Francis J. Poulin

Department of Applied Mathematics, University of Waterloo

††margin:

DOI: 10.21105/joss.01397

Software

•

Review

•

Repository

•

Archive

**Submitted: 02 April 2019

Published: 31 May 2019**

**License

Authors of papers retain copyright and release the work under a Creative Commons Attribution 4.0 International License (CC-BY).**

Summary

The Hankel transform is a one-dimensional functional transform involving a Bessel-function kernel. Moreover, it is the radial solution to an angularly symmetric Fourier transform of any dimension, rendering it very useful in several fields of application. The NASA Astronomical Data Service yields over 700 refereed articles including the term “hankel transform”, in fields as diverse as astronomy, geophysics, fluid mechanics, electrodynamics, thermodynamics and acoustics.

As an example, in cosmology, the density field of the Universe is expected to be isotropic. One of the primary means of describing this field is via its Fourier power spectrum, or equivalently its spatial autocorrelation function. Due to the isotropy of the field, these can be related by an angularly symmetric Fourier transform, which is more simply expressed as a Hankel transform (Szapudi et al. 2005).

Another example that arises in both geophysical and astrophysical contexts is in regards to vortices. The radially-symmetric vortical solution to Laplace’s equation in two, three or even higher dimensions can be performed quickly and accurately via the Hankel transform (Carton 2001).

Conceptually, computation of such problems using the Hankel transform, in contrast to the Fourier transform, has the advantage of reducing the problem’s dimensionality to unity, regardless of the original dimensionality. Analytically, this may be a useful tool in solving the transform. Numerically, it naively promises to enhance efficiency.

Despite these advantages, the Hankel transform introduces some numerical challenges. Most importantly, the Hankel transform is a highly oscillatory integral, especially for large values of the transformation variable, k (henceforth we will use r to denote the magnitude of the real-space co-ordinate). Highly oscillatory integrals are a topic of much interest in applied mathematics, and there does not exist a general optimal solution to numerically evaluate them in general (Huybrechs and Olver 2009). Nevertheless, Ogata (2005) determined that a double-exponential variable transformation based on the zeros of the Bessel function (Ooura and Mori 1999) has the property that the numerical integral converges with many fewer divisions compared to naively computing the transform integral. This procedure is able to efficiently and accurately evaluate the Hankel integral (and hence the Hankel transform) in many cases.

The purpose of hankel is to provide a dead-simple intuitive pure-Python interface for performing Hankel integrals and transforms, written with both Python 2 and 3 compatibility. It enables the accurate determination of the transform of a callable function at any arbitrary value of the transform variable, k, and utilises the proven method of Ogata (2005) to do this efficiently. In addition, it recognizes that the most common application of the Hankel transform is in the context of the radial Fourier transform, and it provides an additional interface dedicated to making the connection between these transforms more transparent.

The chief performance-critical components of hankel are the evaluation of the zeros of the Bessel function, and the sum of terms required for integration. The former is made efficient by a tiered approach – using the efficient scipy for integer-order Bessel functions, directly returning regular arrays for order-1/2, and using the very accurate mpmath for all other orders. The latter is made efficient by utilising numpy, bringing it close to C-level performance.

The hankel package is thoroughly tested to ensure accuracy of transforms, by comparing to known analytic solutions. These tests, supported by continuous integration, are also useful for the user who wishes to explore the numerical limitations of the method. Aside from functions which are theoretically divergent, the method can struggle to transform several classes of functions, including those with very sharp features, especially at small r. The method itself has two free parameters, h and N, which respectively determine the resolution and upper limit of the integration grid. These can be modified to accurately transform any function that theoretically converges. How to choose these values, and the estimated error of the transform under a given choice, are discussed in the hankel’s extensive online documentation (and the reader is referred to Ogata (2005) for more details). Based on the arguments in the documentation, hankel provides an automatic, guided-adaptive algorithm for determination of h and N.

A particularly important limitation of hankel, as currently implemented, is that it does not implement the discrete Hankel transform. That is, it provides no direct means of transforming an array of regular-spaced function values into “radial Fourier space” at regular-spaced k values. It is focused solely on transforming callable functions, so that it can evaluate that function at the non-regular locations required by the double-exponential transform of Ogata (2005). Extensions to discrete Hankel transforms (and even fast Hankel transforms) are envisioned for v2.0 of hankel.

Acknowledgements

The authors acknowledge helpful contributions from Sebastian Mueller during the construction of this code. Parts of this research were supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. FJP would like to thank NSERC for research funding during the time of this research.