]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Copyright (c) 2015 John Maddock |
2 | // Use, modification and distribution are subject to the | |
3 | // Boost Software License, Version 1.0. (See accompanying file | |
4 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
5 | ||
6 | #ifndef BOOST_MATH_ELLINT_HL_HPP | |
7 | #define BOOST_MATH_ELLINT_HL_HPP | |
8 | ||
9 | #ifdef _MSC_VER | |
10 | #pragma once | |
11 | #endif | |
12 | ||
13 | #include <boost/math/special_functions/math_fwd.hpp> | |
14 | #include <boost/math/special_functions/ellint_rj.hpp> | |
15 | #include <boost/math/special_functions/ellint_rj.hpp> | |
16 | #include <boost/math/special_functions/ellint_1.hpp> | |
17 | #include <boost/math/special_functions/jacobi_zeta.hpp> | |
18 | #include <boost/math/constants/constants.hpp> | |
19 | #include <boost/math/policies/error_handling.hpp> | |
20 | #include <boost/math/tools/workaround.hpp> | |
21 | ||
22 | // Elliptic integral the Jacobi Zeta function. | |
23 | ||
24 | namespace boost { namespace math { | |
25 | ||
26 | namespace detail{ | |
27 | ||
28 | // Elliptic integral - Jacobi Zeta | |
29 | template <typename T, typename Policy> | |
30 | T heuman_lambda_imp(T phi, T k, const Policy& pol) | |
31 | { | |
32 | BOOST_MATH_STD_USING | |
33 | using namespace boost::math::tools; | |
34 | using namespace boost::math::constants; | |
35 | ||
36 | const char* function = "boost::math::heuman_lambda<%1%>(%1%, %1%)"; | |
37 | ||
38 | if(fabs(k) > 1) | |
39 | return policies::raise_domain_error<T>(function, "We require |k| <= 1 but got k = %1%", k, pol); | |
40 | ||
41 | T result; | |
42 | T sinp = sin(phi); | |
43 | T cosp = cos(phi); | |
44 | T s2 = sinp * sinp; | |
45 | T k2 = k * k; | |
46 | T kp = 1 - k2; | |
47 | T delta = sqrt(1 - (kp * s2)); | |
48 | if(fabs(phi) <= constants::half_pi<T>()) | |
49 | { | |
50 | result = kp * sinp * cosp / (delta * constants::half_pi<T>()); | |
51 | result *= ellint_rf_imp(T(0), kp, T(1), pol) + k2 * ellint_rj(T(0), kp, T(1), T(1 - k2 / (delta * delta)), pol) / (3 * delta * delta); | |
52 | } | |
53 | else | |
54 | { | |
55 | T rkp = sqrt(kp); | |
56 | T ratio; | |
57 | if(rkp == 1) | |
58 | { | |
59 | return policies::raise_domain_error<T>(function, "When 1-k^2 == 1 then phi must be < Pi/2, but got phi = %1%", phi, pol); | |
60 | } | |
61 | else | |
62 | ratio = ellint_f_imp(phi, rkp, pol) / ellint_k_imp(rkp, pol); | |
63 | result = ratio + ellint_k_imp(k, pol) * jacobi_zeta_imp(phi, rkp, pol) / constants::half_pi<T>(); | |
64 | } | |
65 | return result; | |
66 | } | |
67 | ||
68 | } // detail | |
69 | ||
70 | template <class T1, class T2, class Policy> | |
71 | inline typename tools::promote_args<T1, T2>::type heuman_lambda(T1 k, T2 phi, const Policy& pol) | |
72 | { | |
73 | typedef typename tools::promote_args<T1, T2>::type result_type; | |
74 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
75 | return policies::checked_narrowing_cast<result_type, Policy>(detail::heuman_lambda_imp(static_cast<value_type>(phi), static_cast<value_type>(k), pol), "boost::math::heuman_lambda<%1%>(%1%,%1%)"); | |
76 | } | |
77 | ||
78 | template <class T1, class T2> | |
79 | inline typename tools::promote_args<T1, T2>::type heuman_lambda(T1 k, T2 phi) | |
80 | { | |
81 | return boost::math::heuman_lambda(k, phi, policies::policy<>()); | |
82 | } | |
83 | ||
84 | }} // namespaces | |
85 | ||
86 | #endif // BOOST_MATH_ELLINT_D_HPP | |
87 |