// Copyright (c) 2024-2026 Tarik Moussa.
// SPDX-License-Identifier: MIT

// test_spherical_functional.cpp  (Phase 3c + 3e)
//
// Phase 3c — SphericalFunctional ported to ConformalMesh.
//
// Corresponds to de.varylab.discreteconformal.functional.SphericalFunctionalTest.
//
// Test map (Java → C++)
// ──────────────────────
//   testHessian (Ignored)         → GradientCheck_Hessian     (ported)
//   testGradientWithHyperIdeal…   → GradientCheck_OctaFaceVertex   (ported)
//   testGradientInExtendedDomain  → GradientCheck_SpherTetVertex    (ported)
//   testGradientWithHyperelliptic → GradientCheck_SpherTetAllDofs   (ported)
//   testFunctionalAtNaNValue      → AnglesFiniteAtKnownPoint         (ported)
//
// Energy model
// ────────────
//   The energy is computed as the Schläfli path integral E(x) = ∫₀¹⟨G(tx),x⟩dt
//   using 10-point Gauss-Legendre quadrature. The gradient check therefore
//   verifies that G is curl-free (the integrability / exactness condition of
//   the spherical discrete conformal functional). This is equivalent to the
//   Java FunctionalTest gradient check.

#include "conformal_mesh.hpp"
#include "mesh_builder.hpp"
#include "spherical_functional.hpp"
#include "spherical_hessian.hpp"
#include "spherical_geometry.hpp"
#include "clausen.hpp"
#include <gtest/gtest.h>
#include <cmath>
#include <vector>

using namespace conformallab;

// ════════════════════════════════════════════════════════════════════════════
// Cross-module Hessian check: spherical_gradient() ↔ spherical_hessian()
//
// Java @Ignore reason: "no Hessian implemented" — the Java functional test
// was written before the Hessian existed.  In C++ the analytic spherical
// Hessian (spherical_hessian.hpp, Phase 3f) is complete.
//
// This test verifies cross-module consistency between the functional and
// the Hessian module.  The spherical Hessian is NSD (negative semi-definite)
// because the spherical energy is concave — hessian_check_spherical() uses
// the sign-corrected FD check appropriate for the spherical case.
// ════════════════════════════════════════════════════════════════════════════

TEST(SphericalFunctional, GradientCheck_Hessian)
{
    auto mesh = make_spherical_tetrahedron();
    auto maps = setup_spherical_maps(mesh);
    compute_lambda0_from_mesh(mesh, maps);
    int n = assign_vertex_dof_indices(mesh, maps);

    std::vector<double> x(static_cast<std::size_t>(n), -0.2);

    EXPECT_TRUE(hessian_check_spherical(mesh, x, maps))
        << "Cross-module: spherical_gradient() and spherical_hessian() are inconsistent";
}

// ════════════════════════════════════════════════════════════════════════════
// Angle formula: octahedron-face triangle has all angles = π/2
//
// The triangle (1,0,0)–(0,1,0)–(0,0,1) has l_ij = π/2 for all edges.
// Half-angle formula: s = 3π/4, s_ij = π/4 for all three.
// All angles = π/2 (right-angled spherical triangle).
// ════════════════════════════════════════════════════════════════════════════

TEST(SphericalFunctional, OctaFaceAnglesAreRightAngles)
{
    // l_ij = π/2 for all edges (octahedron face on unit sphere)
    const double l = PI_SPHER / 2.0;
    auto fa = spherical_angles(l, l, l);

    ASSERT_TRUE(fa.valid) << "Equilateral spherical triangle must be valid";
    EXPECT_NEAR(PI_SPHER / 2.0, fa.alpha1, 1e-12);
    EXPECT_NEAR(PI_SPHER / 2.0, fa.alpha2, 1e-12);
    EXPECT_NEAR(PI_SPHER / 2.0, fa.alpha3, 1e-12);
}

// ════════════════════════════════════════════════════════════════════════════
// Angle sum of a spherical triangle exceeds π (positive curvature)
//
// For the spherical tetrahedron face (arccos(−1/3) ≈ 1.9106 per edge):
// The dihedral angle = arccos(1/3) ≈ 70.53°; by symmetry the face angles
// (vertex angles of the spherical triangle) are all equal.
// Angle sum must be > π and equal 3·arccos(1/3) ≈ 3·1.2310 ≈ 3.693 rad.
// ════════════════════════════════════════════════════════════════════════════

