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

// coulomb.hpp
// Coulomb wave functions
//
// Provided functions:
//   coulombF(L, eta, rho) — regular Coulomb wave function F_L(η, ρ)
//   coulombG(L, eta, rho) — irregular Coulomb wave function G_L(η, ρ)
//   coulombCL(L, eta)     — Coulomb normalization constant C_L(η)
//
// Supported types: float, double, long double

#ifndef SANGI_SPECIAL_COULOMB_HPP
#define SANGI_SPECIAL_COULOMB_HPP

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

namespace sangi {
namespace special {

// ================================================================
// Coulomb normalization constant
// C_L(η) = 2^L exp(-πη/2) |Γ(L+1+iη)| / Γ(2L+2)
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T coulombCL(int L, T eta) {
    // |Γ(L+1+iη)|² = π·η / (sinh(πη)) · Π_{s=0}^{L} (s² + η²)  (for L >= 0)
    // For L=0: C_0(η) = sqrt(2πη / (exp(2πη) - 1))
    T pi = std::numbers::pi_v<T>;
    T two_pi_eta = T(2) * pi * eta;

    T p = T(1);
    if (std::abs(eta) < T(1e-15)) {
        // η → 0: C_L → simplification of the form 2^L / (2L+1)!!
        p = T(1);
        for (int s = 1; s <= L; ++s) p *= T(s);
        T denom = T(1);
        for (int s = 1; s <= 2 * L + 1; ++s) denom *= T(s);
        return std::pow(T(2), T(L)) * p / denom;
    }

    // Sommerfeld parameter
    T eta_factor;
    if (two_pi_eta > T(500)) {
        eta_factor = two_pi_eta * std::exp(-two_pi_eta);
    } else if (two_pi_eta < T(-500)) {
        eta_factor = -two_pi_eta;
    } else {
        eta_factor = two_pi_eta / (std::exp(two_pi_eta) - T(1));
    }

    T C0_sq = eta_factor;
    T CL_sq = C0_sq;
    for (int s = 1; s <= L; ++s) {
        CL_sq *= (T(s * s) + eta * eta) / (T(s) * T(2 * s + 1));
    }

    return std::sqrt(CL_sq);
}

// ================================================================
// Regular Coulomb wave function F_L(η, ρ)
// Taylor series: F_L = C_L · ρ^{L+1} Σ A_k ρ^k
// ================================================================

namespace detail {

template<IsNativeFloat T>
void coulomb_fg_series(int L, T eta, T rho,
                       T& F_val, T& G_val)
{
    // F_L via power series: F_L = C_L(η) · ρ^{L+1} · Σ a_k ρ^k
    // a_0 = 1, a_1 = η/(L+1)
    // a_{k+1} = (2ηa_k - a_{k-1}) / ((k+1)(k+2L+2))
    T CL = coulombCL(L, eta);

    const int max_terms = 400;
    T a_prev = T(0);
    T a_curr = T(1);
    T sum = T(1);

    for (int k = 0; k < max_terms; ++k) {
        T a_next = (T(2) * eta * a_curr - a_prev) /
                   (T(k + 1) * T(k + 2 * L + 2));
        a_prev = a_curr;
        a_curr = a_next;

        T term = a_curr * std::pow(rho, T(k + 1));
        sum += term;

        if (k >= 5 && std::abs(term) <
            std::numeric_limits<T>::epsilon() * std::abs(sum))
            break;
    }

    F_val = CL * std::pow(rho, T(L + 1)) * sum;

    // G_L via numerical-derivative approximation (simplified version)
    T dr = rho * T(1e-6);
    if (dr < T(1e-10)) dr = T(1e-10);
    T sum_plus = T(1);
    a_prev = T(0);
    a_curr = T(1);
    T rho_p = rho + dr;
    for (int k = 0; k < max_terms; ++k) {
        T a_next = (T(2) * eta * a_curr - a_prev) /
                   (T(k + 1) * T(k + 2 * L + 2));
        a_prev = a_curr;
        a_curr = a_next;
        T term = a_curr * std::pow(rho_p, T(k + 1));
        sum_plus += term;
        if (k >= 5 && std::abs(term) <
            std::numeric_limits<T>::epsilon() * std::abs(sum_plus))
            break;
    }
    T F_plus = CL * std::pow(rho_p, T(L + 1)) * sum_plus;
    T F_prime = (F_plus - F_val) / dr;

    // Wronskian: F_L · G_L' - F_L' · G_L = 1
    // G_L ≈ (F_L' · ρ - (L+1) · F_L) / ... (Coulomb phase)
    // Simple approximation: G_L from asymptotic for moderate ρ
    if (rho > T(L + 1) + T(10)) {
        T sigma = T(0);
        for (int k = 1; k <= L; ++k) sigma += std::atan(eta / T(k));
        T theta = rho - eta * std::log(T(2) * rho) -
                  T(L) * std::numbers::pi_v<T> / T(2) + sigma;
        G_val = std::cos(theta) / std::cos(theta - std::numbers::pi_v<T> / T(4));
        // Wronskian correction
        if (std::abs(F_val) > std::numeric_limits<T>::epsilon() * T(10))
            G_val = (T(1) + F_prime * G_val) / F_val;
        else
            G_val = T(0);
    } else {
        // Small ρ: G is singular → NaN
        G_val = std::numeric_limits<T>::quiet_NaN();
    }
}

} // namespace detail

template<IsNativeFloat T>
[[nodiscard]] T coulombF(int L, T eta, T rho) {
    if (std::isnan(eta) || std::isnan(rho) || rho < T(0))
        return std::numeric_limits<T>::quiet_NaN();
    if (rho == T(0)) return T(0);

    T F, G;
    detail::coulomb_fg_series(L, eta, rho, F, G);
    return F;
}

template<IsNativeFloat T>
[[nodiscard]] T coulombG(int L, T eta, T rho) {
    if (std::isnan(eta) || std::isnan(rho) || rho <= T(0))
        return std::numeric_limits<T>::quiet_NaN();

    T F, G;
    detail::coulomb_fg_series(L, eta, rho, F, G);
    return G;
}

} // namespace special
} // namespace sangi

#endif // SANGI_SPECIAL_COULOMB_HPP