Dear Aparna,

thank you for the answer.

Unfortunately, it is not helpful. I am not really looking to design an optical system; I am looking to gain insight about (strongly) nonlinear ultrashort pulsed beam propagation. I do not think Zemax could help me with this at the level of detail I am interested in.

What I need is information about how to efficiently compute a Hankel transform for implementation in my own modeling script.

Rimantas,

When calculated in a straightforward way, a 1-Dimensional discrete Fourier transform or Hankel transform requires of order O(N2) operations for a vector of length N.  For an N x N array, the operation count is O(N3)

However, Fast Fourier Transform algorithms exploit symmetries in the sine and cosine functions, re-ordering and combining operations in a hierarchical fashion for cases where N is a multiple of small prime powers.  In particular, the Cooley-Tukey algorithm is an O(N log2N) process when N is a power of 2.  Highly optimised implementations are available.  Operations count for an N x N array should be O(N2 log N).

For moderate size N, 2-dimensional FFT can be comparable with or faster than a 1-dimensional O(N2) Hankel transform, especially if the Bessel functions are re-calculated each time the transforms are called.

In your application you require forward and inverse transforms at each propagation step, so there are multiple calls to transforms of identical size.  Forward and reverse Hankel transforms can be reduced to multiplication of a length N field vector by an (N x N) array, whose values depend only on the chosen dimension and step size.  These can be calculated once at the start of the calculation, avoiding the need for expensive Bessel function calls at each propagation step.

For a moderate or large number of propagation steps, I would expect  1-dimensional Hankel transforms, using pre-calculated coefficient matrices, to out-perform 2-dimensional FFT.

With very high non-linearity there is a risk (depending on the beam quality and dispersion regime) that noise-seeded modulation instability will break the cylindrical symmetry.  If this is a concern, I would be inclined to revert to a 2-dimensional FFT approach.

Hope this helps.

https://en.wikipedia.org/wiki/Fast_Fourier_transform

http://www.fftw.org/

https://en.wikipedia.org/wiki/Hankel_transform