Hankel Transform Package
The Hankel Transform computes the Hankel (Fourier-Bessel) transform of order nu,
F(k) = Integral from 0 to infinity of f(r) J_nu(k r) r dr
where J_nu is the Bessel function of the first kind. This is the radial part of the
two-dimensional Fourier transform of a circularly symmetric function, so it shows up
wherever you have data on a radius: diffraction and beam propagation, small-angle
scattering, tomographic and Abel-type inversions, and radially symmetric PDE work.
The transform is its own inverse: feeding F(k) back through the same routine with the
same order recovers f(r).
The package consists of plain Igor procedure code. There is no XOP to
install. The procedure contains several implementations of the transform using different
approaches. The main transform function:
Function/WAVE HankelTransform(WAVE yw, Variable nu, [WAVE xw, WAVE kw,
Variable numK, Variable kMax, String dest])
- `yw` is the 1D real input f(r). Radii come from the wave's X scaling unless `xw` is
supplied.
- `nu` is the order. It may be non-integer or negative, since the underlying `bessJ`
accepts a real order.
- `xw` supplies the radii explicitly, so non-uniformly sampled data is handled directly.
- `kw` supplies arbitrary output wavenumbers.
- `numK` and `kMax` set the output grid when `kw` is omitted. The default `kMax` is
pi/dr, the Nyquist limit of the input sampling.
- `dest` names the output wave. Without it a free wave is returned, so assign the result
to a wave reference.
When the output grid is uniform the result carries X scaling in k, so it plots and
evaluates directly against k.
The integral is evaluated by direct quadrature. Trapezoidal weights are built from the
actual radial coordinates, and the kernel J_nu(k_i r_j) is formed as an m x n matrix and
applied with `MatrixOP`, so time and memory both scale as m*n. 512 x 512 is essentially
instantaneous; 4096 x 4096 needs about 134 MB for the kernel.
Because the integrand carries a factor of r it vanishes at r = 0, so if the input range
is chosen such that f has also decayed at the far end, trapezoidal quadrature converges
very rapidly. For smooth, well-decayed input the result is accurate to near machine
precision. Against the analytic order-0 Gaussian pair
exp(-a r^2) <-> exp(-k^2/(4a))/(2a)
the forward transform and the round trip both agree to about 1e-9 at n = 512.
Example
Variable a = 2, n = 512, rMax = 6
Make/O/D/N=(n) gaussIn
SetScale/P x, 0, rMax/(n-1), "", gaussIn
gaussIn = exp(-a*x^2)
WAVE ht = HankelTransform(gaussIn, 0, numK=n, kMax=20, dest="gaussOut")
Duplicate/O ht, gaussExact
gaussExact = exp(-x^2/(4*a))/(2*a)
`DemoHankelTransform()` runs exactly this check and then round-trips the result back to
r space, printing both error figures to the history.
Accuracy notes
Truncation (not quadrature) is normally the dominant error. Choose the input range so
that f has genuinely decayed before the last point, because the discarded tail
contributes directly to the answer.
Asking for k beyond the default `kMax` asks for detail the input sampling does not
contain.
Discontinuous inputs converge only as O(dr) and ring. A circular top hat, whose
transform is R J_1(kR)/k, is the standard example. For such functions either oversample
heavily or use a method that samples at the zeros of J_nu.
Also included
`HankelTransformFFT()` implements Siegman's quasi-fast Hankel transform (Opt. Lett. 1,
13 (1977)), an O(N log N) method on logarithmic r and k grids. It is useful when the
O(mn) direct routine is too slow, with the usual caveats of a log grid: r = 0 is
unreachable, kMax/kMin is forced to equal rMax/rMin, and the oscillatory kernel aliases
once its phase step exceeds about 1 radian per sample. The routine warns when that
happens.
`HankelTransformFFTLog()` is an in-progress implementation of FFTLog (Talman 1978;
Hamilton, MNRAS 312, 257 (2000)). It is included for reference along with its diagnostic
routines, but it is not yet verified and should not be relied on. Use
`HankelTransform()` for production work.
References
Bracewell, R.N., *The Fourier Transform and Its Applications*, 3rd ed., McGraw-Hill,
2000. Chapter 13 covers the Hankel transform and its relationship to the 2D Fourier
transform of a circularly symmetric function.
Gradshteyn, I.S., and I.M. Ryzhik, *Table of Integrals, Series, and Products*, 7th ed.,
Academic Press, 2007. Section 6.5 tabulates Hankel transform pairs useful for
validation.
Siegman, A.E., "Quasi fast Hankel transform," *Opt. Lett. 1, 13 (1977).
Hamilton, A.J.S., "Uncorrelated modes of the non-linear power spectrum," MNRAS
312, 257 (2000).