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