boost/math/special_functions/spherical_harmonic.hpp

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