Algorithm 794: Numerical Hankel Transform by the Fortran Program HANKEL.(Statistical Data Included)

1. PROBLEMS WITH THE HANKEL TRANSFORM

The real Hankel transform H[Nu]([Xi], f(x)) of a real function y = f(x) requires the evaluation of the infinite integral

(1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [Nu] [element of] R is the order of the transform; [Xi] [element of] R is the transformation parameter with 0 [is less than] [Xi]; and [J.sub.[Nu]](x) is the Bessel function of the first kind. Precise definitions of [H.sub.[Nu]]([Xi], f(x)) and conditions on f(x) can be found in Erdelyi et al. [1954], Korn and Korn [1968], and Sneddon [1955]. For -1/2 [is less than] v the Hankel transform is self-reciprocal, i.e., with g([Xi]) = [H.sub.[Nu]]([Epsilon], f(x)) given one finds f(x) from [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. The integral (1) will exist if f(x) [element of] [L.sub.1], i.e., if the integral [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] exists. Analytical solutions of (1) are known for a variety of functions f(x) [Erdelyi et al. 1954]. By use of a symbolic computation program further analytical solutions may be found.(1)

Numerical solutions also are of interest and unavoidable, if f(x) is not known analytically. For example, measured data may be given (sampled) as a table of N data pairs ([x.sub.1], f([x.sub.1])), ..., ([x.sub.N], f(x.sub.N])). Several computer programs are available to solve (1) numerically. Siegman used a nonlinear change of variables to convert the one-sided Hankel transform into a two-sided cross-correlation integral [Siegman 1980]. The algorithm is particularly fast and can be applied to sampled data, but requires a sampling at exponentially increasing x-values. Piessens presented a program which solves (1) making use of numerical integration [Piessens 1982]. The approximation for the Bessel function Piessens used allowed him to write the infinite integral as the sum of an integral over a finite interval and of a Fourier-sine and Fourier-cosine transform. The restriction to integer v-values (with 0 [is less than or equal to] [Nu] [is less than or equal to] 10) prevents the application of that program to transforms with noninteger v which in particular arises if spherical Bessel functions are involved. In the same year, Anderson [1982] introduced his algorithm for the Hankel transform by using related and lagged convolutions. Again, the algorithm is restricted to integer [Nu]. The noninteger Hankel transform with [Nu] = j + 1/2 (where j is an integer) is known as spherical Bessel transform. Talman [1983] provided a corresponding program for sampled data, but the x-values have to be distributed uniformly in 1n(x).

In this work, a direct numerical approach to solve Eq. (1) is presented which consists of two parts: (1) the calculation of the infinite integral and (2) the calculation of the Bessel function [J.sub.[Nu]](x).

(1) Several strategies can be found in literature for the numerical evaluation of infinite integrals. None of them is strictly valid, since any numerical procedure has to map the infinite integral onto a finite one. A frequent strategy is …