TEST(SphericalFunctional, SpherTetAngleSumExceedsPi)
{
    // Edge length of spherical tetrahedron face: arccos(−1/3)
    const double l = std::acos(-1.0 / 3.0);
    auto fa = spherical_angles(l, l, l);

    ASSERT_TRUE(fa.valid);
    EXPECT_GT(fa.alpha1 + fa.alpha2 + fa.alpha3, PI_SPHER)
        << "Angle sum of spherical triangle must exceed π";

    // By symmetry all three angles must be equal
    EXPECT_NEAR(fa.alpha1, fa.alpha2, 1e-12);
    EXPECT_NEAR(fa.alpha2, fa.alpha3, 1e-12);

    // For a regular spherical tetrahedron with edge arccos(−1/3):
    // half-angle: tan(α/2) = √(sin(l/2)/sin(3l/2)) = √3 → α/2 = π/3 → α = 2π/3.
    // (arccos(1/3) ≈ 1.231 is the 3D dihedral angle of a Euclidean tetrahedron, not this.)
    double expected = 2.0 * PI_SPHER / 3.0;   // 120°
    EXPECT_NEAR(fa.alpha1, expected, 1e-10);
}

// ════════════════════════════════════════════════════════════════════════════
// Gradient check: octahedron-face triangle, vertex DOFs only
//
// Sets λ° from mesh geometry (unit sphere), all u_i = −0.3 (slightly smaller).
// Mirrors Java testGradientWithHyperIdeal… on a single-triangle mesh.
// ════════════════════════════════════════════════════════════════════════════

TEST(SphericalFunctional, GradientCheck_OctaFaceVertex)
{
    auto mesh = make_octahedron_face();
    auto maps = setup_spherical_maps(mesh);
    compute_lambda0_from_mesh(mesh, maps);
    int n = assign_vertex_dof_indices(mesh, maps);

    // Small uniform conformal factor: shrink the triangle slightly.
    std::vector<double> x(static_cast<std::size_t>(n), -0.3);

    EXPECT_TRUE(gradient_check_spherical(mesh, x, maps))
        << "Gradient check failed on octahedron-face triangle (vertex DOFs)";
}

// ════════════════════════════════════════════════════════════════════════════
// Gradient check: spherical tetrahedron (4 faces), vertex DOFs only
//
// Closed surface; exercises accumulation over multiple faces per vertex.
// Mirrors Java testGradientInTheExtendedDomain.
// ════════════════════════════════════════════════════════════════════════════

TEST(SphericalFunctional, GradientCheck_SpherTetVertex)
{
    auto mesh = make_spherical_tetrahedron();
    auto maps = setup_spherical_maps(mesh);
    compute_lambda0_from_mesh(mesh, maps);
    int n = assign_vertex_dof_indices(mesh, maps);

    std::vector<double> x(static_cast<std::size_t>(n), -0.2);

    EXPECT_TRUE(gradient_check_spherical(mesh, x, maps))
        << "Gradient check failed on spherical tetrahedron (vertex DOFs)";
}

// ════════════════════════════════════════════════════════════════════════════
// Gradient check: spherical tetrahedron, all DOFs (vertex + edge)
//
// Exercises the edge-gradient branch: G_e = α_opp⁺ + α_opp⁻ − π.
// Mirrors Java testGradientWithHyperellipticCurve.
// ════════════════════════════════════════════════════════════════════════════

TEST(SphericalFunctional, GradientCheck_SpherTetAllDofs)
{
    auto mesh = make_spherical_tetrahedron();
    auto maps = setup_spherical_maps(mesh);
    compute_lambda0_from_mesh(mesh, maps);
    int n = assign_all_spherical_dof_indices(mesh, maps);

    // Replacement parameterization (Finding 3): when an edge carries a DOF its
    // value *replaces* λ°_ij + u_i + u_j entirely, so Λ_ij = λ_e.  Here the edge
    // DOFs stay at 0 and only the vertex DOFs are perturbed; this checks that the
    // gradient is curl-free (energy = Schläfli path integral), not Java-faithfulness
    // of the edge formula — that is locked separately by
    // EdgeGradient_RegularTetClosedForm below.
    std::vector<double> x(static_cast<std::size_t>(n), 0.0);
    // Set vertex DOFs (indices 0..3) to -0.2 to keep triangle well-formed.
    for (int i = 0; i < 4; ++i) x[static_cast<std::size_t>(i)] = -0.2;

    EXPECT_TRUE(gradient_check_spherical(mesh, x, maps))
        << "Gradient check failed on spherical tetrahedron (all DOFs)";
}

