@cindex discrete Hankel transforms
@cindex Hankel transforms, discrete
@cindex transforms, Hankel
This chapter describes functions for performing Discrete Hankel
Transforms (DHTs). The functions are declared in the header file
@file{gsl_dht.h}.

@menu
* Discrete Hankel Transform Definition::  
* Discrete Hankel Transform Functions::  
* Discrete Hankel Transform References::  
@end menu

@node Discrete Hankel Transform Definition
@section Definitions

The discrete Hankel transform acts on a vector of sampled data, where
the samples are assumed to have been taken at points related to the
zeros of a Bessel function of fixed order; compare this to the case of
the discrete Fourier transform, where samples are taken at points
related to the zeroes of the sine or cosine function.

Starting with its definition, the Hankel transform (or Bessel transform) of
order @math{\nu} of a function @math{f} with @math{\nu > -1/2} is defined as
(see Johnson, 1987 and Lemoine, 1994)
@tex
\beforedisplay
$$
F_\nu(u) = \int_0^\infty f(t) J_\nu(u t) t dt
$$
\afterdisplay
@end tex
@ifinfo
@example
F_\nu(u) = \int_0^\infty f(t) J_\nu(u t) t dt
@end example
@end ifinfo

@noindent
If the integral exists, @math{F_\nu} is called the Hankel transformation
of @math{f}. The reverse transform is given by
@tex
\beforedisplay
$$
f(t) = \int_0^\infty F_\nu(u) J_\nu(u t) u du ,
$$
\afterdisplay
@end tex
@ifinfo
@example
f(t) = \int_0^\infty F_\nu(u) J_\nu(u t) u du ,
@end example
@end ifinfo

@noindent
where @math{\int_0^\infty f(t) t^{1/2} dt} must exist and be
absolutely convergent, and where @math{f(t)} satisfies Dirichlet's
conditions (of limited total fluctuations) in the interval
@math{[0,\infty]}.

Now the discrete Hankel transform works on a discrete function
@math{f}, which is sampled on points @math{n=1...M} located at
positions @math{t_n=(j_{\nu,n}/j_{\nu,M}) X} in real space and
at @math{u_n=j_{\nu,n}/X} in reciprocal space. Here,
@math{j_{\nu,m}} are the m-th zeros of the Bessel function
@math{J_\nu(x)} arranged in ascending order. Moreover, the
discrete functions are assumed to be band limited, so
@math{f(t_n)=0} and @math{F(u_n)=0} for @math{n>M}. Accordingly,
the function @math{f} is defined on the interval @math{[0,X]}.

Following the work of Johnson, 1987 and
Lemoine, 1994, the discrete Hankel transform is given by
@tex
\beforedisplay
$$
F_\nu(u_m) = {{2 X^2}\over{j_{\nu,M}^2}}
      \sum_{k=1}^{M-1} f\left({{j_{\nu,k} X}\over{j_{\nu,M}}}\right)
          {{J_\nu(j_{\nu,m} j_{\nu,k} / j_{\nu,M})}\over{J_{\nu+1}(j_{\nu,k})^2}}.
$$
\afterdisplay
@end tex
@ifinfo
@example
F_\nu(u_m) = (2 X^2 / j_(\nu,M)^2)
      \sum_@{k=1@}^@{M-1@} f(j_(\nu,k) X/j_(\nu,M))
          (J_\nu(j_(\nu,m) j_(\nu,k) / j_(\nu,M)) / J_(\nu+1)(j_(\nu,k))^2).
@end example

@end ifinfo
@noindent
It is this discrete expression which defines the discrete Hankel
transform calculated by GSL. In GSL, forward and backward transforms
are defined equally and calculate @math{F_\nu(u_m)}.
Following Johnson, the backward transform reads
@tex
\beforedisplay
$$
f(t_k) = {{2}\over{X^2}}
      \sum_{m=1}^{M-1} F\left({{j_{\nu,m}}\over{X}}\right)
          {{J_\nu(j_{\nu,m} j_{\nu,k} / j_{\nu,M})}\over{J_{\nu+1}(j_{\nu,m})^2}}.
