spherical_harmonic.hpp 6.2 KB

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