// ════════════════════════════════════════════════════════════════════════════
// Closed-form oracle for the edge-DOF gradient (Finding 3, missing-test item 4)
//
// The FD gradient check above can only confirm that G is conservative — the
// spherical energy is *defined* as the path integral of G, so the energy↔gradient
// FD agreement is automatic and CANNOT detect a wrong-but-conservative edge
// formula.  This test instead pins the edge gradient against an independent,
// closed-form geometric value, so it would fail if the Finding-3 formula
// (G_e = α_opp⁺ + α_opp⁻ − θ_e, dropping the additive −(S⁺+S⁻)/2 term) ever
// regressed.
//
// Geometry: the regular spherical tetrahedron has all edges a = arccos(−1/3),
// so by the spherical law of cosines every interior corner angle is
//   cos α = (cos a − cos²a)/sin²a = cos a/(1+cos a) = (−1/3)/(2/3) = −1/2
//   ⇒ α = 2π/3.
// Each edge is shared by two faces, so both opposite angles equal 2π/3 and
//   G_e = 2π/3 + 2π/3 − θ_e   with θ_e = π (default)   = π/3.
//
// Setup: all edges carry DOFs, set to their λ⁰ (the replacement convention then
// reproduces the original tetrahedron metric exactly), vertex DOFs left at 0.
// ════════════════════════════════════════════════════════════════════════════

TEST(SphericalFunctional, EdgeGradient_RegularTetClosedForm)
{
    const double PI_ = std::acos(-1.0);

    auto mesh = make_spherical_tetrahedron();
    auto maps = setup_spherical_maps(mesh);
    compute_lambda0_from_mesh(mesh, maps);
    int n = assign_all_spherical_dof_indices(mesh, maps);

    // Edge DOF = λ⁰  → Λ_ij = λ⁰  → reproduces the arccos(−1/3) tetrahedron.
    // Vertex DOFs stay at 0 (ignored by the replacement convention for DOF edges).
    std::vector<double> x(static_cast<std::size_t>(n), 0.0);
    int n_edge_dofs = 0;
    for (auto e : mesh.edges()) {
        int ie = maps.e_idx[e];
        if (ie >= 0) { x[static_cast<std::size_t>(ie)] = maps.lambda0[e]; ++n_edge_dofs; }
    }
    ASSERT_EQ(n_edge_dofs, 6) << "regular tetrahedron must have 6 edge DOFs";

    auto G = spherical_gradient(mesh, x, maps);

    const double expected = PI_ / 3.0;   // 2·(2π/3) − π
    for (auto e : mesh.edges()) {
        int ie = maps.e_idx[e];
        if (ie < 0) continue;
        EXPECT_NEAR(G[static_cast<std::size_t>(ie)], expected, 1e-9)
            << "edge gradient at DOF " << ie
            << " must equal the closed-form value π/3 (Finding 3)";
    }
}

// ════════════════════════════════════════════════════════════════════════════
// Angles are finite at a known interior point
//
// Mirrors Java testFunctionalAtNaNValue: choose DOFs that could hit
// a degenerate branch (l_ij → 0 or triangle inequality fails) and check
// the gradient vector is free of NaN/Inf.
// ════════════════════════════════════════════════════════════════════════════

TEST(SphericalFunctional, AnglesFiniteAtKnownPoint)
{
    auto mesh = make_spherical_tetrahedron();
    auto maps = setup_spherical_maps(mesh);
    compute_lambda0_from_mesh(mesh, maps);
    int n = assign_vertex_dof_indices(mesh, maps);

    // u_i = -1.5: contracts the triangle heavily but stays non-degenerate.
    std::vector<double> x(static_cast<std::size_t>(n), -1.5);
    auto G = spherical_gradient(mesh, x, maps);

    for (std::size_t i = 0; i < G.size(); ++i) {
        EXPECT_FALSE(std::isnan(G[i])) << "Gradient component " << i << " is NaN";
        EXPECT_FALSE(std::isinf(G[i])) << "Gradient component " << i << " is Inf";
    }
}

// ════════════════════════════════════════════════════════════════════════════
// Gradient check: fan-4 mesh on unit sphere, vertex DOFs only
//
// Make a fan of 4 triangles around the north pole (0,0,1);
// rim vertices projected onto the equator.
// Exercises high-valence vertex gradient accumulation.
// ════════════════════════════════════════════════════════════════════════════

