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

// anger_weber.hpp
// Anger function J_ν and Weber function E_ν
//
// Definitions (DLMF 11.10.1, 11.10.2):
//   J_ν(z) = (1/π) ∫₀^π cos(νt - z sin t) dt        [Anger]
//   E_ν(z) = (1/π) ∫₀^π sin(νt - z sin t) dt        [Weber]
//
// Functions provided:
//   angerJ(nu, z)  — Anger function J_ν(z)
//   weberE(nu, z)  — Weber function E_ν(z)
//
// Properties:
//   When ν is an integer n: J_n^Anger(z) = J_n^Bessel(z) (matches Bessel)
//   For non-integer ν it is a distinct, independent function family differing from Bessel
//   Both are entire functions for any finite z
//
// Implementation:
//   Evaluate the defining integral directly via Simpson's rule (N=512).
//   The ν series representation is not used, as the integer-ν limit handling is complex.
//   The integrand is analytic and smooth; with N=512, accuracy is about 1e-13 for |z|≲50, |ν|≲30.
//
// Uses:
//   Bessel variants, special integral evaluations, diffraction theory, non-integer-order oscillation problems
//
// References: DLMF §11.10, A&S §12.3, Watson "Bessel Functions" §10.1

#ifndef SANGI_SPECIAL_ANGER_WEBER_HPP
#define SANGI_SPECIAL_ANGER_WEBER_HPP

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

namespace sangi {
namespace special {

namespace detail_aw {

// Evaluate ∫₀^π cos(νt - z sin t) dt or ∫₀^π sin(νt - z sin t) dt via
// the composite Simpson's rule. is_sine=false selects cos, true selects sin.
template<IsNativeFloat T>
[[nodiscard]] T simpson_integrate(T nu, T z, bool is_sine, int N = 1024) {
    // Force N to be even
    if (N % 2) N++;
    constexpr T PI = std::numbers::pi_v<T>;
    T h = PI / T(N);

    auto integrand = [&](T t) -> T {
        T arg = nu * t - z * std::sin(t);
        return is_sine ? std::sin(arg) : std::cos(arg);
    };

    T sum = integrand(T(0)) + integrand(PI);  // endpoints
    // Odd indices (coefficient 4)
    for (int k = 1; k < N; k += 2) {
        sum += T(4) * integrand(T(k) * h);
    }
    // Even indices (coefficient 2)
    for (int k = 2; k < N; k += 2) {
        sum += T(2) * integrand(T(k) * h);
    }
    return sum * h / T(3);
}

}  // namespace detail_aw

// ================================================================
// Anger function J_ν(z) = (1/π) ∫₀^π cos(νt - z sin t) dt
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T angerJ(T nu, T z) {
    if (std::isnan(nu) || std::isnan(z))
        return std::numeric_limits<T>::quiet_NaN();
    constexpr T PI = std::numbers::pi_v<T>;
    return detail_aw::simpson_integrate(nu, z, /*is_sine=*/false) / PI;
}

// ================================================================
// Weber function E_ν(z) = (1/π) ∫₀^π sin(νt - z sin t) dt
// ================================================================

template<IsNativeFloat T>
[[nodiscard]] T weberE(T nu, T z) {
    if (std::isnan(nu) || std::isnan(z))
        return std::numeric_limits<T>::quiet_NaN();
    constexpr T PI = std::numbers::pi_v<T>;
    return detail_aw::simpson_integrate(nu, z, /*is_sine=*/true) / PI;
}

}  // namespace special
}  // namespace sangi

#endif  // SANGI_SPECIAL_ANGER_WEBER_HPP