math/special_functions/hankel.hpp

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