spherical_harmonic.hpp 6.2 KB

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