
[section:bessel Bessel Functions of the First and Second Kinds]

[h4 Synopsis]

   template <class T1, class T2>
   ``__sf_result`` cyl_bessel_j(T1 v, T2 x);

   template <class T1, class T2, class ``__Policy``>
   ``__sf_result`` cyl_bessel_j(T1 v, T2 x, const ``__Policy``&);

   template <class T1, class T2>
   ``__sf_result`` cyl_neumann(T1 v, T2 x);
   
   template <class T1, class T2, class ``__Policy``>
   ``__sf_result`` cyl_neumann(T1 v, T2 x, const ``__Policy``&);
   
   
[h4 Description]

The functions __cyl_bessel_j and __cyl_neumann return the result of the
Bessel functions of the first and second kinds respectively:

cyl_bessel_j(v, x) = J[sub v](x)

cyl_neumann(v, x) = Y[sub v](x) = N[sub v](x)

where:

[equation bessel2]

[equation bessel3]

The return type of these functions is computed using the __arg_pomotion_rules
when T1 and T2 are different types.  The functions are also optimised for the
relatively common case that T1 is an integer.

[optional_policy]

The functions return the result of __domain_error whenever the result is
undefined or complex.  For __cyl_bessel_j this occurs when `x < 0` and v is not
an integer, or when `x == 0` and `v != 0`.  For __cyl_neumann this occurs
when `x <= 0`.

The following graph illustrates the cyclic nature of J[sub v]:

[graph cyl_bessel_j]

The following graph shows the behaviour of Y[sub v]: this is also
cyclic for large /x/, but tends to -[infin][space] for small /x/:

[graph cyl_neumann]

[h4 Testing]

