boost/math/special_functions/hankel.hpp

   1 // Copyright John Maddock 2012.
   2 // Use, modification and distribution are subject to the
   3 // Boost Software License, Version 1.0.
   4 // (See accompanying file LICENSE_1_0.txt
   5 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   6
   7 #ifndef BOOST_MATH_HANKEL_HPP
   8 #define BOOST_MATH_HANKEL_HPP
   9
  10 #include <boost/math/special_functions/bessel.hpp>
  11
  12 namespace boost{ namespace math{
  13
  14 namespace detail{
  15
  16 template <class T, class Policy>
  17 std::complex<T> hankel_imp(T v, T x, const bessel_no_int_tag&, const Policy& pol, int sign)
  18 {
  19    BOOST_MATH_STD_USING
  20    static const char* function = "boost::math::cyl_hankel_1<%1%>(%1%,%1%)";
  21
  22    if(x < 0)
  23    {
  24       bool isint_v = floor(v) == v;
  25       T j, y;
  26       bessel_jy(v, -x, &j, &y, need_j | need_y, pol);
  27       std::complex<T> cx(x), cv(v);
  28       std::complex<T> j_result, y_result;
  29       if(isint_v)
  30       {
  31          int s = (iround(v) & 1) ? -1 : 1;
  32          j_result = j * s;
  33          y_result = T(s) * (y - (2 / constants::pi<T>()) * (log(-x) - log(cx)) * j);
  34       }
  35       else
  36       {
  37          j_result = pow(cx, v) * pow(-cx, -v) * j;
  38          T p1 = pow(-x, v);
  39          std::complex<T> p2 = pow(cx, v);
  40          y_result = p1 * y / p2
  41             + (p2 / p1 - p1 / p2) * j / tan(constants::pi<T>() * v);
  42       }
  43       // multiply y_result by i:
  44       y_result = std::complex<T>(-sign * y_result.imag(), sign * y_result.real());
  45       return j_result + y_result;
  46    }
  47
  48    if(x == 0)
  49    {
  50       if(v == 0)
  51       {
  52          // J is 1, Y is -INF
  53          return std::complex<T>(1, sign * -policies::raise_overflow_error<T>(function, 0, pol));
  54       }
  55       else
  56       {
  57          // At least one of J and Y is complex infinity:
  58          return std::complex<T>(policies::raise_overflow_error<T>(function, 0, pol), sign * policies::raise_overflow_error<T>(function, 0, pol));
  59       }
  60    }
  61
  62    T j, y;
  63    bessel_jy(v, x, &j, &y, need_j | need_y, pol);
  64    return std::complex<T>(j, sign * y);
  65 }
  66
  67 template <class T, class Policy>
  68 std::complex<T> hankel_imp(int v, T x, const bessel_int_tag&, const Policy& pol, int sign);
  69
  70 template <class T, class Policy>
  71 inline std::complex<T> hankel_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol, int sign)
  72 {
  73    BOOST_MATH_STD_USING  // ADL of std names.
  74    int ival = detail::iconv(v, pol);
  75    if(0 == v - ival)
  76    {
  77       return hankel_imp(ival, x, bessel_int_tag(), pol, sign);
  78    }
  79    return hankel_imp(v, x, bessel_no_int_tag(), pol, sign);
  80 }
  81
  82 template <class T, class Policy>
  83 inline std::complex<T> hankel_imp(int v, T x, const bessel_int_tag&, const Policy& pol, int sign)
  84 {
  85    BOOST_MATH_STD_USING
  86    if((std::abs(v) < 200) && (x > 0))
  87       return std::complex<T>(bessel_jn(v, x, pol), sign * bessel_yn(v, x, pol));
  88    return hankel_imp(static_cast<T>(v), x, bessel_no_int_tag(), pol, sign);
  89 }
  90
  91 template <class T, class Policy>
  92 inline std::complex<T> sph_hankel_imp(T v, T x, const Policy& pol, int sign)
  93 {
  94    BOOST_MATH_STD_USING
  95    return constants::root_half_pi<T>() * hankel_imp(v + 0.5f, x, bessel_no_int_tag(), pol, sign) / sqrt(std::complex<T>(x));
  96 }
  97
  98 } // namespace detail
  99
 100 template <class T1, class T2, class Policy>
 101 inline std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> cyl_hankel_1(T1 v, T2 x, const Policy& pol)
 102 {
 103    BOOST_FPU_EXCEPTION_GUARD
 104    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
 105    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
 106    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 107    return policies::checked_narrowing_cast<std::complex<result_type>, Policy>(detail::hankel_imp<value_type>(v, static_cast<value_type>(x), tag_type(), pol, 1), "boost::math::cyl_hankel_1<%1%>(%1%,%1%)");
 108 }
 109
 110 template <class T1, class T2>
 111 inline std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> cyl_hankel_1(T1 v, T2 x)
 112 {
 113    return cyl_hankel_1(v, x, policies::policy<>());
 114 }
 115
 116 template <class T1, class T2, class Policy>
 117 inline std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> cyl_hankel_2(T1 v, T2 x, const Policy& pol)
 118 {
 119    BOOST_FPU_EXCEPTION_GUARD
 120    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
 121    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
 122    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 123    return policies::checked_narrowing_cast<std::complex<result_type>, Policy>(detail::hankel_imp<value_type>(v, static_cast<value_type>(x), tag_type(), pol, -1), "boost::math::cyl_hankel_1<%1%>(%1%,%1%)");
 124 }
 125
 126 template <class T1, class T2>
 127 inline std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> cyl_hankel_2(T1 v, T2 x)
 128 {
 129    return cyl_hankel_2(v, x, policies::policy<>());
 130 }
 131
 132 template <class T1, class T2, class Policy>
 133 inline std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> sph_hankel_1(T1 v, T2 x, const Policy&)
 134 {
 135    BOOST_FPU_EXCEPTION_GUARD
 136    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
 137    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 138    typedef typename policies::normalise<
 139       Policy,
 140       policies::promote_float<false>,
 141       policies::promote_double<false>,
 142       policies::discrete_quantile<>,
 143       policies::assert_undefined<> >::type forwarding_policy;
 144
 145    return policies::checked_narrowing_cast<std::complex<result_type>, Policy>(detail::sph_hankel_imp<value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy(), 1), "boost::math::sph_hankel_1<%1%>(%1%,%1%)");
 146 }
 147
 148 template <class T1, class T2>
 149 inline std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> sph_hankel_1(T1 v, T2 x)
 150 {
 151    return sph_hankel_1(v, x, policies::policy<>());
 152 }
 153
 154 template <class T1, class T2, class Policy>
 155 inline std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> sph_hankel_2(T1 v, T2 x, const Policy&)
 156 {
 157    BOOST_FPU_EXCEPTION_GUARD
 158    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
 159    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 160    typedef typename policies::normalise<
 161       Policy,
 162       policies::promote_float<false>,
 163       policies::promote_double<false>,
 164       policies::discrete_quantile<>,
 165       policies::assert_undefined<> >::type forwarding_policy;
 166
 167    return policies::checked_narrowing_cast<std::complex<result_type>, Policy>(detail::sph_hankel_imp<value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy(), -1), "boost::math::sph_hankel_1<%1%>(%1%,%1%)");
 168 }
 169
 170 template <class T1, class T2>
 171 inline std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> sph_hankel_2(T1 v, T2 x)
 172 {
 173    return sph_hankel_2(v, x, policies::policy<>());
 174 }
 175
 176 }} // namespaces
 177
 178 #endif // BOOST_MATH_HANKEL_HPP
 179