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

// hyper_ideal_geometry.hpp
//
// Pure-math building blocks for the hyper-ideal discrete conformal map.
// Ported from de.varylab.discreteconformal.functional.HyperIdealUtility
// and the private helpers of HyperIdealFunctional (lij, αij, σi, σij).
//
// All functions are independent of the mesh type.
//
// Notation follows the original Java / paper:
//   b_i, b_j  – vertex variables (log scale factors, hyper-ideal vertices)
//   a_ij      – edge variable (intersection angle between horocycles)
//   l_ij      – effective hyperbolic edge length in the auxiliary triangle
//   β_i       – interior angle of the hyperbolic triangle at vertex i
//   α_ij      – dihedral angle of the tetrahedron at edge ij

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

namespace conformallab {

// ── Length functions ─────────────────────────────────────────────────────────

/// `ζ(x,y,z)` — interior angle (in radians) in a hyperbolic triangle
/// with edge lengths `x`, `y`, `z`, opposite to the side of length `z`.
/// Ports `HyperIdealUtility.ζ(x, y, z)`.
inline double zeta(double x, double y, double z)
{
    double cx = std::cosh(x), cy = std::cosh(y), cz = std::cosh(z);
    double sx = std::sinh(x), sy = std::sinh(y);
    double nbd = (cx*cy - cz) / (sx*sy);
    nbd = std::clamp(nbd, -1.0, 1.0);   // guard floating-point rounding
    return std::acos(nbd);
}

/// `ζ₁₃(x,y,z)` — third edge length in a right-angled hyperbolic hexagon.
/// Ports `HyperIdealUtility.ζ_13(x, y, z)`.
inline double zeta13(double x, double y, double z)
{
    double cx = std::cosh(x), cy = std::cosh(y), cz = std::cosh(z);
    double sx = std::sinh(x), sy = std::sinh(y);
    return std::acosh((cx*cy + cz) / (sx*sy));
}

/// `ζ₁₄(x,y)` — edge length in a hyperbolic pentagon with one ideal vertex.
/// Ports `HyperIdealUtility.ζ_14(x, y)`.
inline double zeta14(double x, double y)
{
    double cy = std::cosh(y), sy = std::sinh(y);
    return std::acosh((std::exp(x) + cy) / sy);
}

/// `ζ₁₅(x)` — length in a hyperbolic quadrilateral with two ideal vertices.
/// Ports `HyperIdealUtility.ζ_15(x)`.
inline double zeta15(double x)
{
    return 2.0 * std::asinh(std::exp(x / 2.0));
}

// ── Effective edge length ─────────────────────────────────────────────────────

/// `l_ij`: effective hyperbolic length of edge ij.
///   * `bi`, `bj` — vertex log scale factors (used only when vertex is hyper-ideal).
///   * `aij`      — edge intersection-angle variable.
///   * `vi_var` / `vj_var` — `true` iff the corresponding vertex is hyper-ideal.
/// Ports `HyperIdealFunctional.lij()`.
inline double lij(double bi, double bj, double aij, bool vi_var, bool vj_var)
{
    if (vi_var && vj_var) return zeta13(bi, bj, aij);
    if (vi_var)           return zeta14(aij, bi);
    if (vj_var)           return zeta14(aij, bj);
    return zeta15(aij);
}

// ── Auxiliary angle functions ─────────────────────────────────────────────────

/// `σᵢ(aᵢⱼ, aₖᵢ, aⱼₖ, vj_var, vk_var)` — intermediate half-length at vertex i.
/// Ports `HyperIdealFunctional.σi()`.
inline double sigma_i(double aij, double aki, double ajk, bool vj_var, bool vk_var)
{
    if (vj_var && vk_var) return zeta13(aij, aki, ajk);
    if (vj_var)           return zeta14(ajk - aki, aij);
    if (vk_var)           return zeta14(ajk - aij, aki);
    return zeta15(ajk - aij - aki);
}

/// `σᵢⱼ(aᵢⱼ, bᵢ, bⱼ, vj_var)` — intermediate half-length for edge ij from vertex i.
/// Ports `HyperIdealFunctional.σij()`.
inline double sigma_ij(double aij, double bi, double bj, bool vj_var)
{
    if (vj_var) return zeta13(aij, bi, bj);
    return zeta14(-aij, bi);
}

/// `α_ij`: computed dihedral angle at edge ij in the face with vertices i, j, k.
///
/// Arguments (cyclic role assignment):
///   * `aij, ajk, aki` — edge variables.
///   * `bi, bj, bk`    — vertex variables.
///   * `beta_i, beta_j, beta_k` — interior angles of the auxiliary hyperbolic triangle.
///   * `vi_var, vj_var, vk_var` — which vertices are hyper-ideal.
///
/// Ports `HyperIdealFunctional.αij()` (the private helper).
/// Note: the `vk_var` case recurses once (never more than one level deep).
inline double alpha_ij(
    double aij, double ajk, double aki,
    double bi,  double bj,  double bk,
    double beta_i, double beta_j, double beta_k,
    bool vi_var, bool vj_var, bool vk_var)
{
    if (vi_var) {
        double si  = sigma_i (aij, aki, ajk, vj_var, vk_var);
        double sij = sigma_ij(aij, bi, bj, vj_var);
        double sik = sigma_ij(aki, bi, bk, vk_var);
        return zeta(si, sij, sik);
    }
    if (vj_var) {
        double sj  = sigma_i (ajk, aij, aki, vk_var, vi_var);
        double sjk = sigma_ij(ajk, bj, bk, vk_var);
        double sji = sigma_ij(aij, bj, bi, vi_var);
        return zeta(sj, sji, sjk);
    }
    if (vk_var) {
        // Derive α_ij from α_jk (one level of recursion).
        double a_jk = alpha_ij(ajk, aki, aij,
                               bj, bk, bi,
                               beta_j, beta_k, beta_i,
                               vj_var, vk_var, vi_var);
        return PI - a_jk - beta_j;
    }
    // All ideal: closed-form formula.
    return 0.5 * (PI + beta_k - beta_i - beta_j);
}

} // namespace conformallab