Line data Source code
1 : /** \file jinc.hpp
2 : * \brief Declares and defines the Jinc and Jinc2 functions
3 : * \ingroup gen_math_files
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2015, 2016, 2017 Jared R. Males (jaredmales@gmail.com)
10 : //
11 : // This file is part of mxlib.
12 : //
13 : // mxlib is free software: you can redistribute it and/or modify
14 : // it under the terms of the GNU General Public License as published by
15 : // the Free Software Foundation, either version 3 of the License, or
16 : // (at your option) any later version.
17 : //
18 : // mxlib is distributed in the hope that it will be useful,
19 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
20 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 : // GNU General Public License for more details.
22 : //
23 : // You should have received a copy of the GNU General Public License
24 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
25 : //***********************************************************************//
26 :
27 : #ifndef math_func_jinc_hpp
28 : #define math_func_jinc_hpp
29 :
30 : #include <limits>
31 : #include <cmath>
32 :
33 : #include "bessel.hpp"
34 : #include "precision.hpp"
35 :
36 : namespace mx
37 : {
38 : namespace math
39 : {
40 : namespace func
41 : {
42 :
43 : /// The Jinc function
44 : /** The Jinc function is defined here as
45 : * \f[
46 : * Ji(x) = \frac{J_1(x)}{x}
47 : * \f]
48 : * where \f$ J_1 \f$ is the cylindrical bessel function of the first kind of order 1.
49 : *
50 : * Follows the technique in boost sinc_pi, using the Taylor series for small arguments. If x
51 : * is smaller than \f$ \epsilon \f$, then it returns 1/2. If x is larger than \f$ \epsilon \f$ but smaller than \f$
52 : * \sqrt{\epsilon} \f$, then this function returns \f[ Ji(x) \approx \frac{1}{2} - \frac{x^2}{16}. \f]
53 : *
54 : * \returns the value of Ji(x)
55 : *
56 : * \tparam T is an floating point type
57 : *
58 : * \ingroup gen_math_bessel
59 : */
60 : template <typename T>
61 330422 : T jinc( const T &x /**< [in] the argument */ )
62 : {
63 330422 : T const taylor_0_bound = std::numeric_limits<T>::epsilon();
64 330422 : T const taylor_2_bound = root_epsilon<T>();
65 :
66 330422 : if( std::fabs( x ) > taylor_2_bound )
67 : {
68 330420 : return bessel_j<int, T>( 1, x ) / ( x );
69 : }
70 : else
71 : {
72 : // approximation by taylor series in x at 0
73 2 : T result = static_cast<T>( 0.5 );
74 2 : if( std::fabs( x ) >= taylor_0_bound )
75 : {
76 0 : T x2 = x * x;
77 :
78 : // approximation by taylor series in x at 0 up to order 2
79 0 : result -= x2 / static_cast<T>( 16 );
80 : }
81 :
82 2 : return result;
83 : }
84 : }
85 :
86 : extern template float jinc<float>( const float &x );
87 :
88 : extern template double jinc<double>( const double &x );
89 :
90 : extern template long double jinc<long double>( const long double &x );
91 :
92 : #ifdef HASQUAD
93 : extern template __float128 jinc<__float128>( const __float128 &x );
94 : #endif
95 :
96 : /// The JincN function
97 : /** The JincN function is defined here as
98 : * \f[
99 : * Ji_N(x) = \frac{J_N(x)}{x}
100 : * \f]
101 : * where \f$ J_N \f$ is the cylindrical bessel function of the first kind of order N, \f$ N \ge 1 \f$.
102 : *
103 : * If \f$ N == 1 \f$ this returns \ref mx::math::func::jinc() "jinc(x)".
104 : *
105 : * Otherwise, if x is smaller than \f$ \sqrt{\epsilon} \f$, returns 0.
106 : *
107 : * \returns the value of JiN(x)
108 : *
109 : * \ingroup gen_math_bessel
110 : */
111 : template <typename T1, typename T2>
112 0 : T2 jincN( const T1 &v, ///< [in] the Bessel function order
113 : const T2 &x ///< [in] the argument
114 : )
115 : {
116 0 : if( v == 1 )
117 0 : return jinc( x );
118 :
119 0 : T2 const taylor_2_bound = root_epsilon<T2>();
120 :
121 0 : if( std::fabs( x ) > taylor_2_bound )
122 : {
123 0 : return bessel_j<T1, T2>( v, x ) / ( x );
124 : }
125 : else
126 : {
127 0 : return static_cast<T2>( 0 );
128 : }
129 : }
130 :
131 : extern template float jincN<float, float>( const float &v, const float &x );
132 :
133 : extern template float jincN<int, float>( const int &v, const float &x );
134 :
135 : extern template double jincN<double, double>( const double &v, const double &x );
136 :
137 : extern template double jincN<int, double>( const int &v, const double &x );
138 :
139 : extern template long double jincN<long double, long double>( const long double &v, const long double &x );
140 :
141 : extern template long double jincN<int, long double>( const int &v, const long double &x );
142 :
143 : #ifdef HASQUAD
144 : extern template
145 :
146 : __float128
147 : jincN<__float128, __float128>( const __float128 &v, const __float128 &x );
148 :
149 : __float128 jincN<int, __float128>( const int &v, const __float128 &x );
150 : #endif
151 :
152 : } // namespace func
153 : } // namespace math
154 : } // namespace mx
155 :
156 : #endif // math_func_jinc_hpp
|
