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

// transport.hpp
// Transport functions and synchrotron functions
//
// Provided functions:
//   transport(n, x)  — transport function J_n(x) = ∫₀ˣ t^n e^t / (e^t - 1)² dt
//   synchrotronF(x)  — synchrotron function F(x) = x ∫ₓ^∞ K_{5/3}(t) dt
//   synchrotronG(x)  — synchrotron function G(x) = x · K_{2/3}(x)
//
// Supported types: float, double, long double

#ifndef SANGI_SPECIAL_TRANSPORT_HPP
#define SANGI_SPECIAL_TRANSPORT_HPP

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

namespace sangi {
namespace special {

// ================================================================
// Transport function J_n(x) = ∫₀ˣ t^n · e^t / (e^t - 1)² dt
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T transport(int n, T x) {
    if (std::isnan(x) || n < 2)
        return std::numeric_limits<T>::quiet_NaN();
    if (x < T(0)) return std::numeric_limits<T>::quiet_NaN();
    if (x == T(0)) return T(0);

    // 15-point Gauss-Legendre quadrature on [0, x]
    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 half_x = x / T(2);
    T sum = T(0);
    for (int i = 0; i < 15; ++i) {
        T t = half_x * (gl_nodes[i] + T(1));
        if (t < T(1e-15)) {
            // t → 0: t^n · e^t / (e^t-1)² ≈ t^{n-2}
            sum += gl_weights[i] * std::pow(t, T(n - 2));
            continue;
        }
        T et = std::exp(t);
        T denom = (et - T(1));
        T integrand = std::pow(t, T(n)) * et / (denom * denom);
        sum += gl_weights[i] * integrand;
    }
    return sum * half_x;
}

// ================================================================
// Synchrotron functions
// F(x) = x ∫ₓ^∞ K_{5/3}(t) dt  (first kind)
// G(x) = x · K_{2/3}(x)         (second kind)
//
// Approximation (Fouka & Ouichaoui 2013):
// F(x) ≈ 2.1495 x^{1/3} (1 + 1.5x^{2/3})^{-1} exp(-x) · H(x)
// G(x) ≈ 1.8084 x^{1/3} (1 + 0.5x^{2/3})^{-1/2} exp(-x)  (simplified)
// ================================================================

namespace detail {

// Approximation of Bessel K_nu(x) (for nu = 2/3, 5/3)
template<IsNativeFloat T>
T besselK_approx(T nu, T x) {
    if (x <= T(0)) return std::numeric_limits<T>::infinity();

    // Large x: K_nu(x) ≈ sqrt(π/(2x)) · exp(-x) · [1 + (4ν²-1)/(8x) + ...]
    if (x > T(2)) {
        T mu = T(4) * nu * nu;
        T sqrtpi2x = std::sqrt(std::numbers::pi_v<T> / (T(2) * x));
        T inv8x = T(1) / (T(8) * x);
        T series = T(1) + (mu - T(1)) * inv8x
                   + (mu - T(1)) * (mu - T(9)) * inv8x * inv8x / T(2)
                   + (mu - T(1)) * (mu - T(9)) * (mu - T(25)) * inv8x * inv8x * inv8x / T(6);
        return sqrtpi2x * std::exp(-x) * series;
    }

    // Small x: numerical integration K_nu(x) = ∫₀^∞ exp(-x·cosh(t))·cosh(νt) dt
    T sum = T(0);
    T dt = T(0.05);
    for (T t = T(0); t < T(20); t += dt) {
        T w = (t == T(0) || t + dt >= T(20)) ? T(0.5) : T(1);
        sum += w * std::exp(-x * std::cosh(t)) * std::cosh(nu * t);
    }
    return sum * dt;
}

} // namespace detail

/// Synchrotron function F(x) = x ∫ₓ^∞ K_{5/3}(t) dt
template<IsNativeFloat T>
[[nodiscard]] T synchrotronF(T x) {
    if (std::isnan(x) || x < T(0))
        return std::numeric_limits<T>::quiet_NaN();
    if (x == T(0)) return T(0);

    // Numerical integration: ∫ₓ^∞ K_{5/3}(t) dt (trapezoidal rule)
    T sum = T(0);
    T dt = x * T(0.02);
    if (dt < T(0.01)) dt = T(0.01);
    if (dt > T(0.5)) dt = T(0.5);

    T prev = detail::besselK_approx(T(5.0 / 3.0), x);
    for (T t = x + dt; t < x + T(100); t += dt) {
        T curr = detail::besselK_approx(T(5.0 / 3.0), t);
        sum += (prev + curr) * T(0.5) * dt;
        prev = curr;
        if (curr < std::numeric_limits<T>::epsilon() * sum)
            break;
    }

    return x * sum;
}

/// Synchrotron function G(x) = x · K_{2/3}(x)
template<IsNativeFloat T>
[[nodiscard]] T synchrotronG(T x) {
    if (std::isnan(x) || x < T(0))
        return std::numeric_limits<T>::quiet_NaN();
    if (x == T(0)) return T(0);

    return x * detail::besselK_approx(T(2.0 / 3.0), x);
}

} // namespace special
} // namespace sangi

#endif // SANGI_SPECIAL_TRANSPORT_HPP