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

// spherical_geometry.hpp
//
// Pure-math building blocks for the spherical discrete conformal map.
// Ported from de.varylab.discreteconformal.functional.SphericalFunctional
// (the geometry helpers embedded there).
//
// Notation:
//   u_i          – vertex conformal factor (DOF)
//   λ°_e         – base log-length of edge e (fixed initial value)
//   λ_ij         – effective log-length = λ°_ij + u_i + u_j
//   l_ij         – spherical arc length = 2·asin(min(exp(λ_ij/2), 1))
//   α_k          – interior angle of the spherical triangle at vertex k

#include "constants.hpp"
#include <cmath>
#include <algorithm>

namespace conformallab {

/// Backward-compatible alias — prefer conformallab::PI in new code.
constexpr double PI_SPHER = PI;

// ── Effective spherical arc length ────────────────────────────────────────────

/// Spherical arc length `l(λ) = 2·asin(min(exp(λ/2), 1))`.
/// Clamps `exp(λ/2)` to `[0, 1]` so `asin` stays in domain.
inline double spherical_l(double lambda)
{
    double half = std::exp(lambda * 0.5);
    if (half >= 1.0) half = ASIN_DOMAIN_GUARD;
    if (half <= 0.0) return 0.0;
    return 2.0 * std::asin(half);
}

// ── Interior angles of a spherical triangle ──────────────────────────────────

/// Interior angles of a spherical triangle, plus a `valid` flag.
struct SphericalFaceAngles {
    double alpha1; ///< Corner angle at vertex v₁.
    double alpha2; ///< Corner angle at vertex v₂.
    double alpha3; ///< Corner angle at vertex v₃.
    bool   valid;  ///< `false` when the three lengths violate the spherical triangle inequality.
};

/// Compute the spherical-triangle corner angles `(α₁, α₂, α₃)` from
/// the three arc lengths `(l₁₂, l₂₃, l₃₁)` using the half-angle form
/// of the spherical law of cosines.  Returns `valid = false` for
/// degenerate or out-of-range triangles.
inline SphericalFaceAngles spherical_angles(double l12, double l23, double l31)
{
    double s   = (l12 + l23 + l31) * 0.5;
    double s12 = s - l12;
    double s23 = s - l23;
    double s31 = s - l31;

    // Degenerate spherical triangle: return the *limiting* angles, matching the
    // Java reference (SphericalFunctional.triangleEnergyAndAlphas).  a1 is the
    // angle opposite l23, a2 opposite l31, a3 opposite l12.  `valid` stays false
    // so the Hessian still skips the face, but the gradient uses these angles
    // (convex C¹ extension onto the infeasible region).
    //   s12<=0 (Δij<=0) → corner opposite l12 = π  → a3 = π
    //   s23<=0 (Δjk<=0) → corner opposite l23 = π  → a1 = π
    //   s31<=0 (Δki<=0) → corner opposite l31 = π  → a2 = π
    //   s>=π  (Δijk>=2π) → all three corners = π
    if (s12 <= 0.0) return {0.0, 0.0, PI_SPHER, false};
    if (s23 <= 0.0) return {PI_SPHER, 0.0, 0.0, false};
    if (s31 <= 0.0) return {0.0, PI_SPHER, 0.0, false};
    if (s >= PI_SPHER) return {PI_SPHER, PI_SPHER, PI_SPHER, false};

    const double ss   = std::sin(s);
    const double ss12 = std::sin(s12);
    const double ss23 = std::sin(s23);
    const double ss31 = std::sin(s31);

    double a1 = 2.0 * std::atan2(std::sqrt(ss12 * ss31), std::sqrt(ss * ss23));
    double a2 = 2.0 * std::atan2(std::sqrt(ss12 * ss23), std::sqrt(ss * ss31));
    double a3 = 2.0 * std::atan2(std::sqrt(ss23 * ss31), std::sqrt(ss * ss12));

    return {a1, a2, a3, true};
}

} // namespace conformallab