206 lines · plain
1 2// (C) Copyright John Maddock 2006.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_MATH_SPECIAL_SPHERICAL_HARMONIC_HPP8#define BOOST_MATH_SPECIAL_SPHERICAL_HARMONIC_HPP9 10#ifdef _MSC_VER11#pragma once12#endif13 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 19namespace boost{20namespace math{21 22namespace detail{23 24//25// Calculates the prefix term that's common to the real26// and imaginary parts. Does *not* fix up the sign of the result27// though.28//29template <class T, class Policy>30inline T spherical_harmonic_prefix(unsigned n, unsigned m, T theta, const Policy& pol)31{32 BOOST_MATH_STD_USING33 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//50template <class T, class Policy>51T spherical_harmonic_r(unsigned n, int m, T theta, T phi, const Policy& pol)52{53 BOOST_MATH_STD_USING // ADL of std functions54 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 77template <class T, class Policy>78T spherical_harmonic_i(unsigned n, int m, T theta, T phi, const Policy& pol)79{80 BOOST_MATH_STD_USING // ADL of std functions81 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 104template <class T, class U, class Policy>105std::complex<T> spherical_harmonic(unsigned n, int m, U theta, U phi, const Policy& pol)106{107 BOOST_MATH_STD_USING108 //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 detail150 151template <class T1, class T2, class Policy>152inline 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 160template <class T1, class T2>161inline 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 167template <class T1, class T2, class Policy>168inline 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), "boost::math::spherical_harmonic_r<%1%>(unsigned, int, %1%, %1%)");174}175 176template <class T1, class T2>177inline 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 183template <class T1, class T2, class Policy>184inline 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 192template <class T1, class T2>193inline 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 math200} // namespace boost201 202#endif // BOOST_MATH_SPECIAL_SPHERICAL_HARMONIC_HPP203 204 205 206