ConformalLabpp/code/include/euclidean_geometry.hpp

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

// euclidean_geometry.hpp
//
// Corner-angle formula for Euclidean triangles in the discrete conformal
// (log-length) parametrisation.
//
// Ported from de.varylab.discreteconformal.functional.EuclideanCyclicFunctional.
//
// In the discrete conformal parametrisation a Euclidean triangle is described by
// its three effective log-lengths Λ̃_ij = λ°_ij + u_i + u_j (+ edge DOF).
// The corresponding side lengths are l_ij = exp(Λ̃_ij / 2).
//
// Vertex ordering convention (matches EuclideanCyclicFunctional.java):
//   v1 is opposite edge l23,  v2 is opposite l31,  v3 is opposite l12.
//
// t-value trick (Springborn 2008 §3):
//   t12 = −l12 + l23 + l31  =  2(s − l12)
//   t23 = +l12 − l23 + l31  =  2(s − l23)
//   t31 = +l12 + l23 − l31  =  2(s − l31)
//   denom = sqrt(t12 · t23 · t31 · l123) = 4 · Area
//
//   α_v = 2 · atan2( product of t-values adjacent to v, denom )
//
// The centering trick (l_ij ← exp((Λ̃_ij − 2·μ)/2), μ = (Λ̃12+Λ̃23+Λ̃31)/6)
// rescales all three sides by the same factor, leaving angles unchanged but
// keeping the arguments of exp in a safe numerical range.

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

namespace conformallab {

/// Interior corner angles of a Euclidean triangle.
struct EuclideanFaceAngles {
    double alpha1;   ///< Corner angle at v₁ (opposite l₂₃).
    double alpha2;   ///< Corner angle at v₂ (opposite l₃₁).
    double alpha3;   ///< Corner angle at v₃ (opposite l₁₂).
    bool   valid;    ///< `false` when the triangle is degenerate.
};

/// Compute the corner angles of a Euclidean triangle from its three
/// side lengths.  Returns `valid = false` when the triangle inequality
/// is violated.
inline EuclideanFaceAngles euclidean_angles_from_lengths(
    double l12, double l23, double l31)
{
    const double t12 = -l12 + l23 + l31;   // 2*(s − l12)
    const double t23 = +l12 - l23 + l31;   // 2*(s − l23)
    const double t31 = +l12 + l23 - l31;   // 2*(s − l31)

    // Degenerate (triangle inequality violated): return the *limiting* angles
    // of the flat-out triangle — the corner opposite the over-long edge is π,
    // the other two are 0.  This is the convex C¹ extension of the BPS energy
    // onto the infeasible region and is what the Java reference
    // (EuclideanCyclicFunctional.triangleEnergyAndAlphas) does.  `valid` stays
    // false so the cotangent Hessian still skips this face.  At most one t can
    // be ≤ 0 (t12+t23 = 2·l31 > 0, etc.), so the order of these checks is moot.
    if (t23 <= 0.0) return {PI, 0.0, 0.0, false};   // l23 too long → α₁ = π
    if (t31 <= 0.0) return {0.0, PI, 0.0, false};   // l31 too long → α₂ = π
    if (t12 <= 0.0) return {0.0, 0.0, PI, false};   // l12 too long → α₃ = π

    const double l123  = l12 + l23 + l31;
    const double denom2 = t12 * t23 * t31 * l123;  // = (4·Area)²
    if (denom2 <= 0.0)
        return {0.0, 0.0, 0.0, false};

    const double denom = std::sqrt(denom2);

    // α at v1 (opposite l23): adjacent t-values are t12 and t31
    // α at v2 (opposite l31): adjacent t-values are t12 and t23
    // α at v3 (opposite l12): adjacent t-values are t23 and t31
    return {
        2.0 * std::atan2(t12 * t31, denom),
        2.0 * std::atan2(t12 * t23, denom),
        2.0 * std::atan2(t23 * t31, denom),
        true
    };
}

/// Compute the corner angles of a Euclidean triangle from its three
/// effective log-lengths `Λ̃ᵢⱼ`.  Internally centres lengths so that
/// `l₁₂·l₂₃·l₃₁ = 1` to avoid float overflow for large `|Λ̃|`.
inline EuclideanFaceAngles euclidean_angles(
    double lam12, double lam23, double lam31)
{
    const double mu   = (lam12 + lam23 + lam31) / 6.0;
    const double l12  = std::exp((lam12 - 2.0 * mu) * 0.5);
    const double l23  = std::exp((lam23 - 2.0 * mu) * 0.5);
    const double l31  = std::exp((lam31 - 2.0 * mu) * 0.5);
    return euclidean_angles_from_lengths(l12, l23, l31);
}

} // namespace conformallab