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