1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
|
/* Copyright (c) 2008-2022 the MRtrix3 contributors.
*
* This Source Code Form is subject to the terms of the Mozilla Public
* License, v. 2.0. If a copy of the MPL was not distributed with this
* file, You can obtain one at http://mozilla.org/MPL/2.0/.
*
* Covered Software is provided under this License on an "as is"
* basis, without warranty of any kind, either expressed, implied, or
* statutory, including, without limitation, warranties that the
* Covered Software is free of defects, merchantable, fit for a
* particular purpose or non-infringing.
* See the Mozilla Public License v. 2.0 for more details.
*
* For more details, see http://www.mrtrix.org/.
*/
#ifndef __math_bessel_h__
#define __math_bessel_h__
#include <limits>
#include "math/chebyshev.h"
namespace MR
{
namespace Math
{
namespace Bessel
{
extern const double coef_aI0[];
extern const double coef_bI0[];
extern const double coef_cI0[];
extern const double coef_aI1[];
extern const double coef_bI1[];
extern const double coef_cI1[];
//* Compute the scaled regular modified cylindrical Bessel function of zeroth order exp(-|x|) I_0(x). */
/** Implementation based on the GSL (http://www.gnu.org/software/gsl/) */
template <typename T> inline T I0_scaled (const T x)
{
assert (x >= 0.0);
if (x*x < 4.0*std::numeric_limits<T>::epsilon()) return (1.0-x);
if (x <= 3.0) return (exp (-x) * (2.75 + Chebyshev::eval (coef_aI0, 11, -1.0, 1.0, x*x/4.5-1.0)));
if (x <= 8.0) return ( (0.375 + Chebyshev::eval (coef_bI0, (sizeof (T) >4?20:13), -1.0, 1.0, (48.0/x-11.0) /5.0)) /sqrt (x));
return ( (0.375 + Chebyshev::eval (coef_cI0, (sizeof (T) >4?21:11), -1.0, 1.0, 16.0/x-1.0)) /sqrt (x));
}
//* Compute the scaled regular modified cylindrical Bessel function of first order exp(-|x|) I_1(x). */
/** Implementation based on the GSL (http://www.gnu.org/software/gsl/) */
template <typename T> inline T I1_scaled (const T x)
{
assert (x >= 0.0);
if (x == 0.0) return (0.0);
if (x*x < 8.0*std::numeric_limits<T>::epsilon()) return (0.5*x);
if (x <= 3.0) return (x * exp (-x) * (0.875 + Chebyshev::eval (coef_aI1, 10, -1.0, 1.0, x*x/4.5-1.0)));
if (x <= 8.0) {
T b = (0.375 + Chebyshev::eval (coef_bI1, (sizeof (T) >4?20:11), -1.0, 1.0, (48.0/x-11.0) /5.0)) / sqrt (x);
return (x > 0.0 ? b : -b);
}
T b = (0.375 + Chebyshev::eval (coef_cI1, (sizeof (T) >4?21:9), -1.0, 1.0, 16.0/x-1.0)) / sqrt (x);
return (x > 0.0 ? b : -b);
}
}
}
}
#endif
|
