Chris@16: 
Chris@16: //  (C) Copyright John Maddock 2006.
Chris@16: //  Use, modification and distribution are subject to the
Chris@16: //  Boost Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_MATH_SPECIAL_SPHERICAL_HARMONIC_HPP
Chris@16: #define BOOST_MATH_SPECIAL_SPHERICAL_HARMONIC_HPP
Chris@16: 
Chris@16: #ifdef _MSC_VER
Chris@16: #pragma once
Chris@16: #endif
Chris@16: 
Chris@101: #include <boost/math/special_functions/math_fwd.hpp>
Chris@16: #include <boost/math/special_functions/legendre.hpp>
Chris@16: #include <boost/math/tools/workaround.hpp>
Chris@16: #include <complex>
Chris@16: 
Chris@16: namespace boost{
Chris@16: namespace math{
Chris@16: 
Chris@16: namespace detail{
Chris@16: 
Chris@16: //
Chris@16: // Calculates the prefix term that's common to the real
Chris@16: // and imaginary parts.  Does *not* fix up the sign of the result
Chris@16: // though.
Chris@16: //
Chris@16: template <class T, class Policy>
Chris@16: inline T spherical_harmonic_prefix(unsigned n, unsigned m, T theta, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    if(m > n)
Chris@16:       return 0;
Chris@16: 
Chris@16:    T sin_theta = sin(theta);
Chris@16:    T x = cos(theta);
Chris@16: 
Chris@16:    T leg = detail::legendre_p_imp(n, m, x, static_cast<T>(pow(fabs(sin_theta), T(m))), pol);
Chris@16:    
Chris@16:    T prefix = boost::math::tgamma_delta_ratio(static_cast<T>(n - m + 1), static_cast<T>(2 * m), pol);
Chris@16:    prefix *= (2 * n + 1) / (4 * constants::pi<T>());
Chris@16:    prefix = sqrt(prefix);
Chris@16:    return prefix * leg;
Chris@16: }
Chris@16: //
Chris@16: // Real Part:
Chris@16: //
Chris@16: template <class T, class Policy>
Chris@16: T spherical_harmonic_r(unsigned n, int m, T theta, T phi, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING  // ADL of std functions
Chris@16: 
Chris@16:    bool sign = false;
Chris@16:    if(m < 0)
Chris@16:    {
Chris@16:       // Reflect and adjust sign if m < 0:
Chris@16:       sign = m&1;
Chris@16:       m = abs(m);
Chris@16:    }
Chris@16:    if(m&1)
Chris@16:    {
Chris@16:       // Check phase if theta is outside [0, PI]:
Chris@16:       T mod = boost::math::tools::fmod_workaround(theta, T(2 * constants::pi<T>()));
Chris@16:       if(mod < 0)
Chris@16:          mod += 2 * constants::pi<T>();
Chris@16:       if(mod > constants::pi<T>())
Chris@16:          sign = !sign;
Chris@16:    }
Chris@16:    // Get the value and adjust sign as required:
Chris@16:    T prefix = spherical_harmonic_prefix(n, m, theta, pol);
Chris@16:    prefix *= cos(m * phi);
Chris@16:    return sign ? T(-prefix) : prefix;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: T spherical_harmonic_i(unsigned n, int m, T theta, T phi, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING  // ADL of std functions
Chris@16: 
Chris@16:    bool sign = false;
Chris@16:    if(m < 0)
Chris@16:    {
Chris@16:       // Reflect and adjust sign if m < 0:
Chris@16:       sign = !(m&1);
Chris@16:       m = abs(m);
Chris@16:    }
Chris@16:    if(m&1)
Chris@16:    {
Chris@16:       // Check phase if theta is outside [0, PI]:
Chris@16:       T mod = boost::math::tools::fmod_workaround(theta, T(2 * constants::pi<T>()));
Chris@16:       if(mod < 0)
Chris@16:          mod += 2 * constants::pi<T>();
Chris@16:       if(mod > constants::pi<T>())
Chris@16:          sign = !sign;
Chris@16:    }
Chris@16:    // Get the value and adjust sign as required:
Chris@16:    T prefix = spherical_harmonic_prefix(n, m, theta, pol);
Chris@16:    prefix *= sin(m * phi);
Chris@16:    return sign ? T(-prefix) : prefix;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class U, class Policy>
Chris@16: std::complex<T> spherical_harmonic(unsigned n, int m, U theta, U phi, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    //
Chris@16:    // Sort out the signs:
Chris@16:    //
Chris@16:    bool r_sign = false;
Chris@16:    bool i_sign = false;
Chris@16:    if(m < 0)
Chris@16:    {
Chris@16:       // Reflect and adjust sign if m < 0:
Chris@16:       r_sign = m&1;
Chris@16:       i_sign = !(m&1);
Chris@16:       m = abs(m);
Chris@16:    }
Chris@16:    if(m&1)
Chris@16:    {
Chris@16:       // Check phase if theta is outside [0, PI]:
Chris@16:       U mod = boost::math::tools::fmod_workaround(theta, U(2 * constants::pi<U>()));
Chris@16:       if(mod < 0)
Chris@16:          mod += 2 * constants::pi<U>();
Chris@16:       if(mod > constants::pi<U>())
Chris@16:       {
Chris@16:          r_sign = !r_sign;
Chris@16:          i_sign = !i_sign;
Chris@16:       }
Chris@16:    }
Chris@16:    //
Chris@16:    // Calculate the value:
Chris@16:    //
Chris@16:    U prefix = spherical_harmonic_prefix(n, m, theta, pol);
Chris@16:    U r = prefix * cos(m * phi);
Chris@16:    U i = prefix * sin(m * phi);
Chris@16:    //
Chris@16:    // Add in the signs:
Chris@16:    //
Chris@16:    if(r_sign)
Chris@16:       r = -r;
Chris@16:    if(i_sign)
Chris@16:       i = -i;
Chris@16:    static const char* function = "boost::math::spherical_harmonic<%1%>(int, int, %1%, %1%)";
Chris@16:    return std::complex<T>(policies::checked_narrowing_cast<T, Policy>(r, function), policies::checked_narrowing_cast<T, Policy>(i, function));
Chris@16: }
Chris@16: 
Chris@16: } // namespace detail
Chris@16: 
Chris@16: template <class T1, class T2, class Policy>
Chris@16: inline std::complex<typename tools::promote_args<T1, T2>::type> 
Chris@16:    spherical_harmonic(unsigned n, int m, T1 theta, T2 phi, const Policy& pol)
Chris@16: {
Chris@16:    typedef typename tools::promote_args<T1, T2>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    return detail::spherical_harmonic<result_type, value_type>(n, m, static_cast<value_type>(theta), static_cast<value_type>(phi), pol);
Chris@16: }
Chris@16: 
Chris@16: template <class T1, class T2>
Chris@16: inline std::complex<typename tools::promote_args<T1, T2>::type> 
Chris@16:    spherical_harmonic(unsigned n, int m, T1 theta, T2 phi)
Chris@16: {
Chris@16:    return boost::math::spherical_harmonic(n, m, theta, phi, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T1, class T2, class Policy>
Chris@16: inline typename tools::promote_args<T1, T2>::type 
Chris@16:    spherical_harmonic_r(unsigned n, int m, T1 theta, T2 phi, const Policy& pol)
Chris@16: {
Chris@16:    typedef typename tools::promote_args<T1, T2>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    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%)");
Chris@16: }
Chris@16: 
Chris@16: template <class T1, class T2>
Chris@16: inline typename tools::promote_args<T1, T2>::type 
Chris@16:    spherical_harmonic_r(unsigned n, int m, T1 theta, T2 phi)
Chris@16: {
Chris@16:    return boost::math::spherical_harmonic_r(n, m, theta, phi, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T1, class T2, class Policy>
Chris@16: inline typename tools::promote_args<T1, T2>::type 
Chris@16:    spherical_harmonic_i(unsigned n, int m, T1 theta, T2 phi, const Policy& pol)
Chris@16: {
Chris@16:    typedef typename tools::promote_args<T1, T2>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    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%)");
Chris@16: }
Chris@16: 
Chris@16: template <class T1, class T2>
Chris@16: inline typename tools::promote_args<T1, T2>::type 
Chris@16:    spherical_harmonic_i(unsigned n, int m, T1 theta, T2 phi)
Chris@16: {
Chris@16:    return boost::math::spherical_harmonic_i(n, m, theta, phi, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: } // namespace math
Chris@16: } // namespace boost
Chris@16: 
Chris@16: #endif // BOOST_MATH_SPECIAL_SPHERICAL_HARMONIC_HPP
Chris@16: 
Chris@16: 
Chris@16: