]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Copyright (c) 2013 Anton Bikineev |
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 | // | |
7 | // This is a partial header, do not include on it's own!!! | |
8 | // | |
9 | // Contains asymptotic expansions for derivatives of Bessel J(v,x) and Y(v,x) | |
10 | // functions, as x -> INF. | |
11 | #ifndef BOOST_MATH_SF_DETAIL_BESSEL_JY_DERIVATIVES_ASYM_HPP | |
12 | #define BOOST_MATH_SF_DETAIL_BESSEL_JY_DERIVATIVES_ASYM_HPP | |
13 | ||
14 | #ifdef _MSC_VER | |
15 | #pragma once | |
16 | #endif | |
17 | ||
18 | namespace boost{ namespace math{ namespace detail{ | |
19 | ||
20 | template <class T> | |
21 | inline T asymptotic_bessel_derivative_amplitude(T v, T x) | |
22 | { | |
23 | // Calculate the amplitude for J'(v,x) and I'(v,x) | |
24 | // for large x: see A&S 9.2.30. | |
25 | BOOST_MATH_STD_USING | |
26 | T s = 1; | |
27 | const T mu = 4 * v * v; | |
28 | T txq = 2 * x; | |
29 | txq *= txq; | |
30 | ||
31 | s -= (mu - 3) / (2 * txq); | |
32 | s -= ((mu - 1) * (mu - 45)) / (txq * txq * 8); | |
33 | ||
34 | return sqrt(s * 2 / (boost::math::constants::pi<T>() * x)); | |
35 | } | |
36 | ||
37 | template <class T> | |
38 | inline T asymptotic_bessel_derivative_phase_mx(T v, T x) | |
39 | { | |
40 | // Calculate the phase of J'(v, x) and Y'(v, x) for large x. | |
41 | // See A&S 9.2.31. | |
42 | // Note that the result returned is the phase less (x - PI(v/2 - 1/4)) | |
43 | // which we'll factor in later when we calculate the sines/cosines of the result: | |
44 | const T mu = 4 * v * v; | |
45 | const T mu2 = mu * mu; | |
46 | const T mu3 = mu2 * mu; | |
47 | T denom = 4 * x; | |
48 | T denom_mult = denom * denom; | |
49 | ||
50 | T s = 0; | |
51 | s += (mu + 3) / (2 * denom); | |
52 | denom *= denom_mult; | |
53 | s += (mu2 + (46 * mu) - 63) / (6 * denom); | |
54 | denom *= denom_mult; | |
55 | s += (mu3 + (185 * mu2) - (2053 * mu) + 1899) / (5 * denom); | |
56 | return s; | |
57 | } | |
58 | ||
59 | template <class T> | |
60 | inline T asymptotic_bessel_y_derivative_large_x_2(T v, T x) | |
61 | { | |
62 | // See A&S 9.2.20. | |
63 | BOOST_MATH_STD_USING | |
64 | // Get the phase and amplitude: | |
65 | const T ampl = asymptotic_bessel_derivative_amplitude(v, x); | |
66 | const T phase = asymptotic_bessel_derivative_phase_mx(v, x); | |
67 | BOOST_MATH_INSTRUMENT_VARIABLE(ampl); | |
68 | BOOST_MATH_INSTRUMENT_VARIABLE(phase); | |
69 | // | |
70 | // Calculate the sine of the phase, using | |
71 | // sine/cosine addition rules to factor in | |
72 | // the x - PI(v/2 - 1/4) term not added to the | |
73 | // phase when we calculated it. | |
74 | // | |
75 | const T cx = cos(x); | |
76 | const T sx = sin(x); | |
77 | const T vd2shifted = (v / 2) - 0.25f; | |
78 | const T ci = cos_pi(vd2shifted); | |
79 | const T si = sin_pi(vd2shifted); | |
80 | const T sin_phase = sin(phase) * (cx * ci + sx * si) + cos(phase) * (sx * ci - cx * si); | |
81 | BOOST_MATH_INSTRUMENT_CODE(sin(phase)); | |
82 | BOOST_MATH_INSTRUMENT_CODE(cos(x)); | |
83 | BOOST_MATH_INSTRUMENT_CODE(cos(phase)); | |
84 | BOOST_MATH_INSTRUMENT_CODE(sin(x)); | |
85 | return sin_phase * ampl; | |
86 | } | |
87 | ||
88 | template <class T> | |
89 | inline T asymptotic_bessel_j_derivative_large_x_2(T v, T x) | |
90 | { | |
91 | // See A&S 9.2.20. | |
92 | BOOST_MATH_STD_USING | |
93 | // Get the phase and amplitude: | |
94 | const T ampl = asymptotic_bessel_derivative_amplitude(v, x); | |
95 | const T phase = asymptotic_bessel_derivative_phase_mx(v, x); | |
96 | BOOST_MATH_INSTRUMENT_VARIABLE(ampl); | |
97 | BOOST_MATH_INSTRUMENT_VARIABLE(phase); | |
98 | // | |
99 | // Calculate the sine of the phase, using | |
100 | // sine/cosine addition rules to factor in | |
101 | // the x - PI(v/2 - 1/4) term not added to the | |
102 | // phase when we calculated it. | |
103 | // | |
104 | BOOST_MATH_INSTRUMENT_CODE(cos(phase)); | |
105 | BOOST_MATH_INSTRUMENT_CODE(cos(x)); | |
106 | BOOST_MATH_INSTRUMENT_CODE(sin(phase)); | |
107 | BOOST_MATH_INSTRUMENT_CODE(sin(x)); | |
108 | const T cx = cos(x); | |
109 | const T sx = sin(x); | |
110 | const T vd2shifted = (v / 2) - 0.25f; | |
111 | const T ci = cos_pi(vd2shifted); | |
112 | const T si = sin_pi(vd2shifted); | |
113 | const T sin_phase = cos(phase) * (cx * ci + sx * si) - sin(phase) * (sx * ci - cx * si); | |
114 | BOOST_MATH_INSTRUMENT_VARIABLE(sin_phase); | |
115 | return sin_phase * ampl; | |
116 | } | |
117 | ||
118 | template <class T> | |
119 | inline bool asymptotic_bessel_derivative_large_x_limit(const T& v, const T& x) | |
120 | { | |
121 | BOOST_MATH_STD_USING | |
122 | // | |
123 | // This function is the copy of math::asymptotic_bessel_large_x_limit | |
124 | // It means that we use the same rules for determining how x is large | |
125 | // compared to v. | |
126 | // | |
127 | // Determines if x is large enough compared to v to take the asymptotic | |
128 | // forms above. From A&S 9.2.28 we require: | |
129 | // v < x * eps^1/8 | |
130 | // and from A&S 9.2.29 we require: | |
131 | // v^12/10 < 1.5 * x * eps^1/10 | |
132 | // using the former seems to work OK in practice with broadly similar | |
133 | // error rates either side of the divide for v < 10000. | |
134 | // At double precision eps^1/8 ~= 0.01. | |
135 | // | |
136 | return (std::max)(T(fabs(v)), T(1)) < x * sqrt(boost::math::tools::forth_root_epsilon<T>()); | |
137 | } | |
138 | ||
139 | }}} // namespaces | |
140 | ||
141 | #endif // BOOST_MATH_SF_DETAIL_BESSEL_JY_DERIVATIVES_ASYM_HPP |