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