Fortran 2008 provides regular Bessel functions.
I’m looking for an implementation of spherical Bessel and Hankel functions. I need them for both first and second kind, and of arbitrary order (accepts order as input).
I have what seems to be a F66 implementation in a scan of an NRL report by Blake (1972). However, I can’t find the code elsewhere and transcribing it would be error prone.
Do you know of a library that implements these?
I don’t know of any libraries. I can see that on p438 (10.1.10) of Abrimowitz and Stegun there is a recursion formula for Spherical Bessel Functions in using trigonometric functions which seems to be very simple.
My thought is not so much that you use that to evaluate the function directly since there could be a serious performance hit doing that, rather that you use might be able to use it in a script to generate the Fortran code for each order n. Might be an interesting case for using parameter packs in C++
Chapter 8 of the book Computation of Special Functions by Shanjie Zhang and Jian-Ming Jin is on Spherical Bessel Functions. John Burkardt has the Fortran 77 code, and his free source form version is here with driver here, giving results
MSPHJ
SPHJ evaluates spherical Bessel J functions
n x jn(x) jn'(x)
0 0.9050000000D+00 0.8689780718D+00 -0.2776712617D+00
1 0.9050000000D+00 0.2776712617D+00 0.2553399244D+00
2 0.9050000000D+00 0.5147914933D-01 0.1070221479D+00
3 0.9050000000D+00 0.6743927972D-02 0.2167173288D-01
4 0.9050000000D+00 0.6838294584D-03 0.2965864666D-02
5 0.9050000000D+00 0.5658597917D-04 0.3086737954D-03
n x jn(x) jn'(x)
0 0.1000000000D+02 -0.5440211109D-01 -0.7846694180D-01
1 0.1000000000D+02 0.7846694180D-01 -0.7009549945D-01
2 0.1000000000D+02 0.7794219363D-01 0.5508428371D-01
3 0.1000000000D+02 -0.3949584498D-01 0.9374053162D-01
4 0.1000000000D+02 -0.1055892851D+00 0.1329879757D-01
5 0.1000000000D+02 -0.5553451162D-01 -0.7226857814D-01
I mean arbitrary values of integer order (i.e. n = 0, 1, 2, etc.)
I need real only.
Fortran needs a maintained version of all commonly used special functions. Either as a standalone library or part of stdlib.
I found that scipy.special provides spherical Bessel functions based on Zhang and Jin. The fixed-form Fortran implementations are here (j_n, first kind) and here (y_n, second kind).
The range is [10^{-2}, 10]. Black (1972), page 13 says that for argument x = 10.88, order n = 83 had the last term in the summation < 10^{-6}, which is more than sufficient for my application.
@kargl thanks a lot, that would be fantastic. It’s not urgent. Do you see a good reason for me to not use the subroutines from specfun in scipy.special?
Last time I used them, they sometimes give incorrect results for certain parameter values. But I can’t find my issue about it anymore.
I have a code that computes j_n ( m chi ) and y_n ( chi ) for double precision chi and complex m. It’s based upon a method of Infeld and a routine DBMIE by J. V. Dave (1968) – used for Mie scattering – and adapted with extensions by C. D. Cantrell (1988). I upgraded it to Fortran 2000 in 2007. I made a simplistic (not comprehensive) comparison to reference codes for arbitrary order Bessel and Neumann functions. The maximum absolute error for j_n (x) with log10(x) = (-1(0.2)2), n = 1…100 and chi=1.0 (i.e., real j_n only) was 1.65e-15. The maximum relative error for y_n was 1.33e-12.
Wrt the library by Zhang and JIn: I had some problems in the past with evaluating the modified Bessel function K (never sure of the exact denomination, I use them too seldom). The routine did not converge properly for some range of the argument and came up with Nan, if I remember correctly. That does not mean that the library is no use, of course, far from that.
FWIW, I found code for the complex gamma function that seems fairly accurate. I checked a few data points using the relation gamma(x) = (x-1)*gamma(x-1). The relative errors were in the order of 1.0e-15.
(Not a function I really need, but I was curious as the Fortran standard does not define the gamma function for complex arguments)
The SLATEC library (SLATEC - A Mathematical Library) contains several routines for the Bessel functions.
Alan Miller’s Fortran Software (Alan Miller's Fortran Software) also contains links for several implementations of the Bessel functions suite.
j_n is unconditionally unstable for forward recursion. y_n is unconditionally unstable for backward recursion. So if you compute Hankel functions by forward (or reverse) recursion, you get (y_n, y_n) or (j_n,j_n).