Quadrature.h

#pragma once
 
#include <array>
#include <cmath>
#include <numbers>
 
namespace environment::field::aggregation
{
 
template <int N>
struct GaussLegendre
{
   std::array<double, N> node {};
   std::array<double, N> weight {};
 
   GaussLegendre()
   {
      for (int i = 0; i < N; ++i) {
         double x = std::cos(std::numbers::pi * (i + 0.75) / (N + 0.5));
         double derivative = 1.0;
         for (int iteration = 0; iteration < 100; ++iteration) {
            double p0 = 1.0;
            double p1 = x;
            for (int k = 2; k <= N; ++k) {
               const double p2 = ((2.0 * k - 1.0) * x * p1 - (k - 1.0) * p0) / k;
               p0              = p1;
               p1              = p2;
            }
            derivative   = N * (x * p1 - p0) / (x * x - 1.0);
            const double step = p1 / derivative;
            x -= step;
            if (std::abs(step) < 1e-15) {
               break;
            }
         }
         node[static_cast<std::size_t>(i)]   = x;
         weight[static_cast<std::size_t>(i)] = 2.0 / ((1.0 - x * x) * derivative * derivative);
      }
   }
};
 
template <int N>
const GaussLegendre<N>& Rule()
{
   static const GaussLegendre<N> rule;
   return rule;
}
 
} // namespace environment::field::aggregation