There are two sets of test values: spot values calculated using
[@http://functions.wolfram.com 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).

[h4 Accuracy]

The following tables show how the accuracy of these functions
varies on various platforms, along with comparisons to the __gsl and
__cephes libraries.  Note that the cyclic nature of these
functions means that they have an infinite number of irrational
roots: in general these functions have arbitrarily large /relative/
errors when the arguments are sufficiently close to a root.  Of
course the absolute error in such cases is always small.
Note that only results for the widest floating-point type on the 
system are given as narrower types have __zero_error.  All values
are relative errors in units of epsilon.

[table Errors Rates in cyl_bessel_j
[[Significand Size] [Platform and Compiler] [J[sub 0][space] and J[sub 1]] [J[sub v]] [J[sub v][space] (large values of x > 1000)] ]
[[53] [Win32 / Visual C++ 8.0] 
      [Peak=2.5 Mean=1.1 

GSL Peak=6.6 

__cephes Peak=2.5 Mean=1.1] 
      [Peak=11 Mean=2.2 

GSL Peak=11 

__cephes Peak=17 Mean=2.5] 
      [Peak=413 Mean=110 

GSL Peak=6x10[super 11] 

__cephes Peak=2x10[super 5] ] ]
[[64] [Red Hat Linux IA64 / G++ 3.4] [Peak=7 Mean=3] [Peak=117 Mean=10] [Peak=2x10[super 4][space] Mean=6x10[super 3]] ]
[[64] [SUSE Linux AMD64 / G++ 4.1] [Peak=7 Mean=3] [Peak=400 Mean=40] [Peak=2x10[super 4][space] Mean=1x10[super 4]] ]
[[113] [HP-UX / HP aCC 6] [Peak=14 Mean=6] [Peak=29 Mean=3] [Peak=2700 Mean=450] ]
]

[table Errors Rates in cyl_neumann
[[Significand Size] [Platform and Compiler] [J[sub 0][space] and J[sub 1]] [J[sub n] (integer orders)] [J[sub v] (fractional orders)] ]
[[53] [Win32 / Visual C++ 8.0] 
      [Peak=330 Mean=54 

GSL Peak=34 Mean=9 

__cephes Peak=330 Mean=54] 
      [Peak=923 Mean=83 

GSL Peak=500 Mean=54 

__cephes Peak=923 Mean=83] 
      [Peak=561 Mean=36 

GSL Peak=1.4x10[super 6][space] Mean\=7x10[super 4][space] 

__cephes Peak=+INF]]
[[64] [Red Hat Linux IA64 / G++ 3.4] [Peak=470 Mean=56] [Peak=843 Mean=51] [Peak=741 Mean=51] ]
[[64] [SUSE Linux AMD64 / G++ 4.1] [Peak=1300 Mean=424] [Peak=2x10[super 4][space] Mean=8x10[super 3]] [Peak=1x10[super 5][space] Mean=6x10[super 3]] ]
[[113] [HP-UX / HP aCC 6] [Peak=180 Mean=63] [Peak=340 Mean=150] [Peak=2x10[super 4][space] Mean=1200] ]
]

Note that for large /x/ these functions are largely dependent on
the accuracy of the `std::sin` and `std::cos` functions.

Comparison to GSL and __cephes is interesting: both __cephes and this library optimise
the integer order case - leading to identical results - simply using the general 
case is for the most part slightly more accurate though, as noted by the 
better accuracy of GSL in the integer argument cases.  This implementation tends to
perform much better when the arguments become large, __cephes in particular produces
some remarkably inaccurate results with some of the test data (no significant figures
correct), and even GSL performs badly with some inputs to J[sub v].  Note that
by way of double-checking these results, the worst performing __cephes and GSL cases
were recomputed using [@http://functions.wolfram.com functions.wolfram.com],
and the result checked against our test data: no errors in the test data were found.

[h4 Implementation]

The implementation is mostly about filtering off various special cases:

When /x/ is negative, then the order /v/ must be an integer or the
result is a domain error.  If the order is an integer then the function
is odd for odd orders and even for even orders, so we reflect to /x > 0/.

When the order /v/ is negative then the reflection formulae can be used to
move to /v > 0/:

[equation bessel9]

[equation bessel10]

Note that if the order is an integer, then these formulae reduce to:

J[sub -n] = (-1)[super n]J[sub n]

Y[sub -n] = (-1)[super n]Y[sub n]

However, in general, a negative order implies that we will need to compute
both J and Y.

When /x/ is large compared to the order /v/ then the asymptotic expansions
for large /x/ in M. Abramowitz and I.A. Stegun, 
['Handbook of Mathematical Functions] 9.2.19 are used 
(these were found to be more reliable
than those in A&S 9.2.5).

When the order /v/ is an integer the method first relates the result
to J[sub 0], J[sub 1], Y[sub 0][space] and Y[sub 1][space] using either
forwards or backwards recurrence (Miller's algorithm) depending upon which is stable.
The values for J[sub 0], J[sub 1], Y[sub 0][space] and Y[sub 1][space] are 
calculated using the rational minimax approximations on
root-bracketing intervals for small ['|x|] and Hankel asymptotic
expansion for large ['|x|]. The coefficients are from:

W.J. Cody, ['ALGORITHM 715: SPECFUN - A Portable FORTRAN Package of 
Special Function Routines and Test Drivers], ACM Transactions on Mathematical 
Software, vol 19, 22 (1993). 

and

J.F. Hart et al, ['Computer Approximations], John Wiley & Sons, New York, 1968. 

These approximations are accurate to around 19 decimal digits: therefore
these methods are not used when type T has more than 64 binary digits.

When /x/ is small, J[sub x][space] is best computed directly from the series:

[equation bessel2]

In the general case we compute J[sub v][space] and 
Y[sub v][space] simultaneously.

To get the initial values, let
[mu][space] = [nu] - floor([nu] + 1/2), then [mu][space] is the fractional part 
of [nu][space] such that
|[mu]| <= 1/2 (we need this for convergence later). The idea is to
calculate J[sub [mu]](x), J[sub [mu]+1](x), Y[sub [mu]](x), Y[sub [mu]+1](x)
and use them to obtain J[sub [nu]](x), Y[sub [nu]](x).

The algorithm is called Steed's method, which needs two
continued fractions as well as the Wronskian:

[equation bessel8]

[equation bessel11]

[equation bessel12]

See: F.S. Acton, ['Numerical Methods that Work],
    The Mathematical Association of America, Washington, 1997.

The continued fractions are computed using the modified Lentz's method
(W.J. Lentz, ['Generating Bessel functions in Mie scattering calculations
using continued fractions], Applied Optics, vol 15, 668 (1976)). 
Their convergence rates depend on ['x], therefore we need
different strategies for large ['x] and small ['x].

['x > v], CF1 needs O(['x]) iterations to converge, CF2 converges rapidly

['x <= v], CF1 converges rapidly, CF2 fails to converge when ['x] [^->] 0

When ['x] is large (['x] > 2), both continued fractions converge (CF1
may be slow for really large ['x]). J[sub [mu]], J[sub [mu]+1],
Y[sub [mu]], Y[sub [mu]+1] can be calculated by

[equation bessel13]

where

[equation bessel14]

J[sub [nu]] and Y[sub [mu]] are then calculated using backward
(Miller's algorithm) and forward recurrence respectively.

When ['x] is small (['x] <= 2), CF2 convergence may fail (but CF1
works very well). The solution here is Temme's series:

[equation bessel15]

where

[equation bessel16]

g[sub k][space] and h[sub k][space]
are also computed by recursions (involving gamma functions), but the
formulas are a little complicated, readers are refered to 
N.M. Temme, ['On the numerical evaluation of the ordinary Bessel function
of the second kind], Journal of Computational Physics, vol 21, 343 (1976).
Note Temme's series converge only for |[mu]| <= 1/2.

As the previous case, Y[sub [nu]][space] is calculated from the forward
recurrence, so is Y[sub [nu]+1]. With these two
values and f[sub [nu]], the Wronskian yields J[sub [nu]](x) directly 
without backward recurrence.

[endsect]

[/ 
  Copyright 2006 John Maddock, Paul A. Bristow and Xiaogang Zhang.
  Distributed under the Boost Software License, Version 1.0.
  (See accompanying file LICENSE_1_0.txt or copy at
  http://www.boost.org/LICENSE_1_0.txt).
]