TEST(SphericalFunctional, GradientCheck_SpherFan4Vertex)
{
    // Build a fan with 4 spherical triangles manually (can't use make_fan
    // directly because those vertices are not on the unit sphere).
    ConformalMesh mesh;
    auto center = mesh.add_vertex(Point3(0, 0, 1));  // north pole
    const int n_rim = 4;
    std::vector<Vertex_index> rim(n_rim);
    const double dtheta = 2.0 * PI_SPHER / n_rim;
    const double phi    = PI_SPHER / 4.0;  // 45° colatitude
    for (int i = 0; i < n_rim; ++i) {
        double theta = i * dtheta;
        rim[i] = mesh.add_vertex(Point3(
            std::sin(phi) * std::cos(theta),
            std::sin(phi) * std::sin(theta),
            std::cos(phi)));
    }
    for (int i = 0; i < n_rim; ++i)
        mesh.add_face(center, rim[i], rim[(i + 1) % n_rim]);

    auto maps = setup_spherical_maps(mesh);
    compute_lambda0_from_mesh(mesh, maps);
    int ndof = assign_vertex_dof_indices(mesh, maps);

    std::vector<double> x(static_cast<std::size_t>(ndof), -0.3);

    EXPECT_TRUE(gradient_check_spherical(mesh, x, maps))
        << "Gradient check failed on spherical fan-4 mesh";
}

// ════════════════════════════════════════════════════════════════════════════
// Gradient check: mixed pinned/variable vertices
//
// One vertex pinned (u_v = 0 fixed), others variable.
// Verifies that the gradient accumulation skips pinned vertices correctly.
// ════════════════════════════════════════════════════════════════════════════

TEST(SphericalFunctional, GradientCheck_MixedPinnedVertices)
{
    auto mesh = make_octahedron_face();
    auto maps = setup_spherical_maps(mesh);
    compute_lambda0_from_mesh(mesh, maps);

    // Pin v0, make v1 and v2 variable.
    auto vit = mesh.vertices().begin();
    Vertex_index v0 = *vit++;
    Vertex_index v1 = *vit++;
    Vertex_index v2 = *vit;

    maps.v_idx[v0] = -1;   // pinned
    maps.v_idx[v1] = 0;
    maps.v_idx[v2] = 1;

    std::vector<double> x = {-0.2, -0.4};

    EXPECT_TRUE(gradient_check_spherical(mesh, x, maps))
        << "Gradient check failed for mixed pinned/variable vertices";
}

// ════════════════════════════════════════════════════════════════════════════
// Phase 3e — Gauge-fix for closed spherical surfaces
//
// On a closed spherical surface, the functional has a gauge mode:
//   E(u + t·1)  is maximised at some t*.
//   At t*, the sum of all vertex gradients equals zero: Σ G_v = 0.
//
// Test: start from a point with non-zero ΣG_v, apply the gauge shift,
// and verify ΣG_v(x + t*·1) ≈ 0.
// ════════════════════════════════════════════════════════════════════════════

TEST(SphericalFunctional, GaugeFix_SpherTetVertexZerosSumGv)
{
    auto mesh = make_spherical_tetrahedron();
    auto maps = setup_spherical_maps(mesh);
    compute_lambda0_from_mesh(mesh, maps);
    int n = assign_vertex_dof_indices(mesh, maps);

    // Off-gauge starting point: all u_i = -0.5
    std::vector<double> x(static_cast<std::size_t>(n), -0.5);

    // Compute ΣG_v before gauge shift.
    {
        auto G = spherical_gradient(mesh, x, maps);
        double sum = 0.0;
        for (auto v : mesh.vertices()) {
            int iv = maps.v_idx[v];
            if (iv >= 0) sum += G[static_cast<std::size_t>(iv)];
        }
        // At -0.5 the surface is compressed; ΣG_v should be non-zero.
        EXPECT_NE(sum, 0.0) << "Pre-gauge ΣG_v should be non-zero";
    }

    // Compute gauge shift and apply.
    double t = spherical_gauge_shift(mesh, x, maps);
    std::vector<double> x_fixed = x;
    for (auto v : mesh.vertices()) {
        int iv = maps.v_idx[v];
        if (iv >= 0)
            x_fixed[static_cast<std::size_t>(iv)] += t;
    }

    // Verify ΣG_v ≈ 0 at the gauge-fixed point.
    {
        auto G = spherical_gradient(mesh, x_fixed, maps);
        double sum = 0.0;
        for (auto v : mesh.vertices()) {
            int iv = maps.v_idx[v];
            if (iv >= 0) sum += G[static_cast<std::size_t>(iv)];
        }
        EXPECT_NEAR(sum, 0.0, 1e-6)
            << "After gauge fix, Σ G_v should vanish; t* = " << t;
    }
}

