// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// debye.hpp
// Debye functions (solid-state physics)
//
// Provided functions:
//   debye(n, x) — Debye function D_n(x) = (n/x^n) ∫₀ˣ t^n/(e^t - 1) dt
//   debye1(x)   — D_1(x)
//   debye2(x)   — D_2(x)
//   debye3(x)   — D_3(x) (specific heat)
//   debye4(x)   — D_4(x) (enthalpy)
//
// Supported types: float, double, long double

#ifndef SANGI_SPECIAL_DEBYE_HPP
#define SANGI_SPECIAL_DEBYE_HPP

#include <math/special/special_common.hpp>
#include <cmath>
#include <limits>
#include <numbers>

namespace sangi {
namespace special {

// ================================================================
// Debye function D_n(x) = (n/x^n) ∫₀ˣ t^n/(e^t - 1) dt
// ================================================================

namespace detail {

template<IsNativeFloat T>
T debye_series(int n, T x) {
    // For small x: Taylor series
    // D_n(x) = 1 - n·x/(2(n+1)) + Σ_{k=1}^∞ B_{2k}·n!·x^{2k} / ((n+2k)·(2k)!)
    // where B_{2k} are the Bernoulli numbers
    //
    // Simplified version: compute ∫₀ˣ t^n/(e^t-1) dt directly via Gauss-Legendre

    // Integral: ∫₀ˣ t^n/(e^t - 1) dt
    // The integrand is asymptotic to t^{n-1} as t→0 (integrable)

    // 15-point Gauss-Legendre on [0, x]
    // Substitution: t = x·(u+1)/2, dt = x/2 du, u ∈ [-1, 1]
    constexpr T gl_nodes[] = {
        T(-0.98799251802048543), T(-0.93727339240070591), T(-0.84820658341042722),
        T(-0.72441773136017005), T(-0.57097217260853885), T(-0.39415134707756337),
        T(-0.20119409399743452), T(0.0),
        T(0.20119409399743452), T(0.39415134707756337), T(0.57097217260853885),
        T(0.72441773136017005), T(0.84820658341042722), T(0.93727339240070591),
        T(0.98799251802048543)
    };
    constexpr T gl_weights[] = {
        T(0.03075324199611727), T(0.07036604748890871), T(0.10715922046717194),
        T(0.13957067792615431), T(0.16626920581699393), T(0.18616100001556222),
        T(0.19843148532711158), T(0.20257824192556127),
        T(0.19843148532711158), T(0.18616100001556222), T(0.16626920581699393),
        T(0.13957067792615431), T(0.10715922046717194), T(0.07036604748890871),
        T(0.03075324199611727)
    };

    T sum = T(0);
    T half_x = x / T(2);
    for (int i = 0; i < 15; ++i) {
        T t = half_x * (gl_nodes[i] + T(1));
        if (t < T(1e-15)) continue;

        T integrand;
        if (t > T(500)) {
            integrand = std::pow(t, T(n)) * std::exp(-t);
        } else {
            integrand = std::pow(t, T(n)) / (std::exp(t) - T(1));
        }
        sum += gl_weights[i] * integrand;
    }
    sum *= half_x;

    // D_n(x) = n / x^n · integral
    return T(n) / std::pow(x, T(n)) * sum;
}

} // namespace detail

/// Debye function D_n(x)
template<IsNativeFloat T>
[[nodiscard]] T debye(int n, T x) {
    if (std::isnan(x)) return std::numeric_limits<T>::quiet_NaN();
    if (n < 1) return std::numeric_limits<T>::quiet_NaN();
    if (x < T(0)) return std::numeric_limits<T>::quiet_NaN();
    if (x == T(0)) return T(1);  // D_n(0) = 1

    // For very large x: D_n(x) → n!/x^n · ζ(n+1) · Γ(n+1) ...
    // Asymptotic: D_n(x) → n · n! · ζ(n+1) / x^n (x → ∞)

    return detail::debye_series(n, x);
}

template<IsNativeFloat T>
[[nodiscard]] T debye1(T x) { return debye(1, x); }

template<IsNativeFloat T>
[[nodiscard]] T debye2(T x) { return debye(2, x); }

template<IsNativeFloat T>
[[nodiscard]] T debye3(T x) { return debye(3, x); }

template<IsNativeFloat T>
[[nodiscard]] T debye4(T x) { return debye(4, x); }

} // namespace special
} // namespace sangi

#endif // SANGI_SPECIAL_DEBYE_HPP