#pragma once

// algoim::TanhSinhQuadrature implements basic tanh-sinh quadrature methods.
// Some of the design choices follow the discussion in the paper
//    R. I. Saye, High-order quadrature on multi-component domains implicitly
//    defined by multivariate polynomials, Journal of Computational Physics,
//    448, 110720 (2022), https://doi.org/10.1016/j.jcp.2021.110720

#include <array>
#include <cassert>
#include <type_traits>

template <typename T>
struct TSIntegrator {
    static_assert(std::is_floating_point_v<T>);

    static constexpr int nMax = 100;

    // Quadrature node, relative to the interval [-1,1], for an n-point scheme
    static inline constexpr T x(int n, int i)
    {
        assert(1 <= n && n <= nMax && 0 <= i && i < n);
        return data[n * (n - 1) + 2 * i + 0];
    }

    // Quadrature weight, relative to the interval [-1,1], for an n-point scheme
    static inline constexpr T w(int n, int i)
    {
        assert(1 <= n && n <= nMax && 0 <= i && i < n);
        return data[n * (n - 1) + 2 * i + 1];
    }

    // Quadrature node, relative to the interval [a,b], for an n-point scheme
    static inline constexpr T x(int n, int i, T a, T b) { return (a + b + (b - a) * x(n, i)) * 0.5; }

    // Quadrature weight, relative to the interval [a,b], for an n-point scheme
    static inline constexpr T w(int n, int i, T a, T b) { return w(n, i) * (b - a) * 0.5; }

    static constexpr std::array<T, nMax *(nMax + 1)> data = {
    // clang-format off
        #include "data_cache/TS_weights.hpp"
        // clang-format on
    };
};