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

/**
 * @file vector_calculus.hpp
 * @brief ベクトル解析 (CAS-37〜39)
 *
 * CAS-37: ベクトルラプラシアン
 * CAS-38: テンソル微分 (発散, 回転)
 * CAS-39: 座標系変換 (直交, 円柱, 球)
 *
 * 数値微分ベースの実装。関数 R³ → R³ を受け取る。
 */

#ifndef CALX_MATH_CORE_VECTOR_CALCULUS_HPP
#define CALX_MATH_CORE_VECTOR_CALCULUS_HPP

#include <array>
#include <functional>
#include <cmath>

namespace calx {

using Vec3 = std::array<double, 3>;
using ScalarField = std::function<double(double, double, double)>;
using VectorField = std::function<Vec3(double, double, double)>;

/// 勾配 ∇f (数値微分)
inline Vec3 gradient(const ScalarField& f, double x, double y, double z,
                      double h = 1e-5) {
    return {{
        (f(x+h,y,z) - f(x-h,y,z)) / (2*h),
        (f(x,y+h,z) - f(x,y-h,z)) / (2*h),
        (f(x,y,z+h) - f(x,y,z-h)) / (2*h)
    }};
}

/// 発散 ∇·F (数値微分)
inline double divergence(const VectorField& F, double x, double y, double z,
                          double h = 1e-5) {
    auto Fxp = F(x+h,y,z), Fxm = F(x-h,y,z);
    auto Fyp = F(x,y+h,z), Fym = F(x,y-h,z);
    auto Fzp = F(x,y,z+h), Fzm = F(x,y,z-h);
    return (Fxp[0]-Fxm[0] + Fyp[1]-Fym[1] + Fzp[2]-Fzm[2]) / (2*h);
}

/// 回転 ∇×F (数値微分)
inline Vec3 curl(const VectorField& F, double x, double y, double z,
                  double h = 1e-5) {
    auto Fyp = F(x,y+h,z), Fym = F(x,y-h,z);
    auto Fzp = F(x,y,z+h), Fzm = F(x,y,z-h);
    auto Fxp = F(x+h,y,z), Fxm = F(x-h,y,z);
    return {{
        (Fyp[2]-Fym[2] - (Fzp[1]-Fzm[1])) / (2*h),
        (Fzp[0]-Fzm[0] - (Fxp[2]-Fxm[2])) / (2*h),
        (Fxp[1]-Fxm[1] - (Fyp[0]-Fym[0])) / (2*h)
    }};
}

/// スカラーラプラシアン ∇²f (数値微分)
inline double laplacian(const ScalarField& f, double x, double y, double z,
                         double h = 1e-5) {
    double fc = f(x,y,z);
    return (f(x+h,y,z) + f(x-h,y,z) + f(x,y+h,z) + f(x,y-h,z)
          + f(x,y,z+h) + f(x,y,z-h) - 6*fc) / (h*h);
}

/// ベクトルラプラシアン ∇²F (CAS-37)
inline Vec3 vectorLaplacian(const VectorField& F, double x, double y, double z,
                             double h = 1e-5) {
    auto Fc = F(x,y,z);
    auto lap = [&](int comp) -> double {
        auto Fxp = F(x+h,y,z), Fxm = F(x-h,y,z);
        auto Fyp = F(x,y+h,z), Fym = F(x,y-h,z);
        auto Fzp = F(x,y,z+h), Fzm = F(x,y,z-h);
        return (Fxp[comp]+Fxm[comp]+Fyp[comp]+Fym[comp]+Fzp[comp]+Fzm[comp]-6*Fc[comp]) / (h*h);
    };
    return {{lap(0), lap(1), lap(2)}};
}

// ========================================================================
// CAS-39: 座標系変換
// ========================================================================

/// 直交座標 → 円柱座標 (r, θ, z)
inline Vec3 cartesianToCylindrical(double x, double y, double z) {
    return {{std::sqrt(x*x + y*y), std::atan2(y, x), z}};
}

/// 円柱座標 → 直交座標
inline Vec3 cylindricalToCartesian(double r, double theta, double z) {
    return {{r * std::cos(theta), r * std::sin(theta), z}};
}

/// 直交座標 → 球座標 (r, θ, φ) — θ=極角, φ=方位角
inline Vec3 cartesianToSpherical(double x, double y, double z) {
    double r = std::sqrt(x*x + y*y + z*z);
    return {{r, std::acos(z / (r + 1e-30)), std::atan2(y, x)}};
}

/// 球座標 → 直交座標
inline Vec3 sphericalToCartesian(double r, double theta, double phi) {
    return {{r * std::sin(theta) * std::cos(phi),
             r * std::sin(theta) * std::sin(phi),
             r * std::cos(theta)}};
}

} // namespace calx

#endif
