patch: algorithm for computation of Bessel function

[slayoo]

Migrated from SF: https://sourceforge.net/p/blitz/patches/7/

Patch proposed by "van den Bosch" (SF user: sedibald) back in 2003 with later modifications to the code by Julian Cummings (final version attached here: besselj_complex.tar.gz).

Extracts from the comments:

van den Bosch:

This tarball gives : 1) the computer code for computation of Bessel functions of complex arguments and integer order 2) the test suite for the algorithm 3) a MATLAB code for comparing bessel functions generated by MATLAB and the algorithm

The Bessel functions algorithm is written with the help of Blitz++, although this is not a requirement.

Julian:

I have converted the original Bessel function that was submitted into a functor, where operator() applies the function. This allows one to use simple syntax for applying the function to a blitz Array of complex numbers. In addition, I have templated the original code on the floating-point number type, so that both single- and double-precision complex numbers will be accepted. At the suggestion of the code author, I modified the original algorithm to ensure convergence when real(z) is negative. I have checked that the results and the speed of the algorithm did not change with these modifications. The besselj_complex function spends approximately 2.95 seconds computing in the test program test_besselj_complex, whereas MATLAB computes the function in about 1.62 seconds on my desktop machine.

van den Bosch:

I have changed the algorithm of the functor in order to compute besselj_complex for the whole complex z plane. Next changes will concern the stability of the algorithm for very high values of n (the order), as well as its profiling for the sake of performance (MATLAB is still 2-3 times faster). The programs needed for the tests are bundled in the tarball. Of course, one can also write its own test routine or simply edit the z_points.ini file.