hankel_matrix
``` hankel_matrix: Generates data to use for Hankel Transforms

s_HT = hankel_matrix(order, rmax, Nr, eps_roots)

s_HT	=	Structure containing data to use for the pQDHT
order	=	Transform order
rmax	=	Radial extent of transform
Nr		=	Number of sample points
eps_roots	=	Error in estimation of roots of Bessel function (optional)

s_HT:
order, rmax, Nr =	As above
J_roots		=	Roots of the pth order Bessel fn.
J_roots_N1	=	(N+1)th root
r			=	Radial co-ordinate vector
v			=	frequency co-ordinate vector
kr			=	Radial wave number co-ordinate vector
vmax		=	Limiting frequency
=	roots_N1 / (2*pi*rmax)
S			=	rmax * 2*pi*vmax (S product)
T			=	Transform matrix
J			=	Scaling vector
=	J_(order+1){roots}

The algorithm used is that from:
"Computation of quasi-discrete Hankel transforms of the integer
order for propagating optical wave fields"
Manuel Guizar-Sicairos and Julio C. Guitierrez-Vega
J. Opt. Soc. Am. A 21(1) 53-58 (2004)

The algorithm also calls the function:
zn = bessel_zeros(1, p, Nr+1, 1e-6),
where p and N are defined above, to calculate the roots of the bessel
function. This algorithm is taken from:
"An Algorithm with ALGOL 60 Program for the Computation of the
zeros of the Ordinary Bessel Functions and those of their
Derivatives".
N. M. Temme
Journal of Computational Physics, 32, 270-279 (1979)

Example: Propagation of radial field

% Note the use of matrix and element products / divisions
H = hankel_matrix(0, 1e-3, 512);
DR0 = 50e-6;
Ur0 = exp(-(H.r/DR0).^2);
Ukr0 = H.T * (Ur0./H.J);
k0 = 2*pi/800e-9;
kz = realsqrt((k0^2 - H.kr.^2).*(k0>H.kr));
z = (-5e-3:1e-5:5e-3);
Ukrz = (Ukr0*ones(1,length(z))).*exp(i*kz*z);
Urz = (H.T * Ukrz) .* (H.J * ones(1,length(z)));```