gfortran和ivf有贝塞尔函数，球贝塞尔函数的话，梁昆淼的数学物理方法里面就有介绍，非常经典的教材
8.26 BESSEL_J0 — Bessel function of the first kind of order 0


Description:
BESSEL_J0(X) computes the Bessel function of the first kind of order 0 of X. This function is available under the name BESJ0 as a GNU extension.
Standard:
Fortran 2008 and later
Class:
Elemental function
Syntax:
RESULT = BESSEL_J0(X)
Arguments:
X The type shall be REAL, and it shall be scalar.


Return value:
The return value is of type REAL and lies in the range - 0.4027... \leq Bessel (0,x) \leq 1. It has the same kind as X.
Example:
program test_besj0
real(8) :: x = 0.0_8
x = bessel_j0(x)
end program test_besj0


Specific names:
Name Argument Return type Standard
DBESJ0(X) REAL(8) X REAL(8) GNU extension