legendre.hpp

Go to the documentation of this file.

// Copyright (c) 2013-2017 Commissariat à l'énergie atomique et aux énergies alternatives (CEA)
// Copyright (c) 2013-2017 Centre national de la recherche scientifique (CNRS)
// Copyright (c) 2019 Simons Foundation
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You may obtain a copy of the License at
//     https://www.gnu.org/licenses/gpl-3.0.txt
//
// Authors: Olivier Parcollet, Nils Wentzell

 
#pragma once
 
#include <array>
#include <complex>
 
#include <nda/macros.hpp>
 
namespace triqs::utility {


  std::complex<double> legendre_T(int n, int l);

  double legendre_t(int l, int p);

  double mod_cyl_bessel_i(int n, double x);

  class legendre_generator {
    public:
    legendre_generator(double x = 0.0) : x_(x), arr_{1.0, x} { EXPECTS(x >= -1.0 and x <= 1.0); }

    double next() {
      if (l_ > 1) {
        auto idx  = static_cast<unsigned int>(l_ % 2);
        arr_[idx] = ((2 * l_ - 1) * x_ * arr_[1 - idx] - (l_ - 1) * arr_[idx]) / l_;
        ++l_;
        return arr_[idx];
      } else {
        ++l_;
        return arr_[l_ - 1];
      }
    }

    void reset(double x) {
      EXPECTS(x >= -1.0 and x <= 1.0);
      x_      = x;
      l_      = 0;
      arr_[0] = 1.0;
      arr_[1] = x;
    }
 
    private:
    double x_{0.0};
    unsigned int l_{0};
    std::array<double, 2> arr_{1.0, 0.0};
  };

 
} // namespace triqs::utility