HANKEL - Functions to compute discrete Hankel transforms
hankelalloc allocate and return a pointer to a Hankel transformer
hankelfree free a Hankel transformer
hankel0 compute the zeroth-order Hankel transform
hankel1 compute the first-order Hankel transform
Function Prototypes:
void *hankelalloc (int nfft);
void hankelfree (void *ht);
void hankel0 (void *ht, float f[], float h[]);
void hankel1 (void *ht, float f[], float h[]);
hankelalloc:
Input:
nfft valid length for real to complex fft (see notes below)
Returned:
pointer to Hankel transformer
hankelfree:
Input:
ht pointer to Hankel transformer (as returned by hankelalloc)
hankel0:
Input:
ht pointer to Hankel transformer (as returned by hankelalloc)
f array[nfft/2+1] to be transformed
Output:
h array[nfft/2+1] transformed
hankel1:
Input:
ht pointer to Hankel transformer (as returned by hankelalloc)
f array[nfft/2+1] to be transformed
Output:
h array[nfft/2+1] transformed
Notes:
The zeroth-order Hankel transform is defined by:
Infinity
h0(k) = Integral dr r j0(k*r) f(r)
0
where j0 denotes the zeroth-order Bessel function.
The first-order Hankel transform is defined by:
Infinity
h1(k) = Integral dr r j1(k*r) f(r)
0
where j1 denotes the first-order Bessel function.
The Hankel transform is its own inverse.
The Hankel transform is computed by an Abel transform followed by
a Fourier transform.
References:
Hansen, E. W., 1985, Fast Hankel transform algorithm: IEEE Trans. on
Acoustics, Speech and Signal Processing, v. ASSP-33, n. 3, p. 666-671.
(Beware of several errors in the equations in this paper!)
Authors: Dave Hale, Colorado School of Mines, 06/04/90