TEST(SphericalFunctional, GaugeFix_ApplyInPlace)
{
    auto mesh = make_spherical_tetrahedron();
    auto maps = setup_spherical_maps(mesh);
    compute_lambda0_from_mesh(mesh, maps);
    int n = assign_vertex_dof_indices(mesh, maps);

    // x = -0.3: compressed but inside the valid spherical domain.
    std::vector<double> x(static_cast<std::size_t>(n), -0.3);

    apply_spherical_gauge(mesh, x, maps);

    // After in-place gauge fix, ΣG_v must be near 0.
    auto G = spherical_gradient(mesh, x, maps);
    double sum = 0.0;
    for (auto v : mesh.vertices()) {
        int iv = maps.v_idx[v];
        if (iv >= 0) sum += G[static_cast<std::size_t>(iv)];
    }
    EXPECT_NEAR(sum, 0.0, 1e-6)
        << "apply_spherical_gauge must drive Σ G_v to zero";
}

TEST(SphericalFunctional, GaugeFix_AlreadyAtGaugeReturnsTNearZero)
{
    // A symmetric, equilateral spherical tetrahedron at x=0 is already
    // at the gauge maximum (by symmetry, ΣG_v = 0).
    auto mesh = make_spherical_tetrahedron();
    auto maps = setup_spherical_maps(mesh);
    compute_lambda0_from_mesh(mesh, maps);
    int n = assign_vertex_dof_indices(mesh, maps);

    std::vector<double> x(static_cast<std::size_t>(n), 0.0);
    double t = spherical_gauge_shift(mesh, x, maps);
    // Symmetric starting point → t* should be very close to 0.
    EXPECT_NEAR(t, 0.0, 1e-5)
        << "Gauge shift from the symmetric point should be ~0; got " << t;
}

// ─────────────────────────────────────────────────────────────────────────────
// Golden-value oracle — pin the spherical law-of-cosines angle formula, the β
// half-angle relations, and the Lobachevsky energy term bit-for-bit against the
// upstream Java reference (SphericalFunctional.triangleEnergyAndAlphas, lines
// 401-453), captured by running the compiled Java library (openjdk 17) with the
// real de.varylab…Clausen.Л on these exact arc lengths.  Unlike the Schläfli
// path-integral gradient check (which only verifies curl-freeness), this locks
// the absolute angle/β/energy values against an independent implementation,
// catching any silent index/sign/convention drift in the spherical port.
//
// To regenerate: /tmp/oracle/SphereOracle.java (recipe in doc/reviewer/
// java-port-audit.md).  Values are Java printf %.17g.  Tolerance 1e-12.
//
// Index map (C++ spherical_angles ↔ Java): with l12=lij, l23=ljk, l31=lki the
// returned alpha1/alpha2/alpha3 are Java αi/αj/αk (angle opposite ljk/lki/lij).
// ─────────────────────────────────────────────────────────────────────────────
TEST(SphericalGoldenJava, AngleBetaEnergyFromLengths)
{
    auto check = [](double lij, double ljk, double lki,
                    double ai_g, double aj_g, double ak_g,
                    double bi_g, double bj_g, double bk_g, double L_g) {
        auto fa = spherical_angles(lij, ljk, lki);
        EXPECT_TRUE(fa.valid);
        EXPECT_NEAR(fa.alpha1, ai_g, 1e-12);
        EXPECT_NEAR(fa.alpha2, aj_g, 1e-12);
        EXPECT_NEAR(fa.alpha3, ak_g, 1e-12);

        const double ai = fa.alpha1, aj = fa.alpha2, ak = fa.alpha3;
        const double bi = 0.5 * (PI + ai - aj - ak);
        const double bj = 0.5 * (PI - ai + aj - ak);
        const double bk = 0.5 * (PI - ai - aj + ak);
        EXPECT_NEAR(bi, bi_g, 1e-12);
        EXPECT_NEAR(bj, bj_g, 1e-12);
        EXPECT_NEAR(bk, bk_g, 1e-12);

        const double Lterm =
            Lobachevsky(ai) + Lobachevsky(aj) + Lobachevsky(ak) +
            Lobachevsky(bi) + Lobachevsky(bj) + Lobachevsky(bk) +
            Lobachevsky(0.5 * (PI - ai - aj - ak));
        EXPECT_NEAR(Lterm, L_g, 1e-12);
    };

    // Scalene spherical triangle (all Δ > 0, Δ_ijk < 2π).
    check(1.0, 1.2, 0.9,
          1.5305813141072122, 0.99684981272794600, 1.1246069519101605,
          1.2753586015294494, 0.74162710015018330, 0.86938423933239760,
          1.3576100185550408);
    // Equilateral spherical triangle.
    check(0.7, 0.7, 0.7,
          1.1225596283199812, 1.1225596283199812, 1.1225596283199812,
          1.0095165126349062, 1.0095165126349060, 1.0095165126349060,
          1.6806828584976297);
}

