...one of the most highly
regarded and expertly designed C++ library projects in the
world.

— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards

`#include <boost/math/special_functions/bessel.hpp>`

template <class T1, class T2>calculated-result-typesph_bessel(unsigned v, T2 x); template <class T1, class T2, class Policy>calculated-result-typesph_bessel(unsigned v, T2 x, const Policy&); template <class T1, class T2>calculated-result-typesph_neumann(unsigned v, T2 x); template <class T1, class T2, class Policy>calculated-result-typesph_neumann(unsigned v, T2 x, const Policy&);

The functions sph_bessel and sph_neumann return the result of the Spherical Bessel functions of the first and second kinds respectively:

sph_bessel(v, x) = j_{v}(x)

sph_neumann(v, x) = y_{v}(x) = n_{v}(x)

where:

The return type of these functions is computed using the *result
type calculation rules* for the single argument type T.

The final Policy argument is optional and can be used to control the behaviour of the function: how it handles errors, what level of precision to use etc. Refer to the policy documentation for more details.

The functions return the result of domain_error
whenever the result is undefined or complex: this occurs when `x < 0`

.

The j_{v} function is cyclic like J_{v} but differs in its behaviour at the origin:

Likewise y_{v} is also cyclic for large x, but tends to -∞
for small *x*:

There are two sets of test values: spot values calculated using functions.wolfram.com, and a much larger set of tests computed using a simplified version of this implementation (with all the special case handling removed).

Other than for some special cases, these functions are computed in terms of cyl_bessel_j and cyl_neumann: refer to these functions for accuracy data.

Other than error handling and a couple of special cases these functions are implemented directly in terms of their definitions:

The special cases occur for:

j_{0} = sinc_pi(x) = sin(x)
/ x

and for small *x < 1*, we can use the series:

which neatly avoids the problem of calculating 0/0 that can occur with the main definition as x → 0.