blitz
blitz copied to clipboard
patch: algorithm for computation of Bessel function
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 :
- the computer code for computation of Bessel functions of complex arguments and integer order
- the test suite for the algorithm
- 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.