This file contains 2 functions:
bessel_zeros  finds zeros of a bessel function
hankel_matrix  generates a structure of data to use for Hankel transforms
The algorithm is taken from:
M. GuizarSicairos and J. C. GutierrezVega, Computation of quasidiscrete Hankel transforms of integer order for propagating optical wave fields,
J. Opt. Soc. Am. A 21, 5358 (2004)
An example of beam propagation with circular symmetry has been written in the hankel_matrix help:
Type "help hankel_matrix" or "doc hankel_matrix" in Matlab.
1.1.0.0  Slight changes to the hankel structure to simplify usage. 

1.0.0.0  Minor modifications. 
Inspired by: Integer order Hankel transform
Create scripts with code, output, and formatted text in a single executable document.
Garima Nagar (view profile)
I want to do discrete Hankel transformation of the 2D image/ 2D matrix. How can I use this function for the 2D matrix?
Thanks!
anping Yu (view profile)
RAVI KUMAR (view profile)
how can we use it for the Hankel transformation of 2D image?
Ahmed Fasih (view profile)
Answered my own question: https://gist.github.com/fasiha/16f796be6a55cc7dae98 but essentially, to produce Fig 1(c):
h = Hankel.hankel_matrix(4, 3, 256, 1e13);
gamma = 5;
f1fun = @(r) sinc(2*gamma*r); % our sinc includes pi
r = h.J_roots / (2*pi*h.vmax);
F1 = f1fun(r) ./ h.J * h.vmax; F2 = h.T * F1;
f2 = F2 .* h.J / h.vmax;
v = h.J_roots / (2*pi*h.rmax);
figure; plot(v, f2)
(See the link for properlycommented version! Don't write code like this!)
Ahmed Fasih (view profile)
Can you show how to use this function to reproduce the sinc plot, Figure 1 in GuizarSicairos & GutierrezVega?
Ismail (view profile)
Thanks for this submission.
I have a question: can one uses this code for non normalized parameters?
Adam Wyatt (view profile)
The example given in the code does not scale the amplitudes of the angular distribution because it is not necessary (since it is only necessary to multiply by a phase factor and so the scaling required will be undone again before transforming back). However, strictly speaking this should be done as follows:
H = hankel_matrix(0, rmax, samples);
HT = @(f) (H.T*(f./H.JR)).*H.JV;
IHT = @(F) (H.T*(F./H.JV)).*H.JR;
Note that two scaling vectors are supplied:
H.JR
H.JV
Therefore I agree that the example is misguiding from a physical point of view, but not computationally, but it is not incomplete.
Note that there is a difference between the Hankel transform (HT) and the quasi discrete Hankel transform (qDHT). This code provides a way of performing the qDHT and the scaling vectors to turn it into an accurate estimation of the HT.
Indeed there is a factor of 2pi, but as stated this does depend on the definition of the HT used  so long as f = IHT(HT(f)), its all good. So using the definition consistent with the code:
HT{exp(br^2)} = pi/b * exp{k^2/(2*b)];
Further details can be found in the references given in the code.
Roy (view profile)
The example given inside hankel_matrix.m is incomplete and misguiding. According to it, one can define the Hankel and inverseHankel transform as follows:
mat_H = hankel_matrix(0, rmax, samples);
ht = @(f) mat_H.T * (f./mat_H.J);
iht = @(F) (mat_H.T * F) .* mat_H.J;
Recall that the Hankel transform of a Gaussian is a Gaussian:
H{e(br^2)} = 1/2b * e(k^2/4b)
However, trying to match this analytic result using the above transform fails. This is because the scaling vector J is defined as Jp1; a possible solution would be to define two vectors, Jp1/rmax and Jp1/vmax, and multiply f before and after the transform, and conversely in the inverse transform.
There is also a factor of 2*pi, but this depends on your formalism of choice.
Kevin (view profile)
This is a nice submission.