// ─────────────────────────────────────────────────────────────────────────────
// FULL-MESH golden oracle (spherical) — drives the REAL upstream SphericalFunctional
// on the shared tetrahedron (scaled so every arc length stays < π) and pins BOTH
// the per-vertex gradient G_v = Θ_v − Σα AND ΔE = E(x) − E(0) bit-for-bit.
//
// IMPORTANT — the oracle calls Java's `conformalEnergyAndGradient` (the RAW
// θ−Σα gradient), NOT `evaluate()`: `evaluate()` first runs a 1-D Brent
// maximization over the global-scale gauge direction (`maximizeInNegativeDirection`),
// which the C++ `spherical_gradient` deliberately does NOT — C++ factors that
// gauge into the Newton solver's `spherical_gauge_shift` instead.  Comparing
// against `evaluate()` would (wrongly) drive every component to ~0.  The raw
// gradient is the piece that corresponds 1:1 to the C++ functional.
//
// Setup parity (UnwrapUtility.prepareInvariantDataHyperbolicAndSpherical, scale):
//   closed mesh, ALL 4 vertices variable, Θ_v = 2π, no edge DOFs,
//   λ°_e = 2·log(SCALE·|p_i − p_j|)  (Java uses chord length × scale, whereas the
//   C++ compute_lambda0_from_mesh helper assumes unit-sphere vertices and uses
//   ARC length — so we set λ° directly here to match Java exactly),
//   per-vertex u(P) = 0.10·X − 0.07·Y + 0.13·Z, SCALE = 0.2.
//
// To regenerate: /tmp/oracle/{tet.obj,SphereMeshOracle.java}.  Values are Java %.17g.
// ─────────────────────────────────────────────────────────────────────────────
TEST(SphericalGoldenJava, FullMeshGradientAndEnergy_Tetrahedron)
{
    constexpr double TWO_PI = 2.0 * 3.14159265358979323846264338328;
    constexpr double SCALE  = 0.2;

    auto mesh = make_tetrahedron();   // same 4 vertices as /tmp/oracle/tet.obj
    auto maps = setup_spherical_maps(mesh);

    // λ°_e = 2·log(SCALE · chord), matching Java's prepareInvariantData(scale).
    for (auto e : mesh.edges()) {
        auto h  = mesh.halfedge(e);
        const auto& p1 = mesh.point(mesh.source(h));
        const auto& p2 = mesh.point(mesh.target(h));
        const double dx = p1.x() - p2.x(), dy = p1.y() - p2.y(), dz = p1.z() - p2.z();
        const double chord = std::sqrt(dx*dx + dy*dy + dz*dz);
        maps.lambda0[e] = 2.0 * std::log(SCALE * chord);
    }

    int idx = 0;
    for (auto v : mesh.vertices()) {
        maps.v_idx[v]   = idx++;
        maps.theta_v[v] = TWO_PI;
    }

    auto u_of = [](const Point3& p) {
        return 0.10 * p.x() - 0.07 * p.y() + 0.13 * p.z();
    };

    std::vector<double> x(static_cast<std::size_t>(idx), 0.0);
    for (auto v : mesh.vertices())
        x[static_cast<std::size_t>(maps.v_idx[v])] = u_of(mesh.point(v));

    auto G = spherical_gradient(mesh, x, maps);

    struct GoldRow { double X, Y, Z, G; };
    const GoldRow gold[4] = {
        { 1,  1,  1, 2.7671034786104927},
        { 1, -1, -1, 2.5216834054857546},
        {-1,  1, -1, 1.5401761310866633},
        {-1, -1,  1, 2.6518569531861100},
    };
    for (auto v : mesh.vertices()) {
        const auto& p = mesh.point(v);
        const double g = G[static_cast<std::size_t>(maps.v_idx[v])];
        bool matched = false;
        for (const auto& row : gold) {
            if (std::abs(p.x() - row.X) < 1e-9 &&
                std::abs(p.y() - row.Y) < 1e-9 &&
                std::abs(p.z() - row.Z) < 1e-9) {
                EXPECT_NEAR(g, row.G, 1e-12)
                    << "gradient mismatch at (" << p.x() << "," << p.y()
                    << "," << p.z() << ")";
                matched = true;
                break;
            }
        }
        EXPECT_TRUE(matched) << "unexpected vertex position";
    }

    std::vector<double> x0(static_cast<std::size_t>(idx), 0.0);
    const double dE = spherical_energy(mesh, x, maps)
                    - spherical_energy(mesh, x0, maps);
    EXPECT_NEAR(dE, 0.16409141487397116, 1e-12);
}