$$
\afterdisplay
@end tex
@ifinfo
@example
f(t_k) = (2 / X^2)
      \sum_@{m=1@}^@{M-1@} F(j_(\nu,m)/X)
          (J_\nu(j_(\nu,m) j_(\nu,k) / j_(\nu,M)) / J_(\nu+1)(j_(\nu,m))^2).
@end example

@end ifinfo
@noindent
Obviously, using the forward transform instead of the backward transform gives an
additional factor @math{X^4/j_{\nu,M}^2=t_m^2/u_m^2}.

The kernel in the summation above defines the matrix of the
@math{\nu}-Hankel transform of size @math{M-1}. The coefficients of
this matrix, being dependent on @math{\nu} and @math{M}, must be
precomputed and stored; the @code{gsl_dht} object encapsulates this
data. The allocation function @code{gsl_dht_alloc} returns a
@code{gsl_dht} object which must be properly initialized with
@code{gsl_dht_init} before it can be used to perform transforms on data
sample vectors, for fixed @math{\nu} and @math{M}, using the
@code{gsl_dht_apply} function. The implementation allows to define the
length @math{X} of the fundamental interval, for convenience, while
discrete Hankel transforms are often defined on the unit interval
instead of @math{[0,X]}.

Notice that by assumption @math{f(t)} vanishes at the endpoints
of the interval, consistent with the inversion formula
and the sampling formula given above. Therefore, this transform
corresponds to an orthogonal expansion in eigenfunctions
of the Dirichlet problem for the Bessel differential equation.

@node Discrete Hankel Transform Functions
@section Functions

@deftypefun {gsl_dht *} gsl_dht_alloc (size_t @var{size})
@tindex gsl_dht
This function allocates a Discrete Hankel transform object of size
@var{size}.
@end deftypefun

@deftypefun int gsl_dht_init (gsl_dht * @var{t}, double @var{nu}, double @var{xmax})
This function initializes the transform @var{t} for the given values of
@var{nu} and @var{xmax}.
@end deftypefun

@deftypefun {gsl_dht *} gsl_dht_new (size_t @var{size}, double @var{nu}, double @var{xmax})
This function allocates a Discrete Hankel transform object of size
@var{size} and initializes it for the given values of @var{nu} and
@var{xmax}.
@end deftypefun

@deftypefun void gsl_dht_free (gsl_dht * @var{t})
This function frees the transform @var{t}.
@end deftypefun

@deftypefun int gsl_dht_apply (const gsl_dht * @var{t}, double * @var{f_in}, double * @var{f_out})
This function applies the transform @var{t} to the array @var{f_in}
whose size is equal to the size of the transform.  The result is stored
in the array @var{f_out} which must be of the same length.   

Applying this function to its output gives the original data
multiplied by @c{$(X^2/j_{\nu,M})^2$}
@math{(1/j_(\nu,M))^2}, up to numerical errors.
@end deftypefun

@deftypefun double gsl_dht_x_sample (const gsl_dht * @var{t}, int @var{n})
This function returns the value of the @var{n}-th sample point in the unit interval,
@c{${({j_{\nu,n+1}} / {j_{\nu,M}}}) X$}
@math{(j_@{\nu,n+1@}/j_@{\nu,M@}) X}. These are the
points where the function @math{f(t)} is assumed to be sampled.
@end deftypefun

@deftypefun double gsl_dht_k_sample (const gsl_dht * @var{t}, int @var{n})
This function returns the value of the @var{n}-th sample point in ``k-space'',
@c{${{j_{\nu,n+1}} / X}$}
@math{j_@{\nu,n+1@}/X}.
@end deftypefun

@node Discrete Hankel Transform References
@section References and Further Reading

The algorithms used by these functions are described in the following papers,

@itemize @w{}
@item
H. Fisk Johnson, Comp.@: Phys.@: Comm.@: 43, 181 (1987).
@end itemize

@itemize @w{}
@item
D. Lemoine, J. Chem.@: Phys.@: 101, 3936 (1994).
@end itemize