Adam Wyatt

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.

I wrote my own program, starting with this one. (You can find it posted under the same title as this one). I found this program to be a nice starting point to use. I found plenty of room for improvements and additions, but I don't remember finding any bugs that would noticeably impact the results.

Hey Andreas, this is what I am searching for.
I runned your package with 'bsxfun'. When the vector is small (or the number of n in your example) it is OK. However when the vector is huge(let's say n>>100) the program ran out of memory. Is there a way to fix it?

This one worked like a charm for a couple of integrals where quadv failed (and iterating over quadgk was too slow).

The file doesn't work if you follow the example though. It should be something like:

Example:
 Y = @(x,n) 1./(bsxfun(@plus,(1:n)',x));
 Qv = quadvgk(@(x) Y(x,10), [0;1], 10)