// ─────────────────────────────────────────────────────────────────────────────
// FULL-MESH EDGE-DOF golden oracle (spherical, Finding 3) — the missing solution-
// level Java oracle for the edge-DOF gradient path.  Makes two opposite edges of
// the shared tetrahedron variable (replacement parameterization: λ_e = x[e_idx],
// no u-shift) and pins BOTH the vertex gradient (Θ−Σα) AND the edge gradient
// (α_opp⁺ + α_opp⁻ − θ_e, θ_e = π) bit-for-bit against the REAL upstream
// SphericalFunctional.conformalEnergyAndGradient (raw gradient, not evaluate()).
//
// This is a pure GRADIENT oracle (no ΔE): with an edge DOF, x = 0 means λ_e = 0
// ⇒ spherical arc length = π (degenerate), so the path-integral-from-origin
// energy reference is ill-defined — the gradient at a fixed non-degenerate x is
// the unambiguous Finding-3 quantity.  It does NOT touch the spherical Hessian,
// so it is unaffected by the Finding-4 edge-DOF Hessian guard.
//
// To regenerate: /tmp/oracle/SphereEdgeOracle.java.  Values are Java %.17g.
// ─────────────────────────────────────────────────────────────────────────────
TEST(SphericalGoldenJava, FullMeshEdgeDofGradient_Tetrahedron)
{
    constexpr double TWO_PI = 2.0 * 3.14159265358979323846264338328;
    constexpr double PI_C   = 3.14159265358979323846264338328;
    constexpr double SCALE  = 0.2;

    auto mesh = make_tetrahedron();
    auto maps = setup_spherical_maps(mesh);

    for (auto e : mesh.edges()) {
        auto h  = mesh.halfedge(e);
        const auto& p1 = mesh.point(mesh.source(h));
        const auto& p2 = mesh.point(mesh.target(h));
        const double dx = p1.x()-p2.x(), dy = p1.y()-p2.y(), dz = p1.z()-p2.z();
        maps.lambda0[e] = 2.0 * std::log(SCALE * std::sqrt(dx*dx + dy*dy + dz*dz));
        maps.theta_e[e] = PI_C;
        maps.e_idx[e]   = -1;
    }

    int idx = 0;
    for (auto v : mesh.vertices()) { maps.v_idx[v] = idx++; maps.theta_v[v] = TWO_PI; }

    auto pos = [&](Vertex_index v){ return mesh.point(v); };
    auto same = [](const Point3& p, double X, double Y, double Z) {
        return std::abs(p.x()-X)<1e-9 && std::abs(p.y()-Y)<1e-9 && std::abs(p.z()-Z)<1e-9;
    };
    // Edge connects positions (aX,aY,aZ)–(bX,bY,bZ) in either order?
    auto connects = [&](Edge_index e, const double a[3], const double b[3]) {
        auto h = mesh.halfedge(e);
        const Point3& s = pos(mesh.source(h));
        const Point3& t = pos(mesh.target(h));
        return (same(s,a[0],a[1],a[2]) && same(t,b[0],b[1],b[2])) ||
               (same(s,b[0],b[1],b[2]) && same(t,a[0],a[1],a[2]));
    };
    const double A[3]={1,1,1}, B[3]={1,-1,-1}, C[3]={-1,1,-1}, D[3]={-1,-1,1};

    // Make edges A–B and C–D variable (replacement parameterization).
    int edof = idx;  // edge DOFs start after the 4 vertex DOFs
    Edge_index eAB, eCD;
    for (auto e : mesh.edges()) {
        if (connects(e, A, B)) { maps.e_idx[e] = edof++; eAB = e; }
        else if (connects(e, C, D)) { maps.e_idx[e] = edof++; eCD = e; }
    }
    ASSERT_EQ(edof, 6);  // 4 vertex + 2 edge DOFs

    auto u_of = [](const Point3& p){ return 0.10*p.x() - 0.07*p.y() + 0.13*p.z(); };

    std::vector<double> x(6, 0.0);
    for (auto v : mesh.vertices())
        x[static_cast<std::size_t>(maps.v_idx[v])] = u_of(mesh.point(v));
    // Edge DOF value = λ⁰_e + 0.1 (same perturbation as the Java oracle).
    x[static_cast<std::size_t>(maps.e_idx[eAB])] = maps.lambda0[eAB] + 0.1;
    x[static_cast<std::size_t>(maps.e_idx[eCD])] = maps.lambda0[eCD] + 0.1;

    auto G = spherical_gradient(mesh, x, maps);

    // Vertex gradients keyed by position.
    struct VRow { double X, Y, Z, G; };
    const VRow vg[4] = {
        { 1,  1,  1, 2.5036751008175546},
        { 1, -1, -1, 2.2303362601050205},
        {-1,  1, -1, 1.8463710314817632},
        {-1, -1,  1, 2.7499719487341570},
    };
    for (auto v : mesh.vertices()) {
        const auto& p = mesh.point(v);
        const double g = G[static_cast<std::size_t>(maps.v_idx[v])];
        bool m = false;
        for (auto& r : vg)
            if (same(p, r.X, r.Y, r.Z)) { EXPECT_NEAR(g, r.G, 1e-12); m = true; break; }
        EXPECT_TRUE(m);
    }
    // Edge gradients: A–B and C–D (Java golden values).
    EXPECT_NEAR(G[static_cast<std::size_t>(maps.e_idx[eAB])], -0.35189517043413690, 1e-12);
    EXPECT_NEAR(G[static_cast<std::size_t>(maps.e_idx[eCD])], -0.44101986058895950, 1e-12);
}

