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