]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | |
2 | // (C) Copyright John Maddock 2006. | |
3 | // Use, modification and distribution are subject to the | |
4 | // Boost Software License, Version 1.0. (See accompanying file | |
5 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
6 | ||
7 | #ifndef BOOST_MATH_SPECIAL_SPHERICAL_HARMONIC_HPP | |
8 | #define BOOST_MATH_SPECIAL_SPHERICAL_HARMONIC_HPP | |
9 | ||
10 | #ifdef _MSC_VER | |
11 | #pragma once | |
12 | #endif | |
13 | ||
14 | #include <boost/math/special_functions/math_fwd.hpp> | |
15 | #include <boost/math/special_functions/legendre.hpp> | |
16 | #include <boost/math/tools/workaround.hpp> | |
17 | #include <complex> | |
18 | ||
19 | namespace boost{ | |
20 | namespace math{ | |
21 | ||
22 | namespace detail{ | |
23 | ||
24 | // | |
25 | // Calculates the prefix term that's common to the real | |
26 | // and imaginary parts. Does *not* fix up the sign of the result | |
27 | // though. | |
28 | // | |
29 | template <class T, class Policy> | |
30 | inline T spherical_harmonic_prefix(unsigned n, unsigned m, T theta, const Policy& pol) | |
31 | { | |
32 | BOOST_MATH_STD_USING | |
33 | ||
34 | if(m > n) | |
35 | return 0; | |
36 | ||
37 | T sin_theta = sin(theta); | |
38 | T x = cos(theta); | |
39 | ||
40 | T leg = detail::legendre_p_imp(n, m, x, static_cast<T>(pow(fabs(sin_theta), T(m))), pol); | |
41 | ||
42 | T prefix = boost::math::tgamma_delta_ratio(static_cast<T>(n - m + 1), static_cast<T>(2 * m), pol); | |
43 | prefix *= (2 * n + 1) / (4 * constants::pi<T>()); | |
44 | prefix = sqrt(prefix); | |
45 | return prefix * leg; | |
46 | } | |
47 | // | |
48 | // Real Part: | |
49 | // | |
50 | template <class T, class Policy> | |
51 | T spherical_harmonic_r(unsigned n, int m, T theta, T phi, const Policy& pol) | |
52 | { | |
53 | BOOST_MATH_STD_USING // ADL of std functions | |
54 | ||
55 | bool sign = false; | |
56 | if(m < 0) | |
57 | { | |
58 | // Reflect and adjust sign if m < 0: | |
59 | sign = m&1; | |
60 | m = abs(m); | |
61 | } | |
62 | if(m&1) | |
63 | { | |
64 | // Check phase if theta is outside [0, PI]: | |
65 | T mod = boost::math::tools::fmod_workaround(theta, T(2 * constants::pi<T>())); | |
66 | if(mod < 0) | |
67 | mod += 2 * constants::pi<T>(); | |
68 | if(mod > constants::pi<T>()) | |
69 | sign = !sign; | |
70 | } | |
71 | // Get the value and adjust sign as required: | |
72 | T prefix = spherical_harmonic_prefix(n, m, theta, pol); | |
73 | prefix *= cos(m * phi); | |
74 | return sign ? T(-prefix) : prefix; | |
75 | } | |
76 | ||
77 | template <class T, class Policy> | |
78 | T spherical_harmonic_i(unsigned n, int m, T theta, T phi, const Policy& pol) | |
79 | { | |
80 | BOOST_MATH_STD_USING // ADL of std functions | |
81 | ||
82 | bool sign = false; | |
83 | if(m < 0) | |
84 | { | |
85 | // Reflect and adjust sign if m < 0: | |
86 | sign = !(m&1); | |
87 | m = abs(m); | |
88 | } | |
89 | if(m&1) | |
90 | { | |
91 | // Check phase if theta is outside [0, PI]: | |
92 | T mod = boost::math::tools::fmod_workaround(theta, T(2 * constants::pi<T>())); | |
93 | if(mod < 0) | |
94 | mod += 2 * constants::pi<T>(); | |
95 | if(mod > constants::pi<T>()) | |
96 | sign = !sign; | |
97 | } | |
98 | // Get the value and adjust sign as required: | |
99 | T prefix = spherical_harmonic_prefix(n, m, theta, pol); | |
100 | prefix *= sin(m * phi); | |
101 | return sign ? T(-prefix) : prefix; | |
102 | } | |
103 | ||
104 | template <class T, class U, class Policy> | |
105 | std::complex<T> spherical_harmonic(unsigned n, int m, U theta, U phi, const Policy& pol) | |
106 | { | |
107 | BOOST_MATH_STD_USING | |
108 | // | |
109 | // Sort out the signs: | |
110 | // | |
111 | bool r_sign = false; | |
112 | bool i_sign = false; | |
113 | if(m < 0) | |
114 | { | |
115 | // Reflect and adjust sign if m < 0: | |
116 | r_sign = m&1; | |
117 | i_sign = !(m&1); | |
118 | m = abs(m); | |
119 | } | |
120 | if(m&1) | |
121 | { | |
122 | // Check phase if theta is outside [0, PI]: | |
123 | U mod = boost::math::tools::fmod_workaround(theta, U(2 * constants::pi<U>())); | |
124 | if(mod < 0) | |
125 | mod += 2 * constants::pi<U>(); | |
126 | if(mod > constants::pi<U>()) | |
127 | { | |
128 | r_sign = !r_sign; | |
129 | i_sign = !i_sign; | |
130 | } | |
131 | } | |
132 | // | |
133 | // Calculate the value: | |
134 | // | |
135 | U prefix = spherical_harmonic_prefix(n, m, theta, pol); | |
136 | U r = prefix * cos(m * phi); | |
137 | U i = prefix * sin(m * phi); | |
138 | // | |
139 | // Add in the signs: | |
140 | // | |
141 | if(r_sign) | |
142 | r = -r; | |
143 | if(i_sign) | |
144 | i = -i; | |
145 | static const char* function = "boost::math::spherical_harmonic<%1%>(int, int, %1%, %1%)"; | |
146 | return std::complex<T>(policies::checked_narrowing_cast<T, Policy>(r, function), policies::checked_narrowing_cast<T, Policy>(i, function)); | |
147 | } | |
148 | ||
149 | } // namespace detail | |
150 | ||
151 | template <class T1, class T2, class Policy> | |
152 | inline std::complex<typename tools::promote_args<T1, T2>::type> | |
153 | spherical_harmonic(unsigned n, int m, T1 theta, T2 phi, const Policy& pol) | |
154 | { | |
155 | typedef typename tools::promote_args<T1, T2>::type result_type; | |
156 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
157 | return detail::spherical_harmonic<result_type, value_type>(n, m, static_cast<value_type>(theta), static_cast<value_type>(phi), pol); | |
158 | } | |
159 | ||
160 | template <class T1, class T2> | |
161 | inline std::complex<typename tools::promote_args<T1, T2>::type> | |
162 | spherical_harmonic(unsigned n, int m, T1 theta, T2 phi) | |
163 | { | |
164 | return boost::math::spherical_harmonic(n, m, theta, phi, policies::policy<>()); | |
165 | } | |
166 | ||
167 | template <class T1, class T2, class Policy> | |
168 | inline typename tools::promote_args<T1, T2>::type | |
169 | spherical_harmonic_r(unsigned n, int m, T1 theta, T2 phi, const Policy& pol) | |
170 | { | |
171 | typedef typename tools::promote_args<T1, T2>::type result_type; | |
172 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
173 | return policies::checked_narrowing_cast<result_type, Policy>(detail::spherical_harmonic_r(n, m, static_cast<value_type>(theta), static_cast<value_type>(phi), pol), "bost::math::spherical_harmonic_r<%1%>(unsigned, int, %1%, %1%)"); | |
174 | } | |
175 | ||
176 | template <class T1, class T2> | |
177 | inline typename tools::promote_args<T1, T2>::type | |
178 | spherical_harmonic_r(unsigned n, int m, T1 theta, T2 phi) | |
179 | { | |
180 | return boost::math::spherical_harmonic_r(n, m, theta, phi, policies::policy<>()); | |
181 | } | |
182 | ||
183 | template <class T1, class T2, class Policy> | |
184 | inline typename tools::promote_args<T1, T2>::type | |
185 | spherical_harmonic_i(unsigned n, int m, T1 theta, T2 phi, const Policy& pol) | |
186 | { | |
187 | typedef typename tools::promote_args<T1, T2>::type result_type; | |
188 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
189 | return policies::checked_narrowing_cast<result_type, Policy>(detail::spherical_harmonic_i(n, m, static_cast<value_type>(theta), static_cast<value_type>(phi), pol), "boost::math::spherical_harmonic_i<%1%>(unsigned, int, %1%, %1%)"); | |
190 | } | |
191 | ||
192 | template <class T1, class T2> | |
193 | inline typename tools::promote_args<T1, T2>::type | |
194 | spherical_harmonic_i(unsigned n, int m, T1 theta, T2 phi) | |
195 | { | |
196 | return boost::math::spherical_harmonic_i(n, m, theta, phi, policies::policy<>()); | |
197 | } | |
198 | ||
199 | } // namespace math | |
200 | } // namespace boost | |
201 | ||
202 | #endif // BOOST_MATH_SPECIAL_SPHERICAL_HARMONIC_HPP | |
203 | ||
204 | ||
205 |