docs(doxygen): 100% public-API coverage (228 → 0 undocumented)

Completes the work begun in the previous commit on this branch.  Every
public symbol under code/include/ now carries a brief Doxygen comment
(0 undocumented per scripts/doxygen-coverage.sh, with the `detail::`
implementation namespaces excluded as before).

Trajectory on this branch:
  start (after Doxyfile fix):  24.0 %  (165 / 437 in the no-detail set
                                       was 105 / 437 when detail counted)
  after PR #17 base commit  :  42.4 %  (165 / 396)
  this commit               : 100.0 %  (396 / 396)

Files touched (all .hpp / .h headers under code/include/):
  * cgal/Conformal_map_traits.h
  * clausen.hpp, conformal_mesh.hpp, constants.hpp (already docd)
  * cp_euclidean_functional.hpp, cut_graph.hpp, discrete_elliptic_utility.hpp
  * euclidean_functional.hpp, euclidean_geometry.hpp, euclidean_hessian.hpp
  * fundamental_domain.hpp, gauss_bonnet.hpp
  * hyper_ideal_{functional,geometry,hessian,utility,visualization_utility}.hpp
  * inversive_distance_functional.hpp, layout.hpp
  * matrix_utility.hpp, mesh_builder.hpp, mesh_io.hpp
  * newton_solver.hpp, p2_utility.hpp, period_matrix.hpp, projective_math.hpp
  * serialization.hpp, spherical_functional.hpp, spherical_geometry.hpp
  * spherical_hessian.hpp, viewer_utils.h

CI:
.gitea/workflows/doxygen-pages.yml now enforces
`scripts/doxygen-coverage.sh --threshold 100`, so any future regression
(a new public function landed without a `///` brief) fails the build
before the Doxygen HTML is published to Codeberg Pages.

Doxygen warnings remain at 0.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

.gitea/workflows/doxygen-pages.yml

@@ -47,8 +47,12 @@ jobs:
             head -30 doc/doxygen/doxygen-warnings.log
           fi
 
-      - name: Report Doxygen coverage
-        run: bash scripts/doxygen-coverage.sh
+      - name: Enforce Doxygen coverage 100%
+        # Coverage is measured against every public symbol under
+        # code/include/ (the `detail::` namespaces are excluded).  As of
+        # the `docs/doxygen-coverage-100` PR the baseline is 100 %, so
+        # the gate fires only on regressions.
+        run: bash scripts/doxygen-coverage.sh --threshold 100
 
       - name: Regenerate doc/api/headers.md from XML
         run: |

code/include/CGAL/Conformal_map_traits.h

