math/special_functions/hankel.hpp

0001 // Copyright John Maddock 2012.
0002 // Use, modification and distribution are subject to the
0003 // Boost Software License, Version 1.0.
0004 // (See accompanying file LICENSE_1_0.txt
0005 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_MATH_HANKEL_HPP
0008 #define BOOST_MATH_HANKEL_HPP
0009
0010 #include <boost/math/special_functions/math_fwd.hpp>
0011 #include <boost/math/special_functions/bessel.hpp>
0012
0013 namespace boost{ namespace math{
0014
0015 namespace detail{
0016
0017 template <class T, class Policy>
0018 std::complex<T> hankel_imp(T v, T x, const bessel_no_int_tag&, const Policy& pol, int sign)
0019 {
0020    BOOST_MATH_STD_USING
0021    static const char* function = "boost::math::cyl_hankel_1<%1%>(%1%,%1%)";
0022
0023    if(x < 0)
0024    {
0025       bool isint_v = floor(v) == v;
0026       T j, y;
0027       bessel_jy(v, -x, &j, &y, need_j | need_y, pol);
0028       std::complex<T> cx(x), cv(v);
0029       std::complex<T> j_result, y_result;
0030       if(isint_v)
0031       {
0032          int s = (iround(v) & 1) ? -1 : 1;
0033          j_result = j * s;
0034          y_result = T(s) * (y - (2 / constants::pi<T>()) * (log(-x) - log(cx)) * j);
0035       }
0036       else
0037       {
0038          j_result = pow(cx, v) * pow(-cx, -v) * j;
0039          T p1 = pow(-x, v);
0040          std::complex<T> p2 = pow(cx, v);
0041          y_result = p1 * y / p2
0042             + (p2 / p1 - p1 / p2) * j / tan(constants::pi<T>() * v);
0043       }
0044       // multiply y_result by i:
0045       y_result = std::complex<T>(-sign * y_result.imag(), sign * y_result.real());
0046       return j_result + y_result;
0047    }
0048
0049    if(x == 0)
0050    {
0051       if(v == 0)
0052       {
0053          // J is 1, Y is -INF
0054          return std::complex<T>(1, sign * -policies::raise_overflow_error<T>(function, nullptr, pol));
0055       }
0056       else
0057       {
0058          // At least one of J and Y is complex infinity:
0059          return std::complex<T>(policies::raise_overflow_error<T>(function, nullptr, pol), sign * policies::raise_overflow_error<T>(function, nullptr, pol));
0060       }
0061    }
0062
0063    T j, y;
0064    bessel_jy(v, x, &j, &y, need_j | need_y, pol);
0065    return std::complex<T>(j, sign * y);
0066 }
0067
0068 template <class T, class Policy>
0069 std::complex<T> hankel_imp(int v, T x, const bessel_int_tag&, const Policy& pol, int sign);
0070
0071 template <class T, class Policy>
0072 inline std::complex<T> hankel_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol, int sign)
0073 {
0074    BOOST_MATH_STD_USING  // ADL of std names.
0075    int ival = detail::iconv(v, pol);
0076    if(0 == v - ival)
0077    {
0078       return hankel_imp(ival, x, bessel_int_tag(), pol, sign);
0079    }
0080    return hankel_imp(v, x, bessel_no_int_tag(), pol, sign);
0081 }
0082
0083 template <class T, class Policy>
0084 inline std::complex<T> hankel_imp(int v, T x, const bessel_int_tag&, const Policy& pol, int sign)
0085 {
0086    BOOST_MATH_STD_USING
0087    if((std::abs(v) < 200) && (x > 0))
0088       return std::complex<T>(bessel_jn(v, x, pol), sign * bessel_yn(v, x, pol));
0089    return hankel_imp(static_cast<T>(v), x, bessel_no_int_tag(), pol, sign);
0090 }
0091
0092 template <class T, class Policy>
0093 inline std::complex<T> sph_hankel_imp(T v, T x, const Policy& pol, int sign)
0094 {
0095    BOOST_MATH_STD_USING
0096    return constants::root_half_pi<T>() * hankel_imp(v + 0.5f, x, bessel_no_int_tag(), pol, sign) / sqrt(std::complex<T>(x));
0097 }
0098
0099 } // namespace detail
0100
0101 template <class T1, class T2, class Policy>
0102 inline std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> cyl_hankel_1(T1 v, T2 x, const Policy& pol)
0103 {
0104    BOOST_FPU_EXCEPTION_GUARD
0105    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
0106    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
0107    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0108    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%)");
0109 }
0110
0111 template <class T1, class T2>
0112 inline std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> cyl_hankel_1(T1 v, T2 x)
0113 {
0114    return cyl_hankel_1(v, x, policies::policy<>());
0115 }
0116
0117 template <class T1, class T2, class Policy>
0118 inline std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> cyl_hankel_2(T1 v, T2 x, const Policy& pol)
0119 {
0120    BOOST_FPU_EXCEPTION_GUARD
0121    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
0122    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
0123    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0124    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%)");
0125 }
0126
0127 template <class T1, class T2>
0128 inline std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> cyl_hankel_2(T1 v, T2 x)
0129 {
0130    return cyl_hankel_2(v, x, policies::policy<>());
0131 }
0132
0133 template <class T1, class T2, class Policy>
0134 inline std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> sph_hankel_1(T1 v, T2 x, const Policy&)
0135 {
0136    BOOST_FPU_EXCEPTION_GUARD
0137    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
0138    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0139    typedef typename policies::normalise<
0140       Policy,
0141       policies::promote_float<false>,
0142       policies::promote_double<false>,
0143       policies::discrete_quantile<>,
0144       policies::assert_undefined<> >::type forwarding_policy;
0145
0146    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%)");
0147 }
0148
0149 template <class T1, class T2>
0150 inline std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> sph_hankel_1(T1 v, T2 x)
0151 {
0152    return sph_hankel_1(v, x, policies::policy<>());
0153 }
0154
0155 template <class T1, class T2, class Policy>
0156 inline std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> sph_hankel_2(T1 v, T2 x, const Policy&)
0157 {
0158    BOOST_FPU_EXCEPTION_GUARD
0159    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
0160    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0161    typedef typename policies::normalise<
0162       Policy,
0163       policies::promote_float<false>,
0164       policies::promote_double<false>,
0165       policies::discrete_quantile<>,
0166       policies::assert_undefined<> >::type forwarding_policy;
0167
0168    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%)");
0169 }
0170
0171 template <class T1, class T2>
0172 inline std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> sph_hankel_2(T1 v, T2 x)
0173 {
0174    return sph_hankel_2(v, x, policies::policy<>());
0175 }
0176
0177 }} // namespaces
0178
0179 #endif // BOOST_MATH_HANKEL_HPP
0180