ceph/src/boost/libs/math/include/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/math_fwd.hpp>
  15 #include <boost/math/special_functions/legendre.hpp>
  16 #include <boost/math/tools/workaround.hpp>
  17 #include <complex>
  18
  19 namespace boost{
  20 namespace math{
  21
  22 namespace detail{
  23
  24 //
  25 // Calculates the prefix term that's common to the real
  26 // and imaginary parts.  Does *not* fix up the sign of the result
  27 // though.
  28 //
  29 template <class T, class Policy>
  30 inline T spherical_harmonic_prefix(unsigned n, unsigned m, T theta, const Policy& pol)
  31 {
  32    BOOST_MATH_STD_USING
  33
  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 //
  50 template <class T, class Policy>
  51 T spherical_harmonic_r(unsigned n, int m, T theta, T phi, const Policy& pol)
  52 {
  53    BOOST_MATH_STD_USING  // ADL of std functions
  54
  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
  77 template <class T, class Policy>
  78 T spherical_harmonic_i(unsigned n, int m, T theta, T phi, const Policy& pol)
  79 {
  80    BOOST_MATH_STD_USING  // ADL of std functions
  81
  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
 104 template <class T, class U, class Policy>
 105 std::complex<T> spherical_harmonic(unsigned n, int m, U theta, U phi, const Policy& pol)
 106 {
 107    BOOST_MATH_STD_USING
 108    //
 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 detail
 150
 151 template <class T1, class T2, class Policy>
 152 inline 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
 160 template <class T1, class T2>
 161 inline 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
 167 template <class T1, class T2, class Policy>
 168 inline 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), "bost::math::spherical_harmonic_r<%1%>(unsigned, int, %1%, %1%)");
 174 }
 175
 176 template <class T1, class T2>
 177 inline 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
 183 template <class T1, class T2, class Policy>
 184 inline 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
 192 template <class T1, class T2>
 193 inline 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 math
 200 } // namespace boost
 201
 202 #endif // BOOST_MATH_SPECIAL_SPHERICAL_HARMONIC_HPP
 203
 204
 205