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

// kelvin.hpp
// Kelvin functions ber, bei, ker, kei (order 0)
//
// Definitions (DLMF §10.61, Wolfram MathWorld):
//   ber(x) + i bei(x) = J_0(x · e^{3iπ/4})    = J_0(x·(-1+i)/√2)
//   ker(x) + i kei(x) = K_0(x · e^{iπ/4})     = K_0(x·(1+i)/√2)
//
// Provided functions:
//   ber(x), bei(x)     — Kelvin functions of the first kind (from Bessel J), entire real domain
//   ker(x), kei(x)     — Kelvin functions of the second kind (from modified Bessel K), x > 0
//
// Applications: EE skin effect, magnetic field calculation of a magnetic core in a solenoid coil,
//      vibrating-beam problems, etc.
//
// Implementation: delegates directly to the existing sangi::special::besselJ(int, Complex<double>) and
//       sangi::special::besselK(int, Complex<double>).
//
// Properties:
//   ber, bei are even functions of x (because J_0 is even)
//   ker, kei diverge as x → 0+ (ker diverges logarithmically), undefined for x < 0 → NaN

#ifndef SANGI_SPECIAL_KELVIN_HPP
#define SANGI_SPECIAL_KELVIN_HPP

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

namespace sangi {
namespace special {

namespace detail_kelvin {

// 1/√2 as a high-precision constant
constexpr double inv_sqrt2 = 0.7071067811865475244008443621048490392848359376884740365883398690;

}  // namespace detail_kelvin

/// ber(x) — Kelvin function of the first kind (real part)
template<IsNativeFloat T>
[[nodiscard]] T ber(T x) {
    if (std::isnan(x)) return std::numeric_limits<T>::quiet_NaN();
    if (x == T(0))     return T(1);

    double xd = static_cast<double>(std::abs(x));  // even function, so the absolute value is OK
    double s  = xd * detail_kelvin::inv_sqrt2;
    Complex<double> z(-s, s);                      // x·e^(3iπ/4)
    return static_cast<T>(besselJ(0, z).re);
}

/// bei(x) — Kelvin function of the first kind (imaginary part)
template<IsNativeFloat T>
[[nodiscard]] T bei(T x) {
    if (std::isnan(x)) return std::numeric_limits<T>::quiet_NaN();
    if (x == T(0))     return T(0);

    double xd = static_cast<double>(std::abs(x));
    double s  = xd * detail_kelvin::inv_sqrt2;
    Complex<double> z(-s, s);
    return static_cast<T>(besselJ(0, z).im);
}

/// ker(x) — Kelvin function of the second kind (real part), defined for x > 0
template<IsNativeFloat T>
[[nodiscard]] T ker(T x) {
    if (std::isnan(x)) return std::numeric_limits<T>::quiet_NaN();
    if (x <= T(0)) {
        if (x == T(0)) return std::numeric_limits<T>::infinity();
        return std::numeric_limits<T>::quiet_NaN();
    }
    if (std::isinf(x)) return T(0);  // decays to 0 as x → ∞

    double xd = static_cast<double>(x);
    double s  = xd * detail_kelvin::inv_sqrt2;
    Complex<double> z(s, s);                       // x·e^(iπ/4)
    return static_cast<T>(besselK(0, z).re);
}

/// kei(x) — Kelvin function of the second kind (imaginary part), defined for x > 0
template<IsNativeFloat T>
[[nodiscard]] T kei(T x) {
    if (std::isnan(x)) return std::numeric_limits<T>::quiet_NaN();
    if (x <= T(0)) {
        // kei(0) = -π/4 (finite value), x < 0 is undefined
        if (x == T(0)) return -std::numbers::pi_v<T> / T(4);
        return std::numeric_limits<T>::quiet_NaN();
    }
    if (std::isinf(x)) return T(0);

    double xd = static_cast<double>(x);
    double s  = xd * detail_kelvin::inv_sqrt2;
    Complex<double> z(s, s);
    return static_cast<T>(besselK(0, z).im);
}

} // namespace special
} // namespace sangi

#endif // SANGI_SPECIAL_KELVIN_HPP