@@ -80,10 +80,13 @@ namespace CGAL {
 // Default_conformal_map_traits  —  primary template (undefined)
 // ════════════════════════════════════════════════════════════════════════════
 //
-// The undefined primary template forces specialisation per mesh type.
-// MVP provides only the `Surface_mesh` specialisation below; further
-// mesh types (Polyhedron_3, OpenMesh, pmp) are deferred to Phase 8a.2.
-
+/*!
+\ingroup PkgConformalMapConcepts
+\brief Primary `ConformalMapTraits` template — undefined, must be
+specialised per mesh type.  The MVP only ships the `Surface_mesh`
+specialisation below; further mesh types (Polyhedron_3, OpenMesh,
+pmp) are deferred to Phase 8a.2.
+*/
 template <typename TriangleMesh,
           typename Kernel_ = CGAL::Simple_cartesian<double>>
 struct Default_conformal_map_traits;

code/include/clausen.hpp

@@ -128,9 +128,9 @@ inline int ncl5pi6() noexcept {
 
 } // namespace detail
 
-// Clausen's integral Cl2(x) = -integral_0^x log|2 sin(t/2)| dt.
-// High-precision Chebyshev implementation.
-// Corresponds to Java Clausen.clausen2().
+/// Clausen's integral `Cl₂(x) = −∫₀ˣ log|2 sin(t/2)| dt`,
+/// computed via a high-precision Chebyshev expansion.
+/// Same as Java `Clausen.clausen2()`.
 inline double clausen2(double x) noexcept {
     using namespace detail;
     constexpr double pi  = 3.14159265358979323846264338328;
@@ -154,8 +154,8 @@ inline double clausen2(double x) noexcept {
     return rh ? -f : f;
 }
 
-// Milnor's Lobachevsky function Л(x) = Cl2(2x)/2.
-// Corresponds to Java Clausen.Л().
+/// Milnor's Lobachevsky function `Л(x) = Cl₂(2x) / 2`.
+/// Same as Java `Clausen.Л()`.
 inline double Lobachevsky(double x) noexcept {
     constexpr double pi = 3.14159265358979323846264338328;
     x = std::fmod(x, pi);
@@ -163,8 +163,8 @@ inline double Lobachevsky(double x) noexcept {
     return clausen2(2.0 * x) / 2.0;
 }
 
-// Imaginary part of the dilogarithm Im(Li2(z)).
-// Corresponds to Java Clausen.ImLi2().
+/// Imaginary part of the dilogarithm `Im(Li₂(z))` for complex `z`.
+/// Same as Java `Clausen.ImLi2()`.
 inline double ImLi2(std::complex<double> z) noexcept {
     auto a = std::log(1.0 - std::conj(z));  // log(1 - conj(z))
     auto b = std::log(1.0 - z);             // log(1 - z)

code/include/cp_euclidean_functional.hpp

@@ -77,12 +77,17 @@ namespace conformallab {
 
 // ── Property-map type aliases ────────────────────────────────────────────────
 
+/// Property map face → `int`    for the CP-Euclidean functional.
 using CPFMapI = ConformalMesh::Property_map<Face_index, int>;
+/// Property map face → `double` for the CP-Euclidean functional.
 using CPFMapD = ConformalMesh::Property_map<Face_index, double>;
+/// Property map edge → `double` for the CP-Euclidean functional.
 using CPEMapD = ConformalMesh::Property_map<Edge_index, double>;
 
 // ── Persistent map bundle ─────────────────────────────────────────────────────
 
+/// Bundle of the three property maps consumed by the CP-Euclidean
+/// (Bobenko-Pinkall-Springborn 2010) circle-packing functional.
 struct CPEuclideanMaps {
     CPFMapI f_idx;     ///< DOF index per face (−1 = pinned)
     CPEMapD theta_e;   ///< intersection angle per edge (default π/2 = orthogonal)
@@ -184,11 +189,8 @@ inline double dof_val(int idx, const std::vector<double>& x) noexcept
 
 } // namespace cp_detail
 
-// ── Energy ────────────────────────────────────────────────────────────────────
-//
-// Mirrors evaluateEnergyAndGradient in the Java code (lines 170-240) for the
-// energy accumulation only.  The gradient is computed in a dedicated function
-// below for clarity.
+/// CP-Euclidean energy value at DOF vector `x` (ρ per face).
+/// Mirrors `evaluateEnergyAndGradient()` in the Java original (lines 170-240).
 inline double cp_euclidean_energy(const ConformalMesh&        mesh,
                                     const std::vector<double>&  x,
                                     const CPEuclideanMaps&      m)
@@ -233,10 +235,8 @@ inline double cp_euclidean_energy(const ConformalMesh&        mesh,
     return E;
 }
 
-// ── Gradient ──────────────────────────────────────────────────────────────────
-//
-// ∂E/∂ρ_f = φ_f − Σ_{h: face(h)=f, !is_border(h)}  (p + θ*)
-//         OR (boundary):                            2 θ*
+/// CP-Euclidean gradient `∂E/∂ρ_f` (per face DOF).  Interior term
+/// `−(p + θ*)`, boundary term `−2 θ*`; see `setup_cp_euclidean_maps`.
 inline std::vector<double> cp_euclidean_gradient(const ConformalMesh&        mesh,
                                                    const std::vector<double>&  x,
                                                    const CPEuclideanMaps&      m)
@@ -280,12 +280,10 @@ inline std::vector<double> cp_euclidean_gradient(const ConformalMesh&        mes
     return G;
 }
 
-// ── Hessian (analytic) ────────────────────────────────────────────────────────
-//
-// Per interior undirected edge e with adjacent faces (j, k):
-//   h_jk = sin θ / (cosh(Δρ) − cos θ)
-// Diagonal contributions on both endpoints; off-diagonal block is −h_jk.
-// Pinned faces are skipped (their DOF index is −1 ⇒ excluded from the matrix).
+/// Analytic CP-Euclidean Hessian, sparse.  Per interior edge `(j,k)`
+/// the contribution is `h_jk = sin θ / (cosh(Δρ) − cos θ)`, added to
+/// diagonals `H_jj`, `H_kk` and subtracted off-diagonals `H_jk = H_kj`.
+/// Pinned faces are excluded (DOF index −1).
 inline Eigen::SparseMatrix<double> cp_euclidean_hessian(const ConformalMesh&        mesh,
                                                           const std::vector<double>&  x,
                                                           const CPEuclideanMaps&      m)
@@ -323,11 +321,8 @@ inline Eigen::SparseMatrix<double> cp_euclidean_hessian(const ConformalMesh&
     return H;
 }
 
-// ── Finite-difference gradient check ─────────────────────────────────────────
-//
-// Mirrors the Java FunctionalTest pattern:
-//   For each DOF i, compare analytic G[i] to (E(x+ε·e_i) − E(x−ε·e_i)) / (2ε).
-// Default tolerance 1e-6 with ε = 1e-5 leaves ~3 digits of margin for sane meshes.
+/// FD gradient check for the CP-Euclidean functional.  Mirrors the
+/// Java `FunctionalTest`; default `eps = 1e-5`, `tol = 1e-6`.
 inline bool gradient_check_cp_euclidean(const ConformalMesh&        mesh,
                                           const std::vector<double>&  x,
                                           const CPEuclideanMaps&      m,
@@ -355,10 +350,8 @@ inline bool gradient_check_cp_euclidean(const ConformalMesh&        mesh,
     return true;
 }
 
-// ── Finite-difference Hessian check ──────────────────────────────────────────
-//
-// Verifies analytic H against ( G(x+ε·e_i) − G(x−ε·e_i) ) / (2ε) column-wise.
-// Symmetry is implicit in the analytic form; we check both off-diagonal entries.
+/// FD Hessian check for the CP-Euclidean functional.  Verifies analytic
+/// `H` column-by-column against `(G(x+εe_j) − G(x−εe_j)) / (2ε)`.
 inline bool hessian_check_cp_euclidean(const ConformalMesh&        mesh,
                                          const std::vector<double>&  x,
                                          const CPEuclideanMaps&      m,

code/include/cut_graph.hpp

@@ -36,6 +36,8 @@ namespace conformallab {
 // CutGraph
 // ─────────────────────────────────────────────────────────────────────────────
 
+/// Cut-graph result of the tree-cotree algorithm: the set of `2g` edges
+/// whose removal turns a closed genus-`g` surface into a topological disk.
 struct CutGraph {
     /// cut_edge_flags[e.idx()] = true ↔ this edge is a cut edge.
     /// Size = mesh.number_of_edges().
@@ -47,6 +49,7 @@ struct CutGraph {
     /// Genus of the surface (0 for topological spheres and open patches).
     int genus = 0;
 
+    /// `true` iff edge `e` is a cut edge of this graph.
     bool is_cut(Edge_index e) const
     {
         return static_cast<std::size_t>(e.idx()) < cut_edge_flags.size()
@@ -54,17 +57,10 @@ struct CutGraph {
     }
 };
 
-// ─────────────────────────────────────────────────────────────────────────────
-// compute_cut_graph
-// ─────────────────────────────────────────────────────────────────────────────
-//
-// Implements the standard tree-cotree algorithm (Erickson–Whittlesey 2005):
-//
-//   Step 1: BFS primal spanning tree T  (V−1 primal tree edges).
-//   Step 2: BFS dual spanning tree T*   (F−1 dual/primal edges, avoiding
-//           edges whose primal crosses T).
-//   Step 3: cut edges = primal edges neither in T nor "used" by T*.
-
+/// Compute the cut graph of `mesh` via the standard tree-cotree
+/// algorithm (Erickson–Whittlesey 2005): primal BFS spanning tree T,
+/// dual BFS spanning tree T* avoiding T-primals, then the `2g` cut
+/// edges are those in neither T nor T*.
 inline CutGraph compute_cut_graph(const ConformalMesh& mesh)
 {
     const std::size_t nv = mesh.number_of_vertices();

code/include/discrete_elliptic_utility.hpp

@@ -17,7 +17,9 @@ namespace conformallab {
 //   3. Re-flip:   Re < 0              →  Re = -Re
 //   4. S-invert:  |tau| < 1           →  tau = 1/tau
 //
-// Corresponds to Java DiscreteEllipticUtility.normalizeModulus(Complex).
+/// Normalise a complex modulus `τ` into the standard fundamental
+/// domain of an elliptic curve (`|τ| ≥ 1`, `0 ≤ Re τ ≤ ½`, `Im τ ≥ 0`).
+/// Same as Java `DiscreteEllipticUtility.normalizeModulus(Complex)`.
 inline std::complex<double> normalizeModulus(std::complex<double> tau) {
     int maxIter = 100;
     while (--maxIter > 0) {

code/include/euclidean_functional.hpp

@@ -48,13 +48,18 @@ namespace conformallab {
 
 // ── Property-map type aliases ─────────────────────────────────────────────────
 
+/// Property map vertex → `double` for the Euclidean functional.
 using EuclVMapD = ConformalMesh::Property_map<Vertex_index, double>;
+/// Property map vertex → `int`    for the Euclidean functional.
 using EuclVMapI = ConformalMesh::Property_map<Vertex_index, int>;
+/// Property map edge   → `double` for the Euclidean functional.
 using EuclEMapD = ConformalMesh::Property_map<Edge_index,   double>;
+/// Property map edge   → `int`    for the Euclidean functional.
 using EuclEMapI = ConformalMesh::Property_map<Edge_index,   int>;
 
 // ── Persistent map bundle ─────────────────────────────────────────────────────
 
+/// Bundle of the five property maps consumed by the Euclidean functional.
 struct EuclideanMaps {
     EuclVMapI v_idx;    ///< DOF index per vertex   (-1 = pinned / u_v = 0)
     EuclEMapI e_idx;    ///< DOF index per edge      (-1 = no edge DOF)
@@ -119,10 +124,9 @@ inline int euclidean_dimension(const ConformalMesh& mesh, const EuclideanMaps& m
     return dim;
 }
 
-// Set lambda0 from mesh vertex positions (Euclidean):
-//   λ°_e = 2·log(|p_i − p_j|)   (natural log of Euclidean edge length squared)
-//
-// This gives exp(Λ̃_ij / 2) = l_ij at x=0.
+/// Set `lambda0` from mesh vertex positions:
+/// `λ°_e = 2·log(|p_i − p_j|)` (natural log of Euclidean edge length²).
+/// This gives `exp(Λ̃_ij / 2) = l_ij` at `x = 0`.
 inline void compute_euclidean_lambda0_from_mesh(ConformalMesh& mesh, EuclideanMaps& m)
 {
     for (auto e : mesh.edges()) {
@@ -142,25 +146,23 @@ inline void compute_euclidean_lambda0_from_mesh(ConformalMesh& mesh, EuclideanMa
 
 // ── Internal helpers ──────────────────────────────────────────────────────────
 
+/// Read DOF value from `x` for index `idx`; return 0 if pinned (idx < 0).
 static inline double eucl_dof_val(int idx, const std::vector<double>& x)
 {
     return idx >= 0 ? x[static_cast<std::size_t>(idx)] : 0.0;
 }
 
+/// Convert a CGAL half-edge index to a plain `std::size_t` for vector indexing.
 static inline std::size_t eucl_hidx(Halfedge_index h)
 {
     return static_cast<std::size_t>(static_cast<std::uint32_t>(h));
 }
 
-// ── Gradient ──────────────────────────────────────────────────────────────────
-//
-// G_v = Θ_v − Σ_{faces adj. v} α_v(face)
-// G_e = α_opp(face⁺) + α_opp(face⁻) − φ_e
-//
-// Corner-angle storage (h_alpha):
-//   h_alpha[h] = corner angle OPPOSITE to the edge of halfedge h in its face.
-//   h_alpha[h0] = α3,  h_alpha[h1] = α1,  h_alpha[h2] = α2
-// (same convention as SphericalFunctional)
+/// Compute the Euclidean-functional gradient G(x):
+///   * `G_v = Θ_v − Σ_faces α_v(face)`
+///   * `G_e = α_opp(face⁺) + α_opp(face⁻) − φ_e`
+///
+/// Same half-edge corner-angle storage convention as `spherical_gradient`.
 inline std::vector<double> euclidean_gradient(
     ConformalMesh&             mesh,
     const std::vector<double>& x,
@@ -238,9 +240,8 @@ inline std::vector<double> euclidean_gradient(
     return G;
 }
 
-// ── Energy via Gauss-Legendre path integral ───────────────────────────────────
-//
-// E(x) = ∫₀¹ ⟨G(tx), x⟩ dt   (10-point GL quadrature, same as SphericalFunctional)
+/// Euclidean energy `E(x) = ∫₀¹ ⟨G(t·x), x⟩ dt`, evaluated with
+/// 10-point Gauss-Legendre quadrature (same as the Spherical functional).
 inline double euclidean_energy(
     ConformalMesh&             mesh,
     const std::vector<double>& x,
@@ -282,11 +283,14 @@ inline double euclidean_energy(
 
 // ── Full evaluation (energy + gradient) ──────────────────────────────────────
 
+/// Output of `evaluate_euclidean()` — energy plus optional gradient.
 struct EuclideanResult {
-    double              energy   = 0.0;
-    std::vector<double> gradient;
+    double              energy   = 0.0; ///< Functional value at input DOFs.
+    std::vector<double> gradient;       ///< Gradient ∇E (empty if not requested).
 };
 
+/// Evaluate the Euclidean functional at DOFs `x`.  Returns energy and
+/// gradient (toggle via `need_energy` / `need_gradient`).
 inline EuclideanResult evaluate_euclidean(
     ConformalMesh&             mesh,
     const std::vector<double>& x,
@@ -302,10 +306,8 @@ inline EuclideanResult evaluate_euclidean(
     return res;
 }
 
-// ── Finite-difference gradient check ─────────────────────────────────────────
-//
-// Tests |G[i] − (E(x+εeᵢ) − E(x−εeᵢ))/(2ε)| / max(1,|G[i]|) < tol
-// for all variable DOFs.
+/// Finite-difference gradient check for the Euclidean functional
+/// (central differences).  Defaults `eps = 1e-5`, `tol = 1e-4`.
 inline bool gradient_check_euclidean(
     ConformalMesh&             mesh,
     const std::vector<double>& x0,

code/include/euclidean_geometry.hpp

@@ -29,19 +29,17 @@
 
 namespace conformallab {
 
+/// Interior corner angles of a Euclidean triangle.
 struct EuclideanFaceAngles {
-    double alpha1;   ///< corner angle at v1 (opposite l23)
-    double alpha2;   ///< corner angle at v2 (opposite l31)
-    double alpha3;   ///< corner angle at v3 (opposite l12)
-    bool   valid;
+    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.
 };
 
-// ── From side lengths ─────────────────────────────────────────────────────────
-//
-// Given three Euclidean side lengths l12, l23, l31 > 0 satisfying the triangle
-// inequality, compute the corner angles.
-//
-// Returns valid=false if the triangle inequality is violated (any t-value ≤ 0).
+/// 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)
 {
@@ -70,14 +68,9 @@ inline EuclideanFaceAngles euclidean_angles_from_lengths(
     };
 }
 
-// ── From effective log-lengths Λ̃ ─────────────────────────────────────────────
-//
-// Converts to side lengths l_ij = exp(Λ̃_ij / 2), applying the centering
-// trick for numerical safety, then delegates to euclidean_angles_from_lengths.
-//
-// The centering constant μ = (Λ̃12 + Λ̃23 + Λ̃31) / 6 ensures
-//   l12 · l23 · l31 = 1  (geometric mean = 1)
-// which keeps all l values near 1 and prevents float overflow for large |Λ̃|.
+/// 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)
 {

code/include/euclidean_hessian.hpp

@@ -52,8 +52,18 @@ namespace conformallab {
 // cot_k = (t_adj·l123 − t_opp·t_other) / (8·Area)
 //
 // Returns {0,0,0} for degenerate faces (triangle inequality violated or Area=0).
-struct EuclCotWeights { double cot1, cot2, cot3; bool valid; };
+/// Three Euclidean cotangent weights `(cot1, cot2, cot3)` for the
+/// vertices opposite to edges (l₂₃, l₃₁, l₁₂) of a triangle, plus a
+/// `valid` flag that is `false` when the triangle is degenerate.
+struct EuclCotWeights {
+    double cot1; ///< Cotangent at vertex 1 (opposite to l₂₃).
+    double cot2; ///< Cotangent at vertex 2 (opposite to l₃₁).
+    double cot3; ///< Cotangent at vertex 3 (opposite to l₁₂).
+    bool   valid;///< `false` when the triangle is degenerate (triangle inequality violated or area = 0).
+};
 
+/// Compute the three Euclidean cotangent weights from edge lengths.
+/// Returns `{0,0,0,false}` for degenerate triangles.
 inline EuclCotWeights euclidean_cot_weights(double l12, double l23, double l31)
 {
     const double t12 = -l12 + l23 + l31;
@@ -81,15 +91,9 @@ inline EuclCotWeights euclidean_cot_weights(double l12, double l23, double l31)
     };
 }
 
-// ── Analytical Hessian (cotangent Laplacian) ──────────────────────────────────
-//
-// Returns the n×n sparse Hessian matrix H where n = euclidean_dimension(mesh, m).
-//
-// Only vertex DOFs are supported.  Edge DOFs (m.e_idx[e] >= 0) produce
-// additional mixed-derivative entries that are not yet implemented; this
-// function asserts they are absent.
-//
-// x   – current DOF vector (used to compute effective log-lengths Λ̃ij).
+/// Analytical Euclidean Hessian (cotangent Laplacian), sparse.
+/// Only vertex DOFs are supported — the function asserts that no edge
+/// DOF is variable.  `x` is used to compute effective log-lengths Λ̃ᵢⱼ.
 inline Eigen::SparseMatrix<double> euclidean_hessian(
     ConformalMesh&             mesh,
     const std::vector<double>& x,
@@ -168,11 +172,9 @@ inline Eigen::SparseMatrix<double> euclidean_hessian(
 }
 
 // ── Finite-difference Hessian check ──────────────────────────────────────────
-//
-// Compares the analytical Hessian column-by-column against
-//   H_fd[:, j] = (G(x + ε·eⱼ) − G(x − ε·eⱼ)) / (2ε).
-//
-// Returns true if max relative error < tol for every entry.
+/// FD Hessian check for the Euclidean functional.  Compares analytic
+/// `H` column-by-column to `(G(x+εeⱼ) − G(x−εeⱼ)) / (2ε)`; returns
+/// `true` iff max relative error is below `tol`.
 inline bool hessian_check_euclidean(
     ConformalMesh&             mesh,
     const std::vector<double>& x0,

code/include/fundamental_domain.hpp

@@ -43,6 +43,9 @@ namespace conformallab {
 // FundamentalDomain
 // ─────────────────────────────────────────────────────────────────────────────
 
+/// Fundamental polygon of a closed surface obtained by cutting along a
+/// `CutGraph`: corner vertices, paired-edge identifications and holonomy
+/// generators.  For genus-1 the polygon is a parallelogram with 4 corners.
 struct FundamentalDomain {
     /// Polygon corners in order (CCW).  Size = 4 for genus-1.
     std::vector<Eigen::Vector2d> vertices;
@@ -56,6 +59,7 @@ struct FundamentalDomain {
     /// For genus-1: generators[0] = ω_1, generators[1] = ω_2.
     std::vector<Eigen::Vector2d> generators;
 
+    /// `true` iff the polygon has at least 3 vertices.
     bool is_valid() const { return vertices.size() >= 3; }
 };
 
@@ -75,6 +79,8 @@ struct FundamentalDomain {
 //   bottom (v0→v1) ≡ top    (v3→v2)   by ω_2
 //   left   (v3→v0) ≡ right  (v2→v1)   by ω_1   (reversed convention)
 // ─────────────────────────────────────────────────────────────────────────────
+/// Build the parallelogram fundamental domain from genus-1 Euclidean
+/// holonomy data (`hol.translations[0] = ω₁`, `hol.translations[1] = ω₂`).
 inline FundamentalDomain compute_fundamental_domain_genus1(
     const HolonomyData& hol)
 {
@@ -148,6 +154,8 @@ inline FundamentalDomain compute_fundamental_domain_genus1(
 //   this is intentionally deferred and NOT implemented here.
 //   See period_matrix.hpp for the genus-1 case (τ = ω_2/ω_1 ∈ ℍ).
 // ─────────────────────────────────────────────────────────────────────────────
+/// Dispatcher: for genus 1 returns `compute_fundamental_domain_genus1`,
+/// for higher genus returns an empty domain (4g-polygon not yet implemented).
 inline FundamentalDomain compute_fundamental_domain(
     const HolonomyData& hol)
 {
@@ -165,6 +173,8 @@ inline FundamentalDomain compute_fundamental_domain(
 // return a translated copy of the layout shifted by m·ω_1 + n·ω_2.
 // Useful for visualising the tiled universal cover.
 // ─────────────────────────────────────────────────────────────────────────────
+/// Return a translated copy of `layout` shifted by `m·ω₁ + n·ω₂`.
+/// Useful for visualising the tiled universal cover.
 inline Layout2D tiling_copy(const Layout2D& layout,
                              const Eigen::Vector2d& w1,
                              const Eigen::Vector2d& w2,
@@ -183,6 +193,8 @@ inline Layout2D tiling_copy(const Layout2D& layout,
 // Returns a vector of tiling copies for (m, n) with |m| ≤ m_max, |n| ≤ n_max.
 // The result includes the original (m=0, n=0) at index (m_max)(2*n_max+1)+n_max.
 // ─────────────────────────────────────────────────────────────────────────────
+/// Build a tiling neighbourhood: all `(m, n)` with `|m| ≤ m_max`,
+/// `|n| ≤ n_max`.  The original tile `(0, 0)` is included.
 inline std::vector<Layout2D> tiling_neighbourhood(
     const Layout2D&    layout,
     const HolonomyData& hol,

code/include/gauss_bonnet.hpp

@@ -56,6 +56,7 @@ inline int genus(const ConformalMesh& mesh)
 
 // ── Left-hand side  Σ(2π − Θ_v) ─────────────────────────────────────────────
 
+/// Sum `Σ_v (2π − Θ_v)` for a raw vertex → angle property map.
 inline double gauss_bonnet_sum(
     const ConformalMesh&                                      mesh,
     const ConformalMesh::Property_map<Vertex_index, double>&  theta)
@@ -66,15 +67,19 @@ inline double gauss_bonnet_sum(
     return s;
 }
 
+/// `gauss_bonnet_sum` for the Euclidean-functional property bundle.
 inline double gauss_bonnet_sum(const ConformalMesh& m, const EuclideanMaps&   mp)
     { return gauss_bonnet_sum(m, mp.theta_v); }
+/// `gauss_bonnet_sum` for the Spherical-functional property bundle.
 inline double gauss_bonnet_sum(const ConformalMesh& m, const SphericalMaps&   mp)
     { return gauss_bonnet_sum(m, mp.theta_v); }
+/// `gauss_bonnet_sum` for the HyperIdeal-functional property bundle.
 inline double gauss_bonnet_sum(const ConformalMesh& m, const HyperIdealMaps&  mp)
     { return gauss_bonnet_sum(m, mp.theta_v); }
 
 // ── Right-hand side  2π · χ(M) ───────────────────────────────────────────────
 
+/// Right-hand side of Gauss-Bonnet: `2π · χ(M)`.
 inline double gauss_bonnet_rhs(const ConformalMesh& mesh)
 {
     return TWO_PI * static_cast<double>(euler_characteristic(mesh));
@@ -82,14 +87,15 @@ inline double gauss_bonnet_rhs(const ConformalMesh& mesh)
 
 // ── Deficit: lhs − rhs  (0 = Gauss–Bonnet satisfied) ─────────────────────────
 
+/// Gauss-Bonnet deficit `lhs − rhs`; zero iff the identity is satisfied.
 template <typename Maps>
 inline double gauss_bonnet_deficit(const ConformalMesh& mesh, const Maps& maps)
 {
     return gauss_bonnet_sum(mesh, maps) - gauss_bonnet_rhs(mesh);
 }
 
-// ── check_gauss_bonnet — throws std::runtime_error if |deficit| > tol ─────────
-
+/// Throws `std::runtime_error` if `|lhs − 2π·χ| > tol`.
+/// Overload accepting a precomputed `lhs`.
 inline void check_gauss_bonnet(const ConformalMesh& mesh,
                                 double lhs,
                                 double tol = 1e-8)
@@ -108,6 +114,7 @@ inline void check_gauss_bonnet(const ConformalMesh& mesh,
     }
 }
 
+/// Throws `std::runtime_error` if Gauss-Bonnet is violated by more than `tol`.
 template <typename Maps>
 inline void check_gauss_bonnet(const ConformalMesh& mesh,
                                 const Maps& maps,
@@ -123,6 +130,9 @@ inline void check_gauss_bonnet(const ConformalMesh& mesh,
 // Only modifies free vertices (v_idx[v] >= 0 for EuclideanMaps / SphericalMaps;
 // always all vertices for the raw property-map overload).
 
+/// Distribute the Gauss-Bonnet deficit uniformly across all `Θ_v`:
+/// add `δ = (lhs − rhs) / V` to every entry so that the identity holds
+/// exactly afterwards.  Overload for a raw property map.
 inline void enforce_gauss_bonnet(
     ConformalMesh&                                           mesh,
     ConformalMesh::Property_map<Vertex_index, double>&       theta)
@@ -136,6 +146,7 @@ inline void enforce_gauss_bonnet(
         theta[v] += delta;
 }
 
+/// Distribute the Gauss-Bonnet deficit uniformly across `maps.theta_v`.
 template <typename Maps>
 inline void enforce_gauss_bonnet(ConformalMesh& mesh, Maps& maps)
 {

code/include/hyper_ideal_functional.hpp

@@ -39,18 +39,23 @@ namespace conformallab {
 
 // ── Property-map type aliases ─────────────────────────────────────────────────
 
+/// Property map vertex → `double` (HyperIdeal scalar-per-vertex data).
 using VMapD = ConformalMesh::Property_map<Vertex_index,   double>;
+/// Property map vertex → `int`    (HyperIdeal DOF indices).
 using VMapI = ConformalMesh::Property_map<Vertex_index,   int>;
+/// Property map edge   → `double` (HyperIdeal scalar-per-edge data).
 using EMapD = ConformalMesh::Property_map<Edge_index,     double>;
+/// Property map edge   → `int`    (HyperIdeal DOF indices).
 using EMapI = ConformalMesh::Property_map<Edge_index,     int>;
 
 // ── Persistent map bundle ─────────────────────────────────────────────────────
 
+/// Bundle of the four property maps consumed by the HyperIdeal functional.
 struct HyperIdealMaps {
-    VMapI v_idx;    // DOF index per vertex    (-1 = pinned / ideal point)
-    EMapI e_idx;    // DOF index per edge      (-1 = fixed)
-    VMapD theta_v;  // target cone angle Θ_v   (parameter, not variable)
-    EMapD theta_e;  // target intersection angle θ_e
+    VMapI v_idx;    ///< DOF index per vertex (−1 = pinned / ideal point).
+    EMapI e_idx;    ///< DOF index per edge   (−1 = fixed).
+    VMapD theta_v;  ///< Target cone angle Θᵥ (parameter, not variable).
+    EMapD theta_e;  ///< Target intersection angle θₑ.
 };
 
 /// Attach the four HyperIdeal property maps to `mesh` and return their
@@ -105,20 +110,22 @@ inline int assign_all_dof_indices(ConformalMesh& mesh, HyperIdealMaps& m)
 
 // ── Evaluation result ─────────────────────────────────────────────────────────
 
+/// Output of `evaluate_hyper_ideal()` — the energy value and (optionally)
+/// its gradient evaluated at the current DOF vector.
 struct HyperIdealResult {
-    double              energy   = 0.0;
-    std::vector<double> gradient;     // empty when gradient was not requested
+    double              energy   = 0.0; ///< Functional value at the input DOFs.
+    std::vector<double> gradient;       ///< Gradient ∇E; empty when not requested.
 };
 
 // ── Internal helpers ──────────────────────────────────────────────────────────
 
-// Get the DOF value from x, or 0.0 if pinned.
+/// Read the DOF value from `x` for index `idx`; return 0 if pinned (idx < 0).
 static inline double dof_val(int idx, const std::vector<double>& x)
 {
     return idx >= 0 ? x[static_cast<std::size_t>(idx)] : 0.0;
 }
 
-// Convert a CGAL halfedge index to a plain std::size_t (for vector indexing).
+/// Convert a CGAL half-edge index to a plain `std::size_t` for vector indexing.
 static inline std::size_t hidx(Halfedge_index h)
 {
     return static_cast<std::size_t>(static_cast<std::uint32_t>(h));
@@ -144,11 +151,21 @@ static inline std::size_t hidx(Halfedge_index h)
 // HyperIdealFunctional.java's defensive behaviour (lines 122-127 of the
 // Java original); this keeps the FD perturbation regime well-defined.
 
+/// Six per-face angle outputs computed from local DOFs (see
+/// `face_angles_from_local_dofs`).  Used by the block-FD Hessian.
 struct FaceAngleOutputs {
-    double beta1, beta2, beta3;       ///< interior angles at v₁,v₂,v₃
-    double alpha12, alpha23, alpha31; ///< dihedral angles at e₁₂,e₂₃,e₃₁
+    double beta1;   ///< Interior angle at v₁.
+    double beta2;   ///< Interior angle at v₂.
+    double beta3;   ///< Interior angle at v₃.
+    double alpha12; ///< Dihedral angle at edge e₁₂.
+    double alpha23; ///< Dihedral angle at edge e₂₃.
+    double alpha31; ///< Dihedral angle at edge e₃₁.
 };
 
+/// Pure-math 6→6 kernel: given the six local DOFs (b₁,b₂,b₃,a₁₂,a₂₃,a₃₁)
+/// of one face plus the per-vertex variability flags, return the six
+/// HyperIdeal angle outputs.  No mesh, no property maps — used by the
+/// per-face block-FD Hessian in `hyper_ideal_hessian.hpp`.
 inline FaceAngleOutputs face_angles_from_local_dofs(
     double b1, double b2, double b3,
     double a12, double a23, double a31,
@@ -196,14 +213,29 @@ inline FaceAngleOutputs face_angles_from_local_dofs(
 
 // ── Per-face angle kernel ─────────────────────────────────────────────────────
 
+/// Per-face angle bundle returned by `compute_face_angles()`.  Carries
+/// the six output angles plus the six input DOFs (so the energy and
+/// gradient kernels can reuse them without re-reading the mesh).
 struct FaceAngles {
-    double alpha12, alpha23, alpha31;  // dihedral angles at each edge
-    double beta1,   beta2,   beta3;    // interior angles at each vertex
-    double a12, a23, a31;              // edge DOF values (used in energy)
-    double b1,  b2,  b3;              // vertex DOF values
-    bool   v1b, v2b, v3b;             // whether each vertex is variable
+    double alpha12; ///< Dihedral angle at edge e₁₂.
+    double alpha23; ///< Dihedral angle at edge e₂₃.
+    double alpha31; ///< Dihedral angle at edge e₃₁.
+    double beta1;   ///< Interior angle at vertex v₁.
+    double beta2;   ///< Interior angle at vertex v₂.
+    double beta3;   ///< Interior angle at vertex v₃.
+    double a12;     ///< Edge DOF value at e₁₂.
+    double a23;     ///< Edge DOF value at e₂₃.
+    double a31;     ///< Edge DOF value at e₃₁.
+    double b1;      ///< Vertex DOF value at v₁.
+    double b2;      ///< Vertex DOF value at v₂.
+    double b3;      ///< Vertex DOF value at v₃.
+    bool   v1b;     ///< `true` iff vertex v₁ is variable (not pinned).
+    bool   v2b;     ///< `true` iff vertex v₂ is variable.
+    bool   v3b;     ///< `true` iff vertex v₃ is variable.
 };
 
+/// Compute the six per-face angles (+ remember the input DOFs) for face
+/// `f` of `mesh`, given the current DOF vector `x` and DOF-index maps.
 static FaceAngles compute_face_angles(
     const ConformalMesh&       mesh,
     Face_index                 f,
@@ -281,7 +313,7 @@ static FaceAngles compute_face_angles(
     return fa;
 }
 
-// Per-face energy contribution U(f) (before subtracting θ·a and Θ·b terms).
+/// Per-face energy contribution U(f) before subtracting the θ·a and Θ·b terms.
 static double face_energy(const FaceAngles& fa)
 {
     double aa = fa.a12*fa.alpha12 + fa.a23*fa.alpha23 + fa.a31*fa.alpha31;
@@ -310,6 +342,14 @@ static double face_energy(const FaceAngles& fa)
 
 // ── Full evaluation ───────────────────────────────────────────────────────────
 
+/// Evaluate the HyperIdeal functional at DOF vector `x`.  Returns the
+/// energy value and (optionally) the gradient in a `HyperIdealResult`.
+///
+/// \param mesh           Triangle mesh carrying the DOF-index property maps.
+/// \param x              Current DOF vector (length = `hyper_ideal_dimension(...)`).
+/// \param m              Property-map bundle from `setup_hyper_ideal_maps(...)`.
+/// \param need_energy    If `true`, fill `result.energy` (default: `true`).
+/// \param need_gradient  If `true`, fill `result.gradient` (default: `true`).
 inline HyperIdealResult evaluate_hyper_ideal(
     ConformalMesh&             mesh,
     const std::vector<double>& x,
@@ -393,10 +433,10 @@ inline HyperIdealResult evaluate_hyper_ideal(
     return res;
 }
 
-// ── Finite-difference gradient check ─────────────────────────────────────────
-//
-// Returns true if  |G[i] − fd[i]| / max(1, |G[i]|)  <  tol  for all DOFs.
-// eps = step size, tol = tolerance (same defaults as Java FunctionalTest).
+/// Finite-difference gradient check (central differences).
+///
+/// Returns `true` iff `|G[i] − fd[i]| / max(1, |G[i]|) < tol` for every
+/// DOF.  Defaults `eps = 1e-5`, `tol = 1e-4` match the Java `FunctionalTest`.
 inline bool gradient_check(
     ConformalMesh&             mesh,
     const std::vector<double>& x0,

code/include/hyper_ideal_geometry.hpp

@@ -22,9 +22,9 @@ namespace conformallab {
 
 // ── Length functions ─────────────────────────────────────────────────────────
 
-// ζ(x,y,z)  — interior angle in a hyperbolic triangle with edge lengths
-//             x, y, z, opposite to the side of length z.
-// Ports HyperIdealUtility.ζ(x, y, z).
+/// `ζ(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);
@@ -34,8 +34,8 @@ inline double zeta(double x, double y, double z)
     return std::acos(nbd);
 }
 
-// ζ₁₃(x,y,z) — third edge length in a right-angled hyperbolic hexagon.
-// Ports HyperIdealUtility.ζ_13(x, y, z).
+/// `ζ₁₃(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);
@@ -43,16 +43,16 @@ inline double zeta13(double x, double y, double z)
     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).
+/// `ζ₁₄(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).
+/// `ζ₁₅(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));
@@ -60,12 +60,11 @@ inline double zeta15(double x)
 
 // ── Effective edge length ─────────────────────────────────────────────────────
 
-// l_ij: effective hyperbolic length of edge ij.
-//   b_i, b_j – vertex log scale factors  (used only if vertex is hyper-ideal)
-//   a_ij     – edge intersection-angle variable
-//   vi_var   – true if vertex i is hyper-ideal (has a DOF b_i)
-//   vj_var   – true if vertex j is hyper-ideal
-// Ports HyperIdealFunctional.lij().
+/// `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);
@@ -76,8 +75,8 @@ inline double lij(double bi, double bj, double aij, bool vi_var, bool vj_var)
 
 // ── Auxiliary angle functions ─────────────────────────────────────────────────
 
-// σᵢ(aᵢⱼ, aₖᵢ, aⱼₖ, vj_var, vk_var) — intermediate half-length at vertex i.
-// Ports HyperIdealFunctional.σi().
+/// `σᵢ(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);
@@ -86,24 +85,24 @@ inline double sigma_i(double aij, double aki, double ajk, bool vj_var, bool vk_v
     return zeta15(ajk - aij - aki);
 }
 
-// σᵢⱼ(aᵢⱼ, bᵢ, bⱼ, vj_var) — intermediate half-length for edge ij from vertex i.
-// Ports HyperIdealFunctional.σij().
+/// `σᵢⱼ(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
-//   βi, βj, β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).
+/// `α_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,

code/include/hyper_ideal_hessian.hpp

@@ -54,13 +54,9 @@
 
 namespace conformallab {
 
-// ── Full finite-difference Hessian (baseline, Phase 4a) ──────────────────────
-//
-// Returns the n×n sparse Hessian, where n = hyper_ideal_dimension(mesh, m).
-// eps: finite-difference step size (default 1e-5 gives ~1e-10 relative error).
-//
-// Cost: n × full-gradient evaluations ≈ O(n·F).  Use for small meshes or
-// as a correctness reference for the block-FD variant below.
+/// Full finite-difference HyperIdeal Hessian (baseline, Phase 4a).
+/// Cost: `n` full-gradient evaluations ≈ `O(n·F)`.  Use for small
+/// meshes or as a correctness reference for the block-FD variant.
 inline Eigen::SparseMatrix<double> hyper_ideal_hessian(
     ConformalMesh&             mesh,
     const std::vector<double>& x,
@@ -96,10 +92,8 @@ inline Eigen::SparseMatrix<double> hyper_ideal_hessian(
     return H;
 }
 
-// ── Symmetrised Hessian (full-FD variant) ────────────────────────────────────
-//
-// The FD Hessian is symmetric in exact arithmetic; floating-point rounding
-// can introduce tiny asymmetries.  This helper returns (H + Hᵀ)/2.
+/// Symmetrised full-FD HyperIdeal Hessian: returns `(H + Hᵀ) / 2` to
+/// scrub the tiny asymmetries introduced by floating-point rounding.
 inline Eigen::SparseMatrix<double> hyper_ideal_hessian_sym(
     ConformalMesh&             mesh,
     const std::vector<double>& x,
@@ -137,6 +131,10 @@ inline Eigen::SparseMatrix<double> hyper_ideal_hessian_sym(
 //                                     vs full-FD ≈ 99,200    → ~33× speed-up.
 // On brezel.obj (F=13824, n≈14000): 165 888 face evaluations
 //                                     vs full-FD ≈ 193 M     → ~1166× speed-up.
+/// Per-face block-FD HyperIdeal Hessian (Phase 9b).  Uses the locality
+/// lemma `∂G_x/∂y = Σ_{f: x,y ∈ local(f)} ∂(β or α)/∂y` to perturb only
+/// the 6 face-local DOFs at a time, giving an `F·12` face-evaluation
+/// budget vs `n·F` for full-FD (~96× speed-up on brezel.obj).
 inline Eigen::SparseMatrix<double> hyper_ideal_hessian_block_fd(
     ConformalMesh&             mesh,
     const std::vector<double>& x,
@@ -213,11 +211,9 @@ inline Eigen::SparseMatrix<double> hyper_ideal_hessian_block_fd(
     return H;
 }
 
-// ── Symmetrised block-FD Hessian ─────────────────────────────────────────────
-//
-// The per-face block is symmetric to within FD rounding, but the
-// accumulation may amplify tiny asymmetries.  This helper returns
-// (H + Hᵀ)/2, identical in spirit to `hyper_ideal_hessian_sym`.
+/// Symmetrised block-FD HyperIdeal Hessian: returns `(H + Hᵀ) / 2` of
+/// `hyper_ideal_hessian_block_fd(...)` for downstream solvers that
+/// require strict symmetry.
 inline Eigen::SparseMatrix<double> hyper_ideal_hessian_block_fd_sym(
     ConformalMesh&             mesh,
     const std::vector<double>& x,

code/include/hyper_ideal_utility.hpp

@@ -12,9 +12,9 @@
 
 namespace conformallab {
 
-// Volume of a generalized hyperbolic tetrahedron with dihedral angles A..F.
-// Formula: Meyerhoff / Ushijima (Springer 2006).
-// Corresponds to Java HyperIdealUtility.calculateTetrahedronVolume().
+/// Volume of a generalized hyperbolic tetrahedron with dihedral
+/// angles `A,…,F` via the Meyerhoff / Ushijima 2006 formula.
+/// Same as Java `HyperIdealUtility.calculateTetrahedronVolume()`.
 inline double calculateTetrahedronVolume(double A, double B, double C,
                                          double D, double E, double F) {
     // PI from constants.hpp (conformallab::PI)
@@ -73,11 +73,9 @@ inline double calculateTetrahedronVolume(double A, double B, double C,
     return (U(z1) - U(z2)) / 2.0;
 }
 
-// Volume of a hyperideal tetrahedron with one ideal vertex (at gamma).
-// Dihedral angles at the ideal vertex: gamma1, gamma2, gamma3.
-// Dihedral angles at opposite edges: alpha23, alpha31, alpha12.
-// Formula: Kolpakov–Mednykh (arxiv math/0603097).
-// Corresponds to Java HyperIdealUtility.calculateTetrahedronVolumeWithIdealVertexAtGamma().
+/// Volume of a hyperideal tetrahedron with one ideal vertex at γ via
+/// the Kolpakov-Mednykh formula (arxiv math/0603097).  Same as Java
+/// `HyperIdealUtility.calculateTetrahedronVolumeWithIdealVertexAtGamma()`.
 inline double calculateTetrahedronVolumeWithIdealVertexAtGamma(
     double gamma1, double gamma2, double gamma3,
     double alpha23, double alpha31, double alpha12)

code/include/hyper_ideal_visualization_utility.hpp

@@ -32,9 +32,7 @@
 
 namespace conformallab {
 
-// ---------------------------------------------------------------------------
-// Circumcenter of three 2-D points
-// ---------------------------------------------------------------------------
+/// Circumcenter of three 2-D points (`a`, `b`, `c`) in the Euclidean plane.
 inline Eigen::Vector2d circumcenter2d(
     const Eigen::Vector2d& a,
     const Eigen::Vector2d& b,
@@ -54,10 +52,8 @@ inline Eigen::Vector2d circumcenter2d(
     return {ux, uy};
 }
 
-// ---------------------------------------------------------------------------
-// 4×4 Lorentz boost: maps the hyperboloid origin e₄=(0,0,0,1) to `center`.
-// `center` must lie on the hyperboloid: center[3]² - ‖center.head<3>()‖² = 1.
-// ---------------------------------------------------------------------------
+/// 4×4 Lorentz boost: maps the hyperboloid origin `e₄ = (0,0,0,1)` to
+/// `center`.  Precondition: `center` lies on the hyperboloid.
 inline Eigen::Matrix4d hyperboloidTranslation(const Eigen::Vector4d& center)
 {
     Eigen::Vector3d p = center.head<3>();
@@ -73,10 +69,8 @@ inline Eigen::Matrix4d hyperboloidTranslation(const Eigen::Vector4d& center)
     return T;
 }
 
-// ---------------------------------------------------------------------------
-// Project a hyperboloid point to the Poincaré disk (jReality convention:
-// add 1 to the w-coordinate, then dehomogenize the spatial part).
-// ---------------------------------------------------------------------------
+/// Project a hyperboloid point `x` onto the Poincaré disk (jReality
+/// convention: add 1 to the w-coordinate, then dehomogenise spatial part).
 inline Eigen::Vector2d toPoincareDisk(const Eigen::Vector4d& x)
 {
     double w = x(3) + 1.0;
@@ -97,6 +91,10 @@ inline Eigen::Vector2d toPoincareDisk(const Eigen::Vector4d& x)
 //
 // Port of HyperIdealVisualizationPlugin.getEuclideanCircleFromHyperbolic()
 // ---------------------------------------------------------------------------
+/// Convert a hyperbolic circle (`center` on the hyperboloid, hyperbolic
+/// `radius`) to the corresponding Euclidean circle in the Poincaré disk;
+/// returns `{cx, cy, r}`.  Port of `HyperIdealVisualizationPlugin
+/// .getEuclideanCircleFromHyperbolic()`.
 inline std::array<double,3> getEuclideanCircleFromHyperbolic(
     const Eigen::Vector4d& center, double radius)
 {

code/include/inversive_distance_functional.hpp

@@ -76,12 +76,17 @@ namespace conformallab {
 
 // ── Property-map type aliases ────────────────────────────────────────────────
 
+/// Property map vertex → `int`    for the Inversive-Distance functional.
 using IDVMapI = ConformalMesh::Property_map<Vertex_index, int>;
+/// Property map vertex → `double` for the Inversive-Distance functional.
 using IDVMapD = ConformalMesh::Property_map<Vertex_index, double>;
+/// Property map edge   → `double` for the Inversive-Distance functional.
 using IDEMapD = ConformalMesh::Property_map<Edge_index, double>;
 
 // ── Persistent map bundle ─────────────────────────────────────────────────────
 
+/// Bundle of the four property maps consumed by the Inversive-Distance
+/// circle-packing functional (Luo 2004 / Bowers-Stephenson 2004).
 struct InversiveDistanceMaps {
     IDVMapI v_idx;       ///< DOF index per vertex (−1 = pinned / u_v = 0)
     IDVMapD theta_v;     ///< target cone angle Θ_v   (default 2π)
@@ -225,12 +230,8 @@ inline double edge_length_squared(double u_i, double u_j, double I_ij) noexcept
 
 } // namespace id_detail
 
-// ── Gradient ──────────────────────────────────────────────────────────────────
-//
-// G_v = Θ_v − Σ_{f ∋ v} α_v(f)
-//
-// Halfedge convention (identical to euclidean_functional.hpp):
-//   h_alpha[h] = corner angle OPPOSITE to the edge of halfedge h in its face.
+/// Inversive-Distance gradient `G_v = Θ_v − Σ_faces α_v(face)`.  Same
+/// half-edge corner-angle storage convention as `euclidean_gradient`.
 inline std::vector<double> inversive_distance_gradient(
     const ConformalMesh&             mesh,
     const std::vector<double>&       x,
@@ -291,9 +292,8 @@ inline std::vector<double> inversive_distance_gradient(
     return G;
 }
 
-// ── Energy via Gauss-Legendre path integral ─────────────────────────────────
-//
-// E(u) = ∫₀¹ ⟨G(tu), u⟩ dt   (10-point GL, identical constants to euclidean_functional)
+/// Inversive-Distance energy `E(u) = ∫₀¹ ⟨G(t·u), u⟩ dt`, evaluated
+/// with 10-point Gauss-Legendre (constants shared with `euclidean_energy`).
 inline double inversive_distance_energy(
     const ConformalMesh&             mesh,
     const std::vector<double>&       x,
@@ -329,7 +329,7 @@ inline double inversive_distance_energy(
     return E;
 }
 
-// ── Finite-difference gradient check ─────────────────────────────────────────
+/// FD gradient check for the Inversive-Distance functional (central diff).
 inline bool gradient_check_inversive_distance(
     const ConformalMesh&             mesh,
     const std::vector<double>&       x,
@@ -358,7 +358,8 @@ inline bool gradient_check_inversive_distance(
     return true;
 }
 
-// ── Newton equilibrium check: gradient vanishes at converged u ──────────────
+/// Newton equilibrium check: returns `true` iff the gradient at `x`
+/// is below `tol` in infinity norm (Σ adj-face angles equal Θ_v).
 inline bool is_inversive_distance_equilibrium(
     const ConformalMesh&             mesh,
     const std::vector<double>&       x,

code/include/layout.hpp

@@ -75,28 +75,38 @@ namespace conformallab {
 // For hyperbolic holonomy the map is an orientation-preserving isometry of
 // the Poincaré disk (SU(1,1) element).
 // ─────────────────────────────────────────────────────────────────────────────
+/// Möbius transformation `T(z) = (a·z + b) / (c·z + d)` of the Riemann
+/// sphere; restricted to SU(1,1) for hyperbolic holonomy on the
+/// Poincaré disk.
 struct MobiusMap {
+    /// Complex scalar used for all entries.
     using C = std::complex<double>;
-    C a{1.0, 0.0};
-    C b{0.0, 0.0};
-    C c{0.0, 0.0};
-    C d{1.0, 0.0};
+    C a{1.0, 0.0}; ///< Top-left coefficient.
+    C b{0.0, 0.0}; ///< Top-right coefficient.
+    C c{0.0, 0.0}; ///< Bottom-left coefficient.
+    C d{1.0, 0.0}; ///< Bottom-right coefficient.
 
+    /// Apply the transformation to a complex point.
     C apply(C z) const { return (a * z + b) / (c * z + d); }
 
+    /// Apply the transformation to a 2-D real point (interpreted as `x + iy`).
     Eigen::Vector2d apply(const Eigen::Vector2d& p) const {
         C w = apply(C(p.x(), p.y()));
         return Eigen::Vector2d(w.real(), w.imag());
     }
 
+    /// Identity map.
     static MobiusMap identity() { return {C(1), C(0), C(0), C(1)}; }
 
+    /// Inverse map.
     MobiusMap inverse()               const { return {d, -b, -c, a}; }
+    /// Composition `(*this) ∘ T`, i.e. apply `T` first then `*this`.
     MobiusMap compose(const MobiusMap& T) const {
         return { a*T.a + b*T.c,  a*T.b + b*T.d,
                  c*T.a + d*T.c,  c*T.b + d*T.d };
     }
 
+    /// `true` iff the map is the identity up to tolerance `tol`.
     bool is_identity(double tol = 1e-9) const {
         if (std::abs(d) < 1e-14) return false;
         C a_ = a/d, b_ = b/d, c_ = c/d;
@@ -121,6 +131,9 @@ struct MobiusMap {
 
 // ── Result types ──────────────────────────────────────────────────────────────
 
+/// Result of a 2-D layout (`euclidean_layout`, `hyper_ideal_layout`):
+/// per-vertex UV coordinates plus a per-half-edge UV atlas for seamed
+/// textures.
 struct Layout2D {
     /// uv[v.idx()] — primary 2-D position (first / shallowest-BFS-depth visit).
     std::vector<Eigen::Vector2d> uv;
@@ -136,14 +149,15 @@ struct Layout2D {
     /// Size = mesh.number_of_halfedges().  Border halfedges = (0,0).
     std::vector<Eigen::Vector2d> halfedge_uv;
 
-    bool success  = false;
-    bool has_seam = false;   ///< true when a vertex was reached via two paths
+    bool success  = false; ///< `true` iff the BFS placed every vertex.
+    bool has_seam = false; ///< `true` when a vertex was reached via two paths.
 };
 
+/// Result of a 3-D layout (`spherical_layout`): per-vertex positions on S².
 struct Layout3D {
-    std::vector<Eigen::Vector3d> pos;
-    bool success  = false;
-    bool has_seam = false;
+    std::vector<Eigen::Vector3d> pos; ///< Per-vertex spherical positions.
+    bool success  = false; ///< `true` iff the BFS placed every vertex.
+    bool has_seam = false; ///< `true` when a vertex was reached via two paths.
 };
 
 /// Per-cut-edge holonomy.
@@ -156,9 +170,9 @@ struct Layout3D {
 ///   trilaterated virtual position obtained by continuing the unfolding across
 ///   the cut.
 struct HolonomyData {
-    std::vector<Eigen::Vector2d> translations;   ///< Euclidean / spherical
-    std::vector<MobiusMap>       mobius_maps;     ///< hyperbolic (Phase 7)
-    std::vector<std::size_t>     cut_edge_indices;
+    std::vector<Eigen::Vector2d> translations;     ///< Euclidean / spherical translation per cut edge.
+    std::vector<MobiusMap>       mobius_maps;      ///< Hyperbolic Möbius isometry per cut edge (Phase 7).
+    std::vector<std::size_t>     cut_edge_indices; ///< Index (in the cut-graph edge list) of each holonomy entry.
 };
 
 // ── Internal helpers ──────────────────────────────────────────────────────────
@@ -343,7 +357,8 @@ inline void center_poincare_disk_weighted(
 
 } // namespace detail
 
-// ── Vertex Voronoi area weights ───────────────────────────────────────────────
+/// Compute per-vertex area weights (sum of 1/3 of each adjacent triangle area).
+/// Used by area-weighted layout normalisation routines.
 inline std::vector<double> compute_vertex_area_weights(const ConformalMesh& mesh)
 {
     std::vector<double> w(mesh.number_of_vertices(), 0.0);
@@ -357,6 +372,8 @@ inline std::vector<double> compute_vertex_area_weights(const ConformalMesh& mesh
 
 // ── Layout normalisation ──────────────────────────────────────────────────────
 
+/// Euclidean canonical normalisation: translate centroid to the origin
+/// and rotate the principal axis of the UV cloud onto the x-axis.
 inline void normalise_euclidean(Layout2D& layout)
 {
     if (!layout.success || layout.uv.empty()) return;
@@ -388,6 +405,8 @@ inline void normalise_hyperbolic(Layout2D& layout, const ConformalMesh& mesh)
     // halfedge_uv follows the same Möbius map
     detail::center_poincare_disk(layout.halfedge_uv);
 }
+/// Hyperbolic canonical normalisation, mesh-free fallback: uniform
+/// (unweighted) iterative Möbius centring of the Poincaré disk.
 inline void normalise_hyperbolic(Layout2D& layout) // fallback without mesh
 {
     if (!layout.success || layout.uv.empty()) return;
@@ -395,6 +414,9 @@ inline void normalise_hyperbolic(Layout2D& layout) // fallback without mesh
     detail::center_poincare_disk(layout.halfedge_uv);
 }
 
+/// Spherical canonical normalisation: rotate the layout so that the
+/// per-vertex centroid (projected back to S²) coincides with the north
+/// pole (Rodrigues rotation).
 inline void normalise_spherical(Layout3D& layout)
 {
     if (!layout.success || layout.pos.empty()) return;
@@ -852,6 +874,8 @@ inline Layout2D hyper_ideal_layout(
 
 // ── Convenience: save layout as OFF ──────────────────────────────────────────
 
+/// Write a 2-D layout to disk in OFF format with z = 0.  Convenience
+/// helper for quickly inspecting the UV result in any OFF viewer.
 inline void save_layout_off(
     const std::string& path, ConformalMesh& mesh, const Layout2D& layout)
 {
@@ -865,6 +889,7 @@ inline void save_layout_off(
     }
 }
 
+/// Write a 3-D (spherical) layout to disk in OFF format.
 inline void save_layout_off(
     const std::string& path, ConformalMesh& mesh, const Layout3D& layout)
 {

code/include/matrix_utility.hpp

@@ -7,12 +7,10 @@
 
 namespace conformallab {
 
-// Find the 4×4 matrix R that maps source points to target points.
-// Each row of `from` / `to` is a homogeneous 4-vector (one point per row).
-// Post-condition: R * from.row(i).T == to.row(i).T for all i.
-//
-// Implementation: R = to^T * (from^T)^{-1}
-// Corresponds to Java MatrixUtility.makeMappingMatrix().
+/// Find the 4×4 matrix `R` that maps each row of `from` (homogeneous
+/// 4-vector) to the corresponding row of `to`: `R · fromᵀ = toᵀ`.
+/// Computed as `R = toᵀ · (fromᵀ)⁻¹`.  Same as Java
+/// `MatrixUtility.makeMappingMatrix()`.
 inline Eigen::Matrix4d makeMappingMatrix(const Eigen::Matrix4d& from,
                                          const Eigen::Matrix4d& to) {
     return to.transpose() * from.transpose().inverse();

code/include/mesh_builder.hpp

@@ -22,7 +22,7 @@ namespace conformallab {
 //   |    \
 //   v0 ─ v1
 //
-// Returns a mesh with 1 face, 3 vertices, 3 edges.
+/// Build a single right-angle triangle in the xy-plane (1 face, 3 vertices, 3 edges).
 // The triangle lies in the xy-plane with a right angle at v0.
 inline ConformalMesh make_triangle(
     double x0=0, double y0=0,
@@ -39,7 +39,7 @@ inline ConformalMesh make_triangle(
 
 // ── Regular tetrahedron ──────────────────────────────────────────────────────
 //
-// 4 vertices, 4 faces, 6 edges.
+/// Build a regular tetrahedron (4 vertices, 4 faces, 6 edges; sphere topology).
 // Euler characteristic: V - E + F = 4 - 6 + 4 = 2  (sphere topology).
 // Used to test closed-surface traversal.
 inline ConformalMesh make_tetrahedron()
@@ -67,8 +67,8 @@ inline ConformalMesh make_tetrahedron()
 //   |  \ |
 //   v0 ─ v1
 //
-// 4 vertices, 2 faces, 5 edges (1 interior edge v1–v2 shared by both faces).
-// Useful for testing edge-interior vs edge-boundary distinction.
+/// Build a two-triangle strip (4 vertices, 2 faces, 5 edges; 1 interior edge).
+/// Useful for testing interior- vs boundary-edge distinction.
 inline ConformalMesh make_quad_strip()
 {
     ConformalMesh mesh;
@@ -85,8 +85,8 @@ inline ConformalMesh make_quad_strip()
 
 // ── Regular flat polygon fan ─────────────────────────────────────────────────
 //
-// n triangles sharing a central vertex; forms a disk topology (boundary).
-// Used to verify valence-n vertex traversal.
+/// Build a regular flat polygon fan: `n` triangles sharing a central
+/// vertex, with rim vertices on the unit circle (disk topology).
 inline ConformalMesh make_fan(int n)
 {
     CGAL_precondition(n >= 3);
@@ -109,10 +109,9 @@ inline ConformalMesh make_fan(int n)
 
 // ── Spherical tetrahedron (vertices on the unit sphere) ───────────────────────
 //
-// The four vertices of a regular tetrahedron projected onto the unit sphere.
-// Starting from (±1,±1,±1), dividing by √3 gives unit-length positions.
-// All edge lengths equal arccos(−1/3) ≈ 1.9106 radians.
-// Used for SphericalFunctional tests (all four faces are valid spherical triangles).
+/// Build a regular tetrahedron with vertices on the unit sphere.
+/// All edge lengths equal `arccos(−1/3) ≈ 1.9106 rad`; used by the
+/// SphericalFunctional tests.
 inline ConformalMesh make_spherical_tetrahedron()
 {
     ConformalMesh mesh;
@@ -133,10 +132,9 @@ inline ConformalMesh make_spherical_tetrahedron()
 
 // ── Octahedron face triangle (vertices on the unit sphere) ────────────────────
 //
-// One face of a regular octahedron: the triangle (1,0,0)→(0,1,0)→(0,0,1).
-// All edge lengths equal arccos(0) = π/2.
-// The corner angles are all π/2 (right-angled spherical triangle).
-// base log-length: λ° = 2·log(sin(π/4)) = 2·log(1/√2) = −log(2) ≈ −0.6931.
+/// Build one face of a regular octahedron `(1,0,0)→(0,1,0)→(0,0,1)`:
+/// a right-angled spherical triangle with edge length `π/2` and base
+/// log-length `λ° = −log 2 ≈ −0.6931`.
 inline ConformalMesh make_octahedron_face()
 {
     ConformalMesh mesh;

code/include/mesh_io.hpp

@@ -28,26 +28,22 @@
 
 namespace conformallab {
 
-// ── Read ──────────────────────────────────────────────────────────────────────
-//
-// Reads a polygon mesh from file into `mesh` (clears any existing content).
-// Returns true on success, false on failure.
+/// Read a polygon mesh from `filename` into `mesh` (clears existing content).
+/// Returns `true` on success, `false` on failure.
 inline bool read_mesh(const std::string& filename, ConformalMesh& mesh)
 {
     mesh.clear();
     return CGAL::IO::read_polygon_mesh(filename, mesh);
 }
 
-// ── Write ─────────────────────────────────────────────────────────────────────
-//
-// Writes `mesh` to `filename`.  Returns true on success.
+/// Write `mesh` to `filename`.  Returns `true` on success.
 inline bool write_mesh(const std::string& filename, const ConformalMesh& mesh)
 {
     return CGAL::IO::write_polygon_mesh(filename, mesh);
 }
 
-// ── Convenience: throwing wrappers ────────────────────────────────────────────
-
+/// Throwing wrapper around `read_mesh`: returns the mesh by value
+/// or throws `std::runtime_error` on read failure.
 inline ConformalMesh load_mesh(const std::string& filename)
 {
     ConformalMesh mesh;
@@ -56,6 +52,8 @@ inline ConformalMesh load_mesh(const std::string& filename)
     return mesh;
 }
 
+/// Throwing wrapper around `write_mesh`; throws `std::runtime_error` on
+/// write failure.
 inline void save_mesh(const std::string& filename, const ConformalMesh& mesh)
 {
     if (!write_mesh(filename, mesh))

code/include/newton_solver.hpp

@@ -44,11 +44,12 @@ namespace conformallab {
 
 // ── Result ────────────────────────────────────────────────────────────────────
 
+/// Result of `newton_solve(...)` — converged DOF vector + diagnostics.
 struct NewtonResult {
-    std::vector<double> x;            ///< DOF vector at termination
-    int                 iterations;   ///< Newton steps taken
-    double              grad_inf_norm;///< max |G_i| at termination
-    bool                converged;    ///< true iff grad_inf_norm < tol
+    std::vector<double> x;            ///< DOF vector at termination.
+    int                 iterations;   ///< Newton steps taken.
+    double              grad_inf_norm;///< max |Gᵢ| at termination.
+    bool                converged;    ///< `true` iff `grad_inf_norm < tol`.
 };
 
 // ── Internal helpers ──────────────────────────────────────────────────────────
@@ -96,6 +97,10 @@ inline Eigen::VectorXd solve_with_fallback(
 //
 //   fallback_used  – if non-null, set to true iff SparseQR was invoked
 //   Returns Eigen::VectorXd::Zero(rhs.size()) if both solvers fail.
+/// Solve `A·x = rhs` with the same SimplicialLDLT → SparseQR fallback
+/// strategy used inside all three Newton solvers.  If `fallback_used`
+/// is non-null, it is set to `true` iff the SparseQR fallback ran.
+/// Returns `Eigen::VectorXd::Zero(rhs.size())` if both solvers fail.
 inline Eigen::VectorXd solve_linear_system(
     const Eigen::SparseMatrix<double>& A,
     const Eigen::VectorXd&             rhs,

code/include/p2_utility.hpp

@@ -13,9 +13,9 @@ namespace conformallab {
 
 // ── Point / line duality ──────────────────────────────────────────────────────
 
-// Intersection of two lines l1, l2  (or  line through two points p1, p2)
-// via the cross product.  Works for any P2 element.
-// Corresponds to Java P2.pointFromLines / P2.lineFromPoints.
+/// Cross-product point–line duality in P²: returns the intersection
+/// of two lines (or the line through two points).  Same as Java
+/// `P2.pointFromLines` / `P2.lineFromPoints`.
 inline Eigen::Vector3d pointFromLines(const Eigen::Vector3d& l1,
                                       const Eigen::Vector3d& l2) {
     return l1.cross(l2);
@@ -23,11 +23,9 @@ inline Eigen::Vector3d pointFromLines(const Eigen::Vector3d& l1,
 
 // ── Euclidean perpendicular bisector ─────────────────────────────────────────
 
-// Returns the homogeneous line coordinates (a, b, c) of the perpendicular
-// bisector of the segment [p, q] in the Euclidean plane.
-// Coordinates: ax + by + c = 0  (after dehomogenizing p and q).
-//
-// Corresponds to Java P2.perpendicularBisector(p, q, Pn.EUCLIDEAN).
+/// Homogeneous line coordinates `(a, b, c)` of the perpendicular
+/// bisector of `[p, q]` in the Euclidean plane (`ax + by + c = 0`).
+/// Same as Java `P2.perpendicularBisector(p, q, Pn.EUCLIDEAN)`.
 inline Eigen::Vector3d perpendicularBisectorEuclidean(const Eigen::Vector3d& p_h,
                                                       const Eigen::Vector3d& q_h) {
     // Dehomogenize
@@ -46,8 +44,7 @@ inline Eigen::Vector3d perpendicularBisectorEuclidean(const Eigen::Vector3d& p_h
     return {d(0), d(1), c};
 }
 
-// ── Euclidean distance between two P2 homogeneous points ─────────────────────
-
+/// Euclidean distance between two P² homogeneous points (dehomogenises both).
 inline double euclideanDistanceP2(const Eigen::Vector3d& p_h,
                                   const Eigen::Vector3d& q_h) {
     Eigen::Vector2d p = p_h.head<2>() / p_h(2);
@@ -57,11 +54,9 @@ inline double euclideanDistanceP2(const Eigen::Vector3d& p_h,
 
 // ── Direct Euclidean isometry from two point-frames ──────────────────────────
 
-// Build the 3×3 projective matrix that represents the coordinate frame
-// anchored at p0 with p1 defining the positive x-direction.
-// Euclidean case: columns are [dehom(p0), unit_dir(p0→p1), perp_dir].
-//
-// Template parameter S allows float / double / long double.
+/// Build the 3×3 projective frame matrix anchored at `p0` with `p1`
+/// defining the positive x-direction (Euclidean case).  Columns:
+/// `[dehom(p0), unit_dir(p0→p1), perp_dir]`.
 template <typename S>
 Eigen::Matrix<S, 3, 3> makeFrameMatrix(Eigen::Matrix<S, 3, 1> p0_h,
                                         Eigen::Matrix<S, 3, 1> p1_h) {
@@ -84,11 +79,9 @@ Eigen::Matrix<S, 3, 3> makeFrameMatrix(Eigen::Matrix<S, 3, 1> p0_h,
     return M;
 }
 
-// Find the 3×3 Euclidean isometry (as a projective matrix) that maps
-// the frame (s1, s2) to the frame (t1, t2).
-//
-// Corresponds to Java P2.makeDirectIsometryFromFrames(s1, s2, t1, t2, Pn.EUCLIDEAN)
-// and P2Big.makeDirectIsometryFromFrames(...) (the BigDecimal / high-precision variant).
+/// 3×3 Euclidean isometry (as a projective matrix) that maps the
+/// frame `(s1, s2)` to the frame `(t1, t2)`.  Same as Java
+/// `P2.makeDirectIsometryFromFrames(..., Pn.EUCLIDEAN)`.
 template <typename S>
 Eigen::Matrix<S, 3, 3> makeDirectIsometryFromFramesEuclidean(
     Eigen::Matrix<S, 3, 1> s1, Eigen::Matrix<S, 3, 1> s2,

code/include/period_matrix.hpp

@@ -50,6 +50,8 @@ namespace conformallab {
 // PeriodData
 // ─────────────────────────────────────────────────────────────────────────────
 
+/// Period-matrix data for a genus-g closed surface.  For genus 1 the
+/// conformal type is fully captured by `τ = ω₂ / ω₁ ∈ ℍ`.
 struct PeriodData {
     /// Lattice generators as complex numbers (one per cut edge).
     /// omega[i] = translations[i].x() + i·translations[i].y()
@@ -63,6 +65,7 @@ struct PeriodData {
     /// True if τ has been reduced to the standard fundamental domain.
     bool in_fundamental_domain = false;
 
+    /// Genus of the surface = `|omega| / 2`.
     int genus() const { return static_cast<int>(omega.size()) / 2; }
 };
 
@@ -74,6 +77,9 @@ struct PeriodData {
 //
 // Returns the reduced τ.  Throws if Im(τ) ≤ 0 (not in upper half-plane).
 // ─────────────────────────────────────────────────────────────────────────────
+/// Reduce `τ ∈ ℍ` to the standard SL(2,ℤ) fundamental domain
+/// `F = { τ ∈ ℍ : |τ| ≥ 1, −½ ≤ Re τ < ½ }` via the generators
+/// `S: τ↦−1/τ` and `T: τ↦τ+1`.  Throws if `Im τ ≤ 0`.
 inline std::complex<double> reduce_to_fundamental_domain(std::complex<double> tau)
 {
     if (tau.imag() <= 0.0) {
@@ -103,6 +109,8 @@ inline std::complex<double> reduce_to_fundamental_domain(std::complex<double> ta
 // ─────────────────────────────────────────────────────────────────────────────
 // is_in_fundamental_domain — check membership in F with tolerance tol.
 // ─────────────────────────────────────────────────────────────────────────────
+/// `true` iff `τ` lies inside the standard SL(2,ℤ) fundamental domain
+/// with tolerance `tol`.
 inline bool is_in_fundamental_domain(std::complex<double> tau, double tol = 1e-9)
 {
     if (tau.imag() <= 0.0) return false;
@@ -117,6 +125,9 @@ inline bool is_in_fundamental_domain(std::complex<double> tau, double tol = 1e-9
 // Computes the period data from the Euclidean holonomy translations.
 // For genus-1 surfaces, also reduces τ to the fundamental domain.
 // ─────────────────────────────────────────────────────────────────────────────
+/// Compute the period data from the Euclidean holonomy translations.
+/// For genus 1, also reduces `τ` to the SL(2,ℤ) fundamental domain
+/// when `reduce` is `true` (default).
 inline PeriodData compute_period_matrix(const HolonomyData& hol, bool reduce = true)
 {
     PeriodData pd;

code/include/projective_math.hpp

@@ -12,17 +12,15 @@
 
 namespace conformallab {
 
-// Divide a homogeneous vector by its last component.
-// Corresponds to Java Pn.dehomogenize().
+/// Dehomogenise: divide a homogeneous vector by its last component.
+/// Same as Java `Pn.dehomogenize()`.
 inline Eigen::VectorXd dehomogenize(const Eigen::VectorXd& p) {
     return p / p(p.size() - 1);
 }
 
-// Hyperbolic distance between two homogeneous vectors of the same dimension.
-// The last component is the "timelike" coordinate (jReality convention).
-// Inner product: <p,q> = -sum_i p_i*q_i + p_last * q_last
-// Distance: arcosh(<p̂, q̂>) where p̂ normalises to the hyperboloid.
-// Corresponds to Java Pn.distanceBetween(p, q, Pn.HYPERBOLIC).
+/// Hyperbolic distance `arcosh(⟨p̂, q̂⟩)` between two homogeneous
+/// vectors (last component = timelike coordinate; jReality convention).
+/// Same as Java `Pn.distanceBetween(p, q, Pn.HYPERBOLIC)`.
 inline double hyperbolicDistance(const Eigen::VectorXd& p,
                                  const Eigen::VectorXd& q) {
     int n = static_cast<int>(p.size());
@@ -34,10 +32,8 @@ inline double hyperbolicDistance(const Eigen::VectorXd& p,
     return std::acosh(std::max(1.0, inner));
 }
 
-// Check whether a homogeneous point p lies on the segment [s[0], s[1]].
-// Works for n-dimensional homogeneous coords; cross product uses the first
-// 3 spatial components after dehomogenization (matching jReality's Rn behaviour).
-// Corresponds to Java SurfaceCurveUtility.isOnSegment().
+/// `true` iff the homogeneous point `p_h` lies on the segment
+/// `[s0_h, s1_h]`.  Same as Java `SurfaceCurveUtility.isOnSegment()`.
 inline bool isOnSegment(const Eigen::VectorXd& p_h,
                         const Eigen::VectorXd& s0_h,
                         const Eigen::VectorXd& s1_h) {
@@ -63,10 +59,10 @@ inline bool isOnSegment(const Eigen::VectorXd& p_h,
     return true;
 }
 
-// Find the point on `target` that corresponds to `p` on `source`.
-// The parameter t is determined by hyperbolic distance ratios on `source`,
-// then applied as a linear interpolation on the dehomogenized `target`.
-// Corresponds to Java SurfaceCurveUtility.getPointOnCorrespondingSegment().
+/// Find the point on the target segment `(tgt0, tgt1)` corresponding
+/// to `p` on the source segment `(src0, src1)`, parametrised by
+/// hyperbolic distance ratios on the source.  Same as Java
+/// `SurfaceCurveUtility.getPointOnCorrespondingSegment()`.
 inline Eigen::VectorXd getPointOnCorrespondingSegment(
     const Eigen::VectorXd& p,
     const Eigen::VectorXd& src0,

code/include/serialization.hpp

@@ -43,7 +43,7 @@ namespace conformallab {
 // JSON
 // ════════════════════════════════════════════════════════════════════════════
 
-// Save solver result (+ optional 2D layout) to a JSON file.
+/// Save the Newton-solver result (+ optional 2-D layout) to a JSON file.
 inline void save_result_json(
     const std::string& path,
     const NewtonResult& res,
@@ -80,8 +80,8 @@ inline void save_result_json(
     ofs << std::setw(2) << j << "\n";
 }
 
-// Load DOF vector from a JSON result file.
-// Returns the DOF vector; fills out res fields if non-null.
+/// Load a DOF vector from a JSON result file written by
+/// `save_result_json`.  If `res` is non-null its fields are filled too.
 inline std::vector<double> load_result_json(
     const std::string& path,
     NewtonResult*       res     = nullptr,
@@ -181,7 +181,7 @@ inline std::vector<double> parse_doubles(const std::string& s)
 
 } // namespace detail_xml
 
-// Save solver result (+ optional layout) to an XML file.
+/// Save the Newton-solver result (+ optional layout) to an XML file.
 inline void save_result_xml(
     const std::string& path,
     const NewtonResult& res,
@@ -229,8 +229,9 @@ inline void save_result_xml(
     ofs << "</ConformalResult>\n";
 }
 
-// Load solver result from an XML file written by save_result_xml.
-// Returns the DOF vector; fills res/geom/layout2d if non-null.
+/// Load a DOF vector from an XML result file written by
+/// `save_result_xml`.  If `res`, `geom`, `layout2d` are non-null they
+/// are filled as well.
 inline std::vector<double> load_result_xml(
     const std::string& path,
     NewtonResult*       res      = nullptr,

code/include/spherical_functional.hpp

@@ -40,19 +40,24 @@ namespace conformallab {
 
 // ── Property-map type aliases ─────────────────────────────────────────────────
 
+/// Property map vertex → `double` for the Spherical functional.
 using SpherVMapD = ConformalMesh::Property_map<Vertex_index, double>;
+/// Property map vertex → `int`    for the Spherical functional.
 using SpherVMapI = ConformalMesh::Property_map<Vertex_index, int>;
+/// Property map edge   → `double` for the Spherical functional.
 using SpherEMapD = ConformalMesh::Property_map<Edge_index,   double>;
+/// Property map edge   → `int`    for the Spherical functional.
 using SpherEMapI = ConformalMesh::Property_map<Edge_index,   int>;
 
 // ── Persistent map bundle ─────────────────────────────────────────────────────
 
+/// Bundle of the five property maps consumed by the Spherical functional.
 struct SphericalMaps {
-    SpherVMapI v_idx;    // DOF index per vertex   (-1 = pinned / u_v = 0)
-    SpherEMapI e_idx;    // DOF index per edge      (-1 = no edge DOF)
-    SpherVMapD theta_v;  // target cone angle Θ_v  (default 2π)
-    SpherEMapD theta_e;  // target edge angle θ_e  (default π)
-    SpherEMapD lambda0;  // base log-length λ°_e   (default 0.0)
+    SpherVMapI v_idx;    ///< DOF index per vertex (−1 = pinned / u_v = 0).
+    SpherEMapI e_idx;    ///< DOF index per edge   (−1 = no edge DOF).
+    SpherVMapD theta_v;  ///< Target cone angle Θᵥ (default 2π).
+    SpherEMapD theta_e;  ///< Target edge angle θₑ (default π).
+    SpherEMapD lambda0;  ///< Base log-length λ⁰ₑ (default 0).
 };
 
 // Defaults: theta_v = 2π, theta_e = π, lambda0 = 0.
@@ -133,18 +138,21 @@ inline void compute_lambda0_from_mesh(ConformalMesh& mesh, SphericalMaps& m)
 
 // ── Evaluation result ─────────────────────────────────────────────────────────
 
+/// Output of `evaluate_spherical()` — energy plus optional gradient.
 struct SphericalResult {
-    double              energy   = 0.0;
-    std::vector<double> gradient;
+    double              energy   = 0.0; ///< Functional value at input DOFs.
+    std::vector<double> gradient;       ///< Gradient ∇E (empty if not requested).
 };
 
 // ── Internal helpers ──────────────────────────────────────────────────────────
 
+/// Read DOF value from `x` for index `idx`; return 0 if pinned (idx < 0).
 static inline double spher_dof_val(int idx, const std::vector<double>& x)
 {
     return idx >= 0 ? x[static_cast<std::size_t>(idx)] : 0.0;
 }
 
+/// Convert a CGAL half-edge index to a plain `std::size_t` for vector indexing.
 static inline std::size_t spher_hidx(Halfedge_index h)
 {
     return static_cast<std::size_t>(static_cast<std::uint32_t>(h));
@@ -152,14 +160,13 @@ static inline std::size_t spher_hidx(Halfedge_index h)
 
 // ── Gradient only (no energy) ─────────────────────────────────────────────────
 
-// Compute gradient G(x).
-// G_v = Θ_v − Σ_faces α_v(face)
-// G_e = α_opp(face+) + α_opp(face−) − θ_e
-//
-// The corner angle α_v is stored on halfedges using the convention:
-//   h_alpha[h] = corner angle at source(prev(h)) = corner angle at the vertex
-//                ACROSS FROM the edge of halfedge h in its face.
-// This convention makes both the vertex and edge gradient accumulators natural.
+/// Compute the Spherical-functional gradient G(x):
+///   * `G_v = Θ_v − Σ_faces α_v(face)`
+///   * `G_e = α_opp(face⁺) + α_opp(face⁻) − θ_e`
+///
+/// The corner angle α_v is stored on half-edges via the convention
+/// `h_alpha[h] = corner angle at the vertex ACROSS FROM the edge of h
+/// in its face`, which makes both gradient accumulators natural.
 inline std::vector<double> spherical_gradient(
     ConformalMesh&             mesh,
     const std::vector<double>& x,
@@ -274,6 +281,9 @@ inline std::vector<double> spherical_gradient(
 //
 // 10-point GL nodes and weights on [0, 1] (transformed from [-1, 1]):
 //   t_k = (1 + s_k) / 2,  w_k = w_GL_k / 2
+/// Spherical energy `E(x) = ∫₀¹ ⟨G(t·x), x⟩ dt`, evaluated with
+/// 10-point Gauss-Legendre quadrature.  This is the correct potential
+/// for any conservative `G = ∇E`; error ≈ O(h²⁰) for smooth G.
 inline double spherical_energy(
     ConformalMesh&             mesh,
     const std::vector<double>& x,
@@ -318,6 +328,8 @@ inline double spherical_energy(
 
 // ── Full evaluation (energy + gradient) ──────────────────────────────────────
 
+/// Evaluate the Spherical functional at DOFs `x`.  Returns energy and
+/// gradient (toggle via `need_energy` / `need_gradient`).
 inline SphericalResult evaluate_spherical(
     ConformalMesh&             mesh,
     const std::vector<double>& x,
@@ -333,10 +345,8 @@ inline SphericalResult evaluate_spherical(
     return res;
 }
 
-// ── Finite-difference gradient check ─────────────────────────────────────────
-//
-// Tests |G[i] − fd[i]| / max(1, |G[i]|) < tol for all DOFs.
-// Same defaults as the hyper-ideal gradient check (Java FunctionalTest).
+/// Finite-difference gradient check for the Spherical functional
+/// (central differences).  Same defaults as the Java `FunctionalTest`.
 inline bool gradient_check_spherical(
     ConformalMesh&             mesh,
     const std::vector<double>& x0,
@@ -390,6 +400,8 @@ inline bool gradient_check_spherical(
 //
 // Returns 0.0 if the zero cannot be bracketed (already at gauge maximum,
 // or open surface — no shift needed).
+/// Find the global-scale gauge shift `t*` for the closed-spherical case
+/// (see comment block above for the maths).  Apply via `apply_spherical_gauge`.
 inline double spherical_gauge_shift(
     ConformalMesh&             mesh,
     const std::vector<double>& x,
@@ -470,7 +482,8 @@ inline double spherical_gauge_shift(
     return t;
 }
 
-// Apply the gauge shift in-place: x_v ← x_v + t* for all variable vertices.
+/// Apply the spherical gauge shift in-place: `x_v ← x_v + t*` for every
+/// variable vertex, where `t* = spherical_gauge_shift(mesh, x, m, ...)`.
 inline void apply_spherical_gauge(
     ConformalMesh&       mesh,
     std::vector<double>& x,

code/include/spherical_geometry.hpp

@@ -23,8 +23,8 @@ constexpr double PI_SPHER = PI;
 
 // ── Effective spherical arc length ────────────────────────────────────────────
 
-// l(λ) = 2·asin(min(exp(λ/2), 1)).
-// Clamps exp(λ/2) to [0, 1] so the arcsin stays in domain.
+/// 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);
@@ -35,29 +35,18 @@ inline double spherical_l(double lambda)
 
 // ── Interior angles of a spherical triangle ──────────────────────────────────
 
+/// Interior angles of a spherical triangle, plus a `valid` flag.
 struct SphericalFaceAngles {
-    double alpha1, alpha2, alpha3;  // corner angles at v1, v2, v3
-    bool   valid;                   // false when the three lengths fail the
-                                    // spherical triangle inequality
+    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 corner angles from spherical arc lengths using the half-angle formula.
-//
-// Convention (matching the halfedge cycle h0→v1→v2, h1→v2→v3, h2→v3→v1):
-//   l12 – arc length of edge opposite v3 (edge e12)
-//   l23 – arc length of edge opposite v1 (edge e23)
-//   l31 – arc length of edge opposite v2 (edge e31)
-//
-// Half-angle formula (spherical law of cosines):
-//   α_k = 2·atan2(sqrt(sin(s-a)·sin(s-b)), sqrt(sin(s)·sin(s-c)))
-// where a,b are the two edges ADJACENT to vertex k, c is the opposite edge.
-//
-// Equivalently (in terms of s-deficiencies):
-//   α1 = 2·atan2( sqrt(sin(s12)·sin(s31)), sqrt(sin(s)·sin(s23)) )
-//   α2 = 2·atan2( sqrt(sin(s12)·sin(s23)), sqrt(sin(s)·sin(s31)) )
-//   α3 = 2·atan2( sqrt(sin(s23)·sin(s31)), sqrt(sin(s)·sin(s12)) )
-//
-// where s = (l12+l23+l31)/2 and s_ij = s - l_ij.
+/// 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;

code/include/spherical_hessian.hpp

@@ -49,8 +49,18 @@ namespace conformallab {
 //   h2: edge v3-v1 → opposite v2 → w = cot(β2), β2=(π-α3-α1+α2)/2
 //
 // Returns valid=false if any β_k is out of range (degenerate face).
-struct SpherCotWeights { double w12, w23, w31; bool valid; };
+/// Three spherical "cotangent" weights for the three edges of a face,
+/// derived from the per-vertex interior angles `α₁, α₂, α₃` via
+/// `w_ij = cot(β_k)` with `β_k = (π − α_i − α_j + α_k) / 2`.
+struct SpherCotWeights {
+    double w12; ///< Weight for edge v₁-v₂ (opposite vertex v₃).
+    double w23; ///< Weight for edge v₂-v₃ (opposite vertex v₁).
+    double w31; ///< Weight for edge v₃-v₁ (opposite vertex v₂).
+    bool   valid; ///< `false` when any β_k is out of `(0, π/2]` (degenerate face).
+};
 
+/// Compute the three spherical cot weights from the three interior
+/// angles `(α₁, α₂, α₃)` of a spherical triangle.  See `SpherCotWeights`.
 inline SpherCotWeights spherical_cot_weights(double alpha1, double alpha2, double alpha3)
 {
     // β for each edge:
@@ -79,25 +89,10 @@ inline SpherCotWeights spherical_cot_weights(double alpha1, double alpha2, doubl
     return {1.0 / tb3, 1.0 / tb1, 1.0 / tb2, true};
 }
 
-// ── Analytical Hessian ────────────────────────────────────────────────────────
-//
-// Returns the n×n sparse Hessian matrix H where n = spherical_dimension(mesh, m).
-// x – current DOF vector.
-//
-// Derivation: G_v = θ_v − Σ_f α_v^f  →  H[i,j] = −Σ_f ∂α_i^f/∂u_j
-//
-// For a face (v1,v2,v3) with arc-lengths l12,l23,l31 and angles α1,α2,α3,
-// differentiating the spherical law of cosines
-//   cos(l_opp) = cos(l_a)cos(l_b) + sin(l_a)sin(l_b)cos(α)
-// gives:
-//   ∂α1/∂l12 = [cot(l12)cos(α1) − cot(l31)] / sin(α1)     (adjacent side)
-//   ∂α1/∂l31 = [cot(l31)cos(α1) − cot(l12)] / sin(α1)     (adjacent side)
-//   ∂α1/∂l23 = sin(l23) / [sin(l12)sin(l31)sin(α1)]        (opposite side)
-//
-// Chain rule with ∂l_ij/∂u_k = tan(l_ij/2)  (from l = 2·asin(exp(λ/2))):
-//   ∂α1/∂u1 = ∂α1/∂l12·t12 + ∂α1/∂l31·t31
-//   ∂α1/∂u2 = ∂α1/∂l12·t12 + ∂α1/∂l23·t23
-//   ∂α1/∂u3 = ∂α1/∂l23·t23 + ∂α1/∂l31·t31
+/// Analytical Spherical Hessian via `∂α/∂u` from the spherical law of
+/// cosines + chain rule `∂l/∂u = tan(l/2)`; returns an n×n sparse
+/// matrix with `n = spherical_dimension(mesh, m)`.  See block comment
+/// inside the body for the per-face derivation.
 inline Eigen::SparseMatrix<double> spherical_hessian(
     ConformalMesh&             mesh,
     const std::vector<double>& x,
@@ -214,10 +209,8 @@ inline Eigen::SparseMatrix<double> spherical_hessian(
     return H;
 }
 
-// ── Finite-difference Hessian check ──────────────────────────────────────────
-//
-// Compares the analytical Hessian column-by-column against
-//   H_fd[:, j] = (G(x + ε·eⱼ) − G(x − ε·eⱼ)) / (2ε).
+/// FD Hessian check for the Spherical functional.  Compares analytic
+/// `H` column-by-column to `(G(x+εeⱼ) − G(x−εeⱼ)) / (2ε)`.
 inline bool hessian_check_spherical(
     ConformalMesh&             mesh,
     const std::vector<double>& x0,

code/include/viewer_utils.h

@@ -5,7 +5,9 @@
 
 namespace viewer_utils {
 
-// Deklaration (Implementation in viewer.cpp)
+/// Open an interactive libigl OpenGL viewer window showing the mesh
+/// `(V, F)`.  Built only when `WITH_VIEWER=ON`; declaration here, body
+/// in `viewer.cpp`.
 void simple_visualize(Eigen::MatrixXd& V, Eigen::MatrixXi& F);
 
 }