Line data    Source code

       1             : /*
       2             :  * Copyright (c) 2016 Sebastian Weber, Henri Menke. All rights reserved.
       3             :  *
       4             :  * This file is part of the pairinteraction library.
       5             :  *
       6             :  * The pairinteraction library is free software: you can redistribute it and/or modify
       7             :  * it under the terms of the GNU Lesser General Public License as published by
       8             :  * the Free Software Foundation, either version 3 of the License, or
       9             :  * (at your option) any later version.
      10             :  *
      11             :  * The pairinteraction library is distributed in the hope that it will be useful,
      12             :  * but WITHOUT ANY WARRANTY; without even the implied warranty of
      13             :  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      14             :  * GNU Lesser General Public License for more details.
      15             :  *
      16             :  * You should have received a copy of the GNU Lesser General Public License
      17             :  * along with the pairinteraction library. If not, see <http://www.gnu.org/licenses/>.
      18             :  */
      19             : 
      20             : #include "WignerD.hpp"
      21             : 
      22          53 : WignerD::WignerD() = default;
      23             : 
      24        4540 : double WignerD::operator()(float j, float m, float mp, double beta) {
      25        4540 :     double tolerance = 1e-16;
      26        4540 :     if (std::abs(beta - M_PI / 2) < tolerance) {
      27         156 :         return this->evalWignerdPiHalf(j, m, mp);
      28             :     }
      29        4384 :     return this->evalWignerd(j, m, mp, beta);
      30             : }
      31             : 
      32        4536 : std::complex<double> WignerD::operator()(float j, float m, float mp, double alpha, double beta,
      33             :                                          double gamma) {
      34           0 :     return std::complex<double>(std::cos(-m * alpha), std::sin(-m * alpha)) *
      35        4536 :         this->operator()(j, m, mp, beta) *
      36       13608 :         std::complex<double>(std::cos(-mp * gamma), std::sin(-mp * gamma));
      37             : }
      38             : 
      39       52628 : double WignerD::evalWignerdPiHalf(float j, float m, float mp) {
      40       52628 :     double r = 0;
      41      138294 :     for (unsigned int k = std::max(0, static_cast<int>(mp - m)); k <= j + std::min(mp, -m); ++k) {
      42       85666 :         r += std::pow(-1, k) * boost::math::binomial_coefficient<double>(j + mp, k) *
      43       85666 :             boost::math::binomial_coefficient<double>(j - mp, k + m - mp);
      44             :     }
      45       52628 :     r *= std::pow(-1., m - mp) / std::pow(2., j) *
      46       52628 :         std::sqrt(
      47       52628 :              boost::math::factorial<double>(j + m) * boost::math::factorial<double>(j - m) /
      48       52628 :              (boost::math::factorial<double>(j + mp) * boost::math::factorial<double>(j - mp)));
      49       52628 :     return r;
      50             : }
      51             : 
      52        4384 : double WignerD::evalWignerd(float j, float m, float mp, double beta) {
      53        4384 :     std::complex<double> r = 0;
      54       30620 :     for (float mpp = j; mpp >= -j; --mpp) { // TODO parallelize
      55       26236 :         r += this->evalWignerdPiHalf(j, m, mpp) *
      56       26236 :             std::complex<double>(std::cos(-mpp * beta), std::sin(-mpp * beta)) *
      57       52472 :             this->evalWignerdPiHalf(j, mpp, -mp);
      58             :     }
      59        4384 :     r *= std::pow(std::complex<double>(0, 1), 2. * j - m - mp) * std::pow(-1., 2. * m);
      60        8768 :     return r.real();
      61             : }