jinc.hpp

Go to the documentation of this file.

/** \file jinc.hpp
  * \brief Declares and defines the Jinc and Jinc2 functions
  * \ingroup gen_math_files
  * \author Jared R. Males (jaredmales@gmail.com)
  *
  */
 
//***********************************************************************//
// Copyright 2015, 2016, 2017 Jared R. Males (jaredmales@gmail.com)
//
// This file is part of mxlib.
//
// mxlib is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// mxlib is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with mxlib.  If not, see <http://www.gnu.org/licenses/>.
//***********************************************************************//
 
#ifndef math_func_jinc_hpp
#define math_func_jinc_hpp
 
#include <limits>
#include <cmath>
 
#include "bessel.hpp"
#include "precision.hpp"
 
namespace mx
{
namespace math
{
namespace func 
{
 
/// The Jinc function
/** The Jinc function is defined here as
  * \f[ 
  *     Ji(x) = \frac{J_1(x)}{x} 
  * \f]
  * where \f$ J_1 \f$ is the cylindrical bessel function of the first kind of order 1.
  * 
  * Follows the technique in boost sinc_pi, using the Taylor series for small arguments. If x
  * is smaller than \f$ \epsilon \f$, then it returns 1/2.  If x is larger than \f$ \epsilon \f$ but smaller than \f$ \sqrt{\epsilon} \f$, then
  * this function returns
  * \f[
  *    Ji(x) \approx \frac{1}{2} - \frac{x^2}{16}.
  * \f]
  * 
  * \returns the value of Ji(x)
  * 
  * \tparam T is an floating point type 
  * 
  * \ingroup gen_math_bessel
  */
template<typename T> 
T jinc( const T & x /**< [in] the argument */)
{
   T const taylor_0_bound = std::numeric_limits<T>::epsilon();
   T const taylor_2_bound = root_epsilon<T>();
   
   if( std::fabs(x) > taylor_2_bound)
   {
      return bessel_j<int, T>(1,x)/(x);
   }
   else
   {
      // approximation by taylor series in x at 0
      T  result = static_cast<T>(0.5);
      if (std::fabs(x) >= taylor_0_bound)
      {
         T x2 = x*x;
 
         // approximation by taylor series in x at 0 up to order 2
         result -= x2/static_cast<T>(16);
      }
      
      return result;
   }
            
}
 
extern template
float jinc<float>(const float & x);
 
extern template
double jinc<double>(const double & x);
 
extern template
long double jinc<long double>(const long double & x);
 
#ifdef HASQUAD
extern template
__float128 jinc<__float128>(const __float128 & x);
#endif
 
/// The JincN function
/** The JincN function is defined here as
  * \f[ 
  *     Ji_N(x) = \frac{J_N(x)}{x} 
  * \f]
  * where \f$ J_N \f$ is the cylindrical bessel function of the first kind of order N, \f$ N \ge 1 \f$.
  * 
  * If \f$ N == 1 \f$ this returns \ref mx::math::func::jinc() "jinc(x)".
  * 
  * Otherwise, if x is smaller than \f$ \sqrt{\epsilon} \f$, returns 0.
  * 
  * \returns the value of JiN(x)
  * 
  * \ingroup gen_math_bessel
  */
template<typename T1, typename T2> 
T2 jincN( const T1 & v, ///< [in] the Bessel function order
          const T2 & x  ///< [in] the argument 
        )
{
   if(v == 1) return jinc(x);
 
   T2 const taylor_2_bound = root_epsilon<T2>();
      
   if( std::fabs(x) > taylor_2_bound )
   {
      return bessel_j<T1, T2>(v,x)/(x);
   }
   else
   {
      return static_cast<T2>(0);
   }
}
 
extern template
float jincN<float, float>( const float & v,
                           const float & x
                         );
 
extern template
float jincN<int, float>( const int & v,
                         const float & x
                       );
 
extern template
double jincN<double, double>( const double & v,
                              const double & x
                            );
 
extern template
double jincN<int, double>( const int & v,
                           const double & x
                         );
 
 
extern template
long double jincN<long double, long double>( const long double & v,
                                             const long double & x
                                           );
 
extern template
long double jincN<int, long double>( const int & v,
                                     const long double & x
                                   );
 
#ifdef HASQUAD
extern template
 
__float128 jincN<__float128, __float128>( const __float128 & v,
                                          const __float128 & x
                                        );
 
__float128 jincN<int, __float128>( const int & v,
                                   const __float128 & x
                                 );
#endif
 
} //namespace func 
} //namespace math
} //namespace mx
 
#endif //math_func_jinc_hpp