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.

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

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?

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).

No one that knows what they are doing will compute Hankel functions with strictly forward or backwards recursion. Fortunately, the Hankel function of first kind is h_n = cmplx(j_n,y_n), and so is trivial to compute.

Also note the claim of unconditionally unstable is a bit strong. For example, one can use forward recursion for j_n provided n < x. Why would one want to use forward recursion if possible, well speed of course.

program van
implicit none
integer, parameter :: n = 20, nmax = 10
integer i
real jfwd(0:n), jbwd(0:n), x, ix, jmx1, jmx, j0, j1
x = 10
ix = 1 / x
j0 = sin(x) * ix
j1 = (j0 - cos(x)) * ix
!
! Forward recursion.
!
jfwd(0) = j0
jfwd(1) = j1
do i = 1, n - 1
jfwd(i+1) = (2 * i + 1) * ix * jfwd(i) - jfwd(i-1)
end do
!
! Backward recursion.
!
jmx1 = 0
jmx = 1
do i = n + nmax, n + 2, -1
jbwd(0) = (2 * i + 1) * (ix * jmx - jmx1 / (2 * i + 1))
jmx1 = jmx
jmx = jbwd(0)
end do
do i = n + 1, 1, -1
jbwd(i - 1) = (2 * i + 1) * ix * jmx - jmx1
jmx1 = jmx
jmx = jbwd(i - 1)
end do
if (abs(jmx) >= abs(jmx1)) then
j0 = j0 / jbwd(0)
else
j0 = j1 / jbwd(1)
end if
jbwd = j0 * jbwd
!
! https://keisan.casio.com/exec/system/1180573464
!
! j(0) = -5.44021110889369813405E-2
! j(1) = 7.84669417987515470918E-2
! ...
! j(9) = 1.00096409548490576567E-2
! j(10) = 6.46051544925642642714E-2
! j(11) = 3.55744148858943784028E-2
! ...
! j(20) = 2.30837196131946871671E-6
!
do i = 0, n
write(*,'(I2, 2ES15.6)') i, jfwd(i), jbwd(i)
end do
end program van

I leave it as an exercise to others to copy-n-paste to test.