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