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

[h4 Synopsis]

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

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

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

The functions __cyl_bessel_i and __cyl_bessel_k return the result of the
modified Bessel functions of the first and second kind respectively:

cyl_bessel_i(v, x) = I[sub v](x)

cyl_bessel_k(v, x) = K[sub v](x)

where:

[equation mbessel2]

[equation mbessel3]

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 exponential behaviour of I[sub v].

[graph cyl_bessel_i]

The following graph illustrates the exponential decay of K[sub v].

[graph cyl_bessel_k]

[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 a comparison to the __gsl library.  
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_i
[[Significand Size] [Platform and Compiler] [I[sub v]] ]
[[53] [Win32 / Visual C++ 8.0] [Peak=10 Mean=3.4
      GSL Peak=6000]  ]
[[64] [Red Hat Linux IA64 / G++ 3.4] [Peak=11 Mean=3] ]
[[64] [SUSE Linux AMD64 / G++ 4.1] [Peak=11 Mean=4] ]
[[113] [HP-UX / HP aCC 6] [Peak=15 Mean=4]  ]
]

[table Errors Rates in cyl_bessel_k
[[Significand Size] [Platform and Compiler] [K[sub v]] ]
[[53] [Win32 / Visual C++ 8.0] [Peak=9 Mean=2 

GSL Peak=9] ]
[[64] [Red Hat Linux IA64 / G++ 3.4] [Peak=10 Mean=2] ]
[[64] [SUSE Linux AMD64 / G++ 4.1] [Peak=10 Mean=2] ]
[[113] [HP-UX / HP aCC 6] [Peak=12 Mean=5] ]
]

[h4 Implementation]

The following are handled as special cases first:

When computing I[sub v][space] for ['x < 0], then [nu][space] must be an integer
or a domain error occurs.  If [nu][space] is an integer, then the function is
odd if [nu][space] is odd and even if [nu][space] is even, and we can reflect to
['x > 0].

For I[sub v][space] with v equal to 0, 1 or 0.5 are handled as special cases.

The 0 and 1 cases use minimax rational approximations on
finite and infinite intervals. The coefficients are from:

* J.M. Blair and C.A. Edwards, ['Stable rational minimax approximations
    to the modified Bessel functions I_0(x) and I_1(x)], Atomic Energy of Canada
    Limited Report 4928, Chalk River, 1974.
* S. Moshier, ['Methods and Programs for Mathematical Functions],
    Ellis Horwood Ltd, Chichester, 1989.    

While the 0.5 case is a simple trigonometric function:

I[sub 0.5](x) = sqrt(2 / [pi]x) * sinh(x)

For K[sub v][space] with /v/ an integer, the result is calculated using the
recurrence relation:

[equation mbessel5]

starting from K[sub 0][space] and K[sub 1][space] which are calculated
using rational the approximations above.  These rational approximations are
accurate to around 19 digits, and are therefore only used when T has 
no more than 64 binary digits of precision.

In the general case, we first normalize [nu][space] to \[[^0, [inf]])
with the help of the reflection formulae:

[equation mbessel9]

[equation mbessel10]

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 K[sub [mu]](x) and K[sub [mu]+1](x), and use them to obtain
I[sub [nu]](x) and K[sub [nu]](x).

The algorithm is proposed by Temme in 
N.M. Temme, ['On the numerical evaluation of the modified bessel function
    of the third kind], Journal of Computational Physics, vol 19, 324 (1975), 
which needs two continued fractions as well as the Wronskian:

[equation mbessel11]

[equation mbessel12]

[equation mbessel8]

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]). K[sub [mu]][space] and K[sub [mu]+1][space]
can be calculated by

[equation mbessel13]

where

[equation mbessel14]

['S] is also a series that is summed along with CF2, see 
I.J. Thompson and A.R. Barnett, ['Modified Bessel functions I_v and K_v
    of real order and complex argument to selected accuracy], Computer Physics
    Communications, vol 47, 245 (1987).

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

[equation mbessel15]

where

[equation mbessel16]

f[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 referred to 
N.M. Temme, ['On the numerical evaluation of the modified Bessel function
    of the third kind], Journal of Computational Physics, vol 19, 324 (1975).
Note: Temme's series converge only for |[mu]| <= 1/2.

K[sub [nu]](x) is then calculated from the forward 
recurrence, as is K[sub [nu]+1](x). With these two values and
f[sub [nu]], the Wronskian yields I[sub [nu]](x) directly.

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