// ════════════════════════════════════════════════════════════════════════════
// Degenerate spherical triangle — limiting angles (Finding-F, java-port-audit item 1)
//
// spherical_angles() must return the limiting angles (π opposite the
// over-long edge, 0/0 elsewhere) when the spherical triangle inequality
// is violated, with valid = false.  This mirrors the Euclidean behaviour
// and is required for the convex C¹ BPS extension.
// ════════════════════════════════════════════════════════════════════════════

TEST(SphericalFunctional, DegenerateTriangle_LimitingAngles_S12TooLong)
{
    // s12 > s23 + s31: s12 = 2.5, s23 = s31 = 0.5  (all < π so valid arc lengths)
    // s23 < 0 → actually use s-based check
    // Easier: use s12 = π − ε (nearly degenerate hemisphere edge)
    // and very short s23, s31 so s12 > s23 + s31.
    const double s12 = 2.0, s23 = 0.4, s31 = 0.4;   // s12 > s23+s31 = 0.8
    auto fa = spherical_angles(s12, s23, s31);
    EXPECT_FALSE(fa.valid);
    // s12 is the edge opposite v3 → α3 = π
    EXPECT_NEAR(fa.alpha3, PI, 1e-12)
        << "corner opposite over-long s12 must be π";
    EXPECT_NEAR(fa.alpha1, 0.0, 1e-12);
    EXPECT_NEAR(fa.alpha2, 0.0, 1e-12);
}

TEST(SphericalFunctional, DegenerateTriangle_LimitingAngles_S23TooLong)
{
    // s23 > s12 + s31 → α1 = π
    const double s12 = 0.4, s23 = 2.0, s31 = 0.4;
    auto fa = spherical_angles(s12, s23, s31);
    EXPECT_FALSE(fa.valid);
    EXPECT_NEAR(fa.alpha1, PI, 1e-12)
        << "corner opposite over-long s23 must be π";
    EXPECT_NEAR(fa.alpha2, 0.0, 1e-12);
    EXPECT_NEAR(fa.alpha3, 0.0, 1e-12);
}

TEST(SphericalFunctional, DegenerateTriangle_LimitingAngles_S31TooLong)
{
    // s31 > s12 + s23 → α2 = π
    const double s12 = 0.4, s23 = 0.4, s31 = 2.0;
    auto fa = spherical_angles(s12, s23, s31);
    EXPECT_FALSE(fa.valid);
    EXPECT_NEAR(fa.alpha2, PI, 1e-12)
        << "corner opposite over-long s31 must be π";
    EXPECT_NEAR(fa.alpha1, 0.0, 1e-12);
    EXPECT_NEAR(fa.alpha3, 0.0, 1e-12);
}