math/special_functions/spherical_harmonic.hpp

0001
0002 //  (C) Copyright John Maddock 2006.
0003 //  Use, modification and distribution are subject to the
0004 //  Boost Software License, Version 1.0. (See accompanying file
0005 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_MATH_SPECIAL_SPHERICAL_HARMONIC_HPP
0008 #define BOOST_MATH_SPECIAL_SPHERICAL_HARMONIC_HPP
0009
0010 #ifdef _MSC_VER
0011 #pragma once
0012 #endif
0013
0014 #include <boost/math/special_functions/math_fwd.hpp>
0015 #include <boost/math/special_functions/legendre.hpp>
0016 #include <boost/math/tools/workaround.hpp>
0017 #include <complex>
0018
0019 namespace boost{
0020 namespace math{
0021
0022 namespace detail{
0023
0024 //
0025 // Calculates the prefix term that's common to the real
0026 // and imaginary parts.  Does *not* fix up the sign of the result
0027 // though.
0028 //
0029 template <class T, class Policy>
0030 inline T spherical_harmonic_prefix(unsigned n, unsigned m, T theta, const Policy& pol)
0031 {
0032    BOOST_MATH_STD_USING
0033
0034    if(m > n)
0035       return 0;
0036
0037    T sin_theta = sin(theta);
0038    T x = cos(theta);
0039
0040    T leg = detail::legendre_p_imp(n, m, x, static_cast<T>(pow(fabs(sin_theta), T(m))), pol);
0041
0042    T prefix = boost::math::tgamma_delta_ratio(static_cast<T>(n - m + 1), static_cast<T>(2 * m), pol);
0043    prefix *= (2 * n + 1) / (4 * constants::pi<T>());
0044    prefix = sqrt(prefix);
0045    return prefix * leg;
0046 }
0047 //
0048 // Real Part:
0049 //
0050 template <class T, class Policy>
0051 T spherical_harmonic_r(unsigned n, int m, T theta, T phi, const Policy& pol)
0052 {
0053    BOOST_MATH_STD_USING  // ADL of std functions
0054
0055    bool sign = false;
0056    if(m < 0)
0057    {
0058       // Reflect and adjust sign if m < 0:
0059       sign = m&1;
0060       m = abs(m);
0061    }
0062    if(m&1)
0063    {
0064       // Check phase if theta is outside [0, PI]:
0065       T mod = boost::math::tools::fmod_workaround(theta, T(2 * constants::pi<T>()));
0066       if(mod < 0)
0067          mod += 2 * constants::pi<T>();
0068       if(mod > constants::pi<T>())
0069          sign = !sign;
0070    }
0071    // Get the value and adjust sign as required:
0072    T prefix = spherical_harmonic_prefix(n, m, theta, pol);
0073    prefix *= cos(m * phi);
0074    return sign ? T(-prefix) : prefix;
0075 }
0076
0077 template <class T, class Policy>
0078 T spherical_harmonic_i(unsigned n, int m, T theta, T phi, const Policy& pol)
0079 {
0080    BOOST_MATH_STD_USING  // ADL of std functions
0081
0082    bool sign = false;
0083    if(m < 0)
0084    {
0085       // Reflect and adjust sign if m < 0:
0086       sign = !(m&1);
0087       m = abs(m);
0088    }
0089    if(m&1)
0090    {
0091       // Check phase if theta is outside [0, PI]:
0092       T mod = boost::math::tools::fmod_workaround(theta, T(2 * constants::pi<T>()));
0093       if(mod < 0)
0094          mod += 2 * constants::pi<T>();
0095       if(mod > constants::pi<T>())
0096          sign = !sign;
0097    }
0098    // Get the value and adjust sign as required:
0099    T prefix = spherical_harmonic_prefix(n, m, theta, pol);
0100    prefix *= sin(m * phi);
0101    return sign ? T(-prefix) : prefix;
0102 }
0103
0104 template <class T, class U, class Policy>
0105 std::complex<T> spherical_harmonic(unsigned n, int m, U theta, U phi, const Policy& pol)
0106 {
0107    BOOST_MATH_STD_USING
0108    //
0109    // Sort out the signs:
0110    //
0111    bool r_sign = false;
0112    bool i_sign = false;
0113    if(m < 0)
0114    {
0115       // Reflect and adjust sign if m < 0:
0116       r_sign = m&1;
0117       i_sign = !(m&1);
0118       m = abs(m);
0119    }
0120    if(m&1)
0121    {
0122       // Check phase if theta is outside [0, PI]:
0123       U mod = boost::math::tools::fmod_workaround(theta, U(2 * constants::pi<U>()));
0124       if(mod < 0)
0125          mod += 2 * constants::pi<U>();
0126       if(mod > constants::pi<U>())
0127       {
0128          r_sign = !r_sign;
0129          i_sign = !i_sign;
0130       }
0131    }
0132    //
0133    // Calculate the value:
0134    //
0135    U prefix = spherical_harmonic_prefix(n, m, theta, pol);
0136    U r = prefix * cos(m * phi);
0137    U i = prefix * sin(m * phi);
0138    //
0139    // Add in the signs:
0140    //
0141    if(r_sign)
0142       r = -r;
0143    if(i_sign)
0144       i = -i;
0145    static const char* function = "boost::math::spherical_harmonic<%1%>(int, int, %1%, %1%)";
0146    return std::complex<T>(policies::checked_narrowing_cast<T, Policy>(r, function), policies::checked_narrowing_cast<T, Policy>(i, function));
0147 }
0148
0149 } // namespace detail
0150
0151 template <class T1, class T2, class Policy>
0152 inline std::complex<typename tools::promote_args<T1, T2>::type>
0153    spherical_harmonic(unsigned n, int m, T1 theta, T2 phi, const Policy& pol)
0154 {
0155    typedef typename tools::promote_args<T1, T2>::type result_type;
0156    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0157    return detail::spherical_harmonic<result_type, value_type>(n, m, static_cast<value_type>(theta), static_cast<value_type>(phi), pol);
0158 }
0159
0160 template <class T1, class T2>
0161 inline std::complex<typename tools::promote_args<T1, T2>::type>
0162    spherical_harmonic(unsigned n, int m, T1 theta, T2 phi)
0163 {
0164    return boost::math::spherical_harmonic(n, m, theta, phi, policies::policy<>());
0165 }
0166
0167 template <class T1, class T2, class Policy>
0168 inline typename tools::promote_args<T1, T2>::type
0169    spherical_harmonic_r(unsigned n, int m, T1 theta, T2 phi, const Policy& pol)
0170 {
0171    typedef typename tools::promote_args<T1, T2>::type result_type;
0172    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0173    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%)");
0174 }
0175
0176 template <class T1, class T2>
0177 inline typename tools::promote_args<T1, T2>::type
0178    spherical_harmonic_r(unsigned n, int m, T1 theta, T2 phi)
0179 {
0180    return boost::math::spherical_harmonic_r(n, m, theta, phi, policies::policy<>());
0181 }
0182
0183 template <class T1, class T2, class Policy>
0184 inline typename tools::promote_args<T1, T2>::type
0185    spherical_harmonic_i(unsigned n, int m, T1 theta, T2 phi, const Policy& pol)
0186 {
0187    typedef typename tools::promote_args<T1, T2>::type result_type;
0188    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0189    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%)");
0190 }
0191
0192 template <class T1, class T2>
0193 inline typename tools::promote_args<T1, T2>::type
0194    spherical_harmonic_i(unsigned n, int m, T1 theta, T2 phi)
0195 {
0196    return boost::math::spherical_harmonic_i(n, m, theta, phi, policies::policy<>());
0197 }
0198
0199 } // namespace math
0200 } // namespace boost
0201
0202 #endif // BOOST_MATH_SPECIAL_SPHERICAL_HARMONIC_HPP
0203
0204
0205