Merge pull request 'docs: drive Doxygen coverage to 100 percent + auto-generate headers.md' (#17) from docs/doxygen-coverage-100 into main

2026-05-26 09:14:15 +00:00
parent 8869ead3c9 8e9ec9eccf
commit 63f2b9799d
47 changed files with 1354 additions and 482 deletions
--- a/.gitea/workflows/doxygen-pages.yml
+++ b/.gitea/workflows/doxygen-pages.yml
@@ -47,6 +47,28 @@ jobs:
            head -30 doc/doxygen/doxygen-warnings.log
          fi

+      - 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: |
+          python3 scripts/gen-headers-md.py
+          # If the auto-generated headers.md drifted from main, note it.
+          # This job runs on every main push so a drift only persists
+          # for the duration of one push — the next push that lands
+          # will fold the new headers.md back into main (via the
+          # codeberg pages branch).  For deterministic regeneration
+          # within main itself, run `bash scripts/regen-docs.sh`
+          # locally before pushing.
+          if ! git diff --quiet -- doc/api/headers.md; then
+            echo "::warning::doc/api/headers.md drifted — run scripts/regen-docs.sh locally and commit before next push"
+            git --no-pager diff -- doc/api/headers.md | head -30
+          fi
+
      - name: Publish HTML to codeberg pages branch
        env:
          CODEBERG_TOKEN: ${{ secrets.CODEBERG_TOKEN }}
--- a/42
+++ b/42
@@ -31,7 +31,15 @@ EXCLUDE_PATTERNS       = */build*/* \
                         */deps/* \
                         */.git/* \
                         */test-reports/* \
-                         */* 2.hpp
+                         *\ 2.hpp \
+                         *\ 2.h
+# Research-quality LaTeX notes use raw \sinh / \cosh / \frac / \beta /
+# \cdot / \partial / \zeta macros which are valid LaTeX but unknown to
+# Doxygen.  These files are intended to be read as PDF or in a LaTeX-
+# aware markdown viewer, not as Doxygen pages.  Excluding them removes
+# ~500 spurious "unknown command" warnings while keeping the .md files
+# discoverable on GitHub.
+EXCLUDE                = doc/math/hyperideal-hessian-derivation.md
 EXCLUDE_SYMBOLS        = Eigen::* boost::* std::*

 # Markdown filter: rewrites repo-relative links like [x](doc/api/tests.md)
@@ -70,7 +78,7 @@ EXTENSION_MAPPING      = h=C++ hpp=C++
 # ── Warnings ─────────────────────────────────────────────────────────────────
 QUIET                  = NO
 WARNINGS               = YES
-WARN_IF_UNDOCUMENTED   = NO
+WARN_IF_UNDOCUMENTED   = YES
 WARN_IF_DOC_ERROR      = YES
 WARN_IF_INCOMPLETE_DOC = YES
 WARN_NO_PARAMDOC       = NO
@@ -80,13 +88,24 @@ WARN_LOGFILE           = doc/doxygen/doxygen-warnings.log

 # ── HTML output ──────────────────────────────────────────────────────────────
 GENERATE_HTML          = YES
+
+# MathJax — render LaTeX math in markdown ($...$ and $$...$$) and in
+# code-comment `\f$ ... \f$` blocks via MathJax in the generated HTML.
+# Required for the conformal-mapping math notation (\Theta, \omega, \tau,
+# \mathbb{H}, …) in doc/architecture/overall_pipeline.md and the
+# header docstrings.
+USE_MATHJAX            = YES
+MATHJAX_VERSION        = MathJax_3
+MATHJAX_FORMAT         = HTML-CSS
+MATHJAX_RELPATH        = https://cdn.jsdelivr.net/npm/mathjax@3/es5/
 HTML_OUTPUT            = html
 HTML_FILE_EXTENSION    = .html
 HTML_COLORSTYLE        = LIGHT
 HTML_COLORSTYLE_HUE    = 220
 HTML_COLORSTYLE_SAT    = 100
 HTML_COLORSTYLE_GAMMA  = 80
-HTML_TIMESTAMP         = NO
+# HTML_TIMESTAMP was removed in Doxygen 1.10; use TIMESTAMP=NO instead.
+TIMESTAMP              = NO
 HTML_DYNAMIC_SECTIONS  = YES
 GENERATE_TREEVIEW      = YES
 DISABLE_INDEX          = NO
@@ -100,7 +119,9 @@ SERVER_BASED_SEARCH    = NO
 GENERATE_LATEX         = NO
 GENERATE_RTF           = NO
 GENERATE_MAN           = NO
-GENERATE_XML           = NO
+GENERATE_XML           = YES
+XML_OUTPUT             = xml
+XML_PROGRAMLISTING     = NO
 GENERATE_DOCBOOK       = NO
 GENERATE_AUTOGEN_DEF   = NO
 GENERATE_PERLMOD       = NO
@@ -130,3 +151,16 @@ ALIASES               += "concept{1}=\xrefitem concept \"Concept\" \"Concepts\"
 ALIASES               += "models{1}=\xrefitem models \"Models\" \"Models\" \1"
 ALIASES               += "cgalRequires{1}=\par Requirements: \n\1"
 ALIASES               += "cgalParam{2}=\param \1 \2"
+# CGAL named-parameter block aliases — replicates the upstream
+# ${CGAL}/Documentation/doc/Documentation/Doxyfile_common conventions
+# so that \cgalParamNBegin{name} … \cgalParamNEnd blocks render as
+# nested HTML lists in our Doxygen output.
+ALIASES               += "cgalNamedParamsBegin=<dl class=\"params\"><dt>Optional named parameters</dt><dd><table class=\"params\">"
+ALIASES               += "cgalNamedParamsEnd=</table></dd></dl>"
+ALIASES               += "cgalParamNBegin{1}=<tr><td class=\"paramname\"><code>\1</code></td><td>"
+ALIASES               += "cgalParamNEnd=</td></tr>"
+ALIASES               += "cgalParamDescription{1}=<b>Description:</b> \1<br/>"
+ALIASES               += "cgalParamType{1}=<b>Type:</b> \1<br/>"
+ALIASES               += "cgalParamDefault{1}=<b>Default:</b> \1<br/>"
+ALIASES               += "cgalParamPrecondition{1}=<b>Precondition:</b> \1<br/>"
+ALIASES               += "cgalParamExtra{1}=<i>\1</i><br/>"
--- a/code/include/CGAL/Conformal_map/doxygen_groups.h
+++ b/code/include/CGAL/Conformal_map/doxygen_groups.h
@@ -0,0 +1,44 @@
+// Copyright (c) 2024-2026 Tarik Moussa.
+// SPDX-License-Identifier: MIT
+//
+// This header contains only Doxygen \defgroup commands.  It is included
+// nowhere in the build; its sole purpose is to register the package's
+// Doxygen group hierarchy so that `@ingroup Pkg...` references in the
+// other public headers resolve cleanly.
+
+#ifndef CGAL_CONFORMAL_MAP_DOXYGEN_GROUPS_H
+#define CGAL_CONFORMAL_MAP_DOXYGEN_GROUPS_H
+
+/*!
+\defgroup PkgConformalMap CGAL Discrete Conformal Map package
+\brief Discrete conformal maps on triangulated surfaces — five DCE models
+(Euclidean, Spherical, Hyper-Ideal, Circle-Packing Euclidean,
+Inversive-Distance).
+
+This package provides the C++ implementation of the variational discrete
+conformal equivalence solvers from Springborn 2020, Bobenko/Pinkall/Springborn
+2010, and Luo 2004, together with a CGAL-style named-parameter API.
+
+See `doc/api/cgal-package.md` for the full design rationale.
+*/
+
+/*!
+\defgroup PkgConformalMapRef Reference manual
+\ingroup PkgConformalMap
+\brief Public C++ API: entry functions, traits, layout helpers.
+*/
+
+/*!
+\defgroup PkgConformalMapConcepts Concepts
+\ingroup PkgConformalMap
+\brief C++ concepts and traits classes consumed by the entry functions.
+*/
+
+/*!
+\defgroup PkgConformalMapNamedParameters Named function parameters
+\ingroup PkgConformalMap
+\brief Package-specific named-parameter helpers in `CGAL::parameters::*`,
+plus the pipe-operator chaining convention.
+*/
+
+#endif  // CGAL_CONFORMAL_MAP_DOXYGEN_GROUPS_H
--- a/code/include/CGAL/Conformal_map/doxygen_namespaces.h
+++ b/code/include/CGAL/Conformal_map/doxygen_namespaces.h
@@ -0,0 +1,88 @@
+// Copyright (c) 2024-2026 Tarik Moussa.
+// SPDX-License-Identifier: MIT
+//
+// Doxygen namespace documentation only — no declarations.  Centralised
+// here so that each namespace gets a single, consistent description in
+// the generated HTML, regardless of which header is parsed first.
+
+#ifndef CGAL_CONFORMAL_MAP_DOXYGEN_NAMESPACES_H
+#define CGAL_CONFORMAL_MAP_DOXYGEN_NAMESPACES_H
+
+/*!
+\namespace CGAL
+\brief Root namespace of the CGAL library; conformallab++ adds its
+public entry points (`discrete_conformal_map_*`, `Conformal_map_traits`,
+…) directly into this namespace, matching CGAL package conventions.
+*/
+
+/*!
+\namespace CGAL::Conformal_map
+\brief Implementation-detail namespace for the Discrete Conformal Map
+package.  Users normally do not need to enter this namespace; all
+public entry points are re-exported into `CGAL::`.
+*/
+
+/*!
+\namespace CGAL::Conformal_map::internal_np
+\brief Tag types backing the package-local named-function parameters
+(`vertex_curvature_map_t`, `gradient_tolerance_t`, `output_uv_map_t`, …).
+Users invoke them via the helpers in `CGAL::parameters::*`.
+*/
+
+/*!
+\namespace CGAL::parameters
+\brief CGAL named-function-parameter helpers — both upstream CGAL's and
+the conformallab++ package extensions (`vertex_curvature_map(...)`,
+`gradient_tolerance(...)`, `output_uv_map(...)`, `normalise_layout(...)`).
+Also home of the pipe-operator chaining convention; see
+`doc/tutorials/add-output-uv-map.md` §3.4.
+*/
+
+/*!
+\namespace conformallab
+\brief Core math/algorithm namespace of conformallab++.  Holds the five
+DCE functionals (Euclidean / Spherical / HyperIdeal / CP-Euclidean /
+Inversive-Distance), the Newton solver, layout helpers, mesh-property
+typedefs, and serialisation utilities.  Lives under
+`code/include/*.hpp` and is consumed both by the standalone CLI and
+by the thin CGAL wrappers under `CGAL::`.
+*/
+
+/*!
+\namespace conformallab::detail
+\brief Implementation-private helpers for the `conformallab` namespace.
+Not part of the stable public API.
+*/
+
+/*!
+\namespace conformallab::cp_detail
+\brief Implementation-private helpers for the Circle-Packing Euclidean
+functional (see `cp_euclidean_functional.hpp`).
+*/
+
+/*!
+\namespace conformallab::id_detail
+\brief Implementation-private helpers for the Inversive-Distance
+functional (see `inversive_distance_functional.hpp`).
+*/
+
+/*!
+\namespace conformallab::detail_xml
+\brief Implementation-private XML helpers for the (de)serialisation
+layer (see `serialization.hpp`).
+*/
+
+/*!
+\namespace mesh_utils
+\brief Small, opinion-free mesh utilities (loaders, validators,
+property-map registration) used by both the standalone tools and the
+CGAL wrappers.
+*/
+
+/*!
+\namespace viewer_utils
+\brief libigl-based interactive viewer helpers; built only when
+`WITH_VIEWER=ON`.  Not part of the headless / CGAL public surface.
+*/
+
+#endif  // CGAL_CONFORMAL_MAP_DOXYGEN_NAMESPACES_H
--- a/code/include/CGAL/Conformal_map/internal/parameters.h
+++ b/code/include/CGAL/Conformal_map/internal/parameters.h
@@ -198,11 +198,8 @@ inline auto normalise_layout(bool flag)
 // operator on the same type.
 // ════════════════════════════════════════════════════════════════════════════

-/*!
-\addtogroup PkgConformalMapNamedParameters
-\{
-*/
-
+// Close the PkgConformalMapNamedParameters group block that was opened
+// above the helper functions (see \addtogroup at the top of this section).
 /// \}

 }  // namespace parameters
--- a/code/include/CGAL/Conformal_map_traits.h
+++ b/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;
@@ -110,20 +113,33 @@ selected automatically when `TriangleMesh = CGAL::Surface_mesh<...>`.
 template <typename K>
 struct Default_conformal_map_traits<CGAL::Surface_mesh<typename K::Point_3>, K>
 {
+    /// The CGAL kernel parameter; defaults to `Simple_cartesian<double>`.
    using Kernel        = K;
+    /// Field type used for all scalar conformal-map data (lengths, λ, Θ, …).
    using FT            = typename K::FT;
+    /// 3-D point type used for vertex coordinates.
    using Point_3       = typename K::Point_3;
+    /// The triangle-mesh type this specialisation targets.
    using Triangle_mesh = CGAL::Surface_mesh<Point_3>;

+    /// Boost-graph vertex descriptor for `Triangle_mesh`.
    using Vertex_descriptor   = typename boost::graph_traits<Triangle_mesh>::vertex_descriptor;
+    /// Boost-graph half-edge descriptor for `Triangle_mesh`.
    using Halfedge_descriptor = typename boost::graph_traits<Triangle_mesh>::halfedge_descriptor;
+    /// Boost-graph edge descriptor for `Triangle_mesh`.
    using Edge_descriptor     = typename boost::graph_traits<Triangle_mesh>::edge_descriptor;
+    /// Boost-graph face descriptor for `Triangle_mesh`.
    using Face_descriptor     = typename boost::graph_traits<Triangle_mesh>::face_descriptor;

    // Property-map types — match the names used by setup_euclidean_maps().
+
+    /// Property map vertex → `Point_3` (the mesh's geometric embedding).
    using Vertex_point_map      = typename Triangle_mesh::template Property_map<Vertex_descriptor, Point_3>;
+    /// Property map vertex → target cone angle Θᵥ in radians (legacy name `ev:theta`).
    using Theta_pmap            = typename Triangle_mesh::template Property_map<Vertex_descriptor, FT>;
+    /// Property map vertex → contiguous integer index (legacy name `ev:idx`).
    using Vertex_index_pmap     = typename Triangle_mesh::template Property_map<Vertex_descriptor, int>;
+    /// Property map edge → log of original edge length λ⁰ (legacy name `ee:lam0`).
    using Lambda0_pmap          = typename Triangle_mesh::template Property_map<Edge_descriptor, FT>;

    // ─── Property-map accessors ───────────────────────────────────────────
@@ -133,10 +149,12 @@ struct Default_conformal_map_traits<CGAL::Surface_mesh<typename K::Point_3>, K>
    // calling either `setup_euclidean_maps(m)` first or the accessor first
    // is equivalent.

+    /// Return the built-in vertex-point map of `m` (the geometric embedding).
    static Vertex_point_map vertex_points(Triangle_mesh& m) {
        return m.points();
    }

+    /// Return (or create with default 2π) the target-angle property map.
    static Theta_pmap theta_map(Triangle_mesh& m) {
        auto [pm, created] = m.template add_property_map<Vertex_descriptor, FT>(
            "ev:theta", FT(2.0 * 3.141592653589793238));
@@ -144,6 +162,7 @@ struct Default_conformal_map_traits<CGAL::Surface_mesh<typename K::Point_3>, K>
        return pm;
    }

+    /// Return (or create with default −1) the vertex-index property map.
    static Vertex_index_pmap vertex_index_map(Triangle_mesh& m) {
        auto [pm, created] = m.template add_property_map<Vertex_descriptor, int>(
            "ev:idx", -1);
@@ -151,6 +170,7 @@ struct Default_conformal_map_traits<CGAL::Surface_mesh<typename K::Point_3>, K>
        return pm;
    }

+    /// Return (or create with default 0) the λ⁰ (initial log-length) property map.
    static Lambda0_pmap lambda0_map(Triangle_mesh& m) {
        auto [pm, created] = m.template add_property_map<Edge_descriptor, FT>(
            "ee:lam0", FT(0));
--- a/code/include/CGAL/Discrete_circle_packing.h
+++ b/code/include/CGAL/Discrete_circle_packing.h
@@ -38,26 +38,51 @@ namespace CGAL {

 // ── Default traits for CP-Euclidean ───────────────────────────────────────────

+/*!
+\ingroup PkgConformalMapConcepts
+\brief Traits class for `discrete_circle_packing_euclidean()` —
+declares the kernel, mesh and property-map types used by the
+BPS-2010 face-based circle-packing functional.
+
+Primary template; specialise it for non-`Surface_mesh` triangle meshes.
+*/
 template <typename TriangleMesh,
          typename Kernel_ = CGAL::Simple_cartesian<double>>
 struct Default_cp_euclidean_traits;

+/*!
+\ingroup PkgConformalMapConcepts
+\brief Specialisation for `CGAL::Surface_mesh<P>`; the only one shipped
+in Phase 8b-Lite.
+*/
 template <typename K>
 struct Default_cp_euclidean_traits<CGAL::Surface_mesh<typename K::Point_3>, K>
 {
+    /// CGAL kernel parameter (defaults to `Simple_cartesian<double>`).
    using Kernel        = K;
+    /// Scalar field type used for all CP-Euclidean DOFs (`ρ_f`, `θ_e`, `φ_f`).
    using FT            = typename K::FT;
+    /// 3-D point type (vertex coordinates).
    using Point_3       = typename K::Point_3;
+    /// Triangle-mesh type this specialisation targets.
    using Triangle_mesh = CGAL::Surface_mesh<Point_3>;

+    /// Boost-graph vertex descriptor for `Triangle_mesh`.
    using Vertex_descriptor   = typename boost::graph_traits<Triangle_mesh>::vertex_descriptor;
+    /// Boost-graph half-edge descriptor for `Triangle_mesh`.
    using Halfedge_descriptor = typename boost::graph_traits<Triangle_mesh>::halfedge_descriptor;
+    /// Boost-graph edge descriptor for `Triangle_mesh`.
    using Edge_descriptor     = typename boost::graph_traits<Triangle_mesh>::edge_descriptor;
+    /// Boost-graph face descriptor for `Triangle_mesh`.
    using Face_descriptor     = typename boost::graph_traits<Triangle_mesh>::face_descriptor;

    // CP-Euclidean property maps — note the *face* DOF index map.
+
+    /// Property map face → contiguous integer DOF index (legacy `cf:idx`).
    using Face_index_pmap = typename Triangle_mesh::template Property_map<Face_descriptor,   int>;
+    /// Property map edge → intersection angle θₑ (legacy `ce:theta`).
    using Theta_e_pmap    = typename Triangle_mesh::template Property_map<Edge_descriptor,   FT>;
+    /// Property map face → target angle sum φ_f (legacy `cf:phi`).
    using Phi_f_pmap      = typename Triangle_mesh::template Property_map<Face_descriptor,   FT>;
 };

@@ -75,8 +100,11 @@ struct Circle_packing_result
    /// Face DOFs `ρ_f = log R_f` (length = num_faces(mesh); pinned face = 0).
    std::vector<FT> rho_per_face;

+    /// Newton iterations actually performed (≤ `max_iterations`).
    int  iterations    = 0;
+    /// Final infinity-norm of the gradient (Newton stopping criterion).
    FT   gradient_norm = FT(0);
+    /// `true` iff `gradient_norm < gradient_tolerance` at exit.
    bool converged     = false;
 };

--- a/code/include/CGAL/Discrete_conformal_map.h
+++ b/code/include/CGAL/Discrete_conformal_map.h
@@ -449,8 +449,11 @@ struct Hyper_ideal_map_result
    /// Edge DOFs `a_e` (length = num_edges(mesh); pinned edges = 0).
    std::vector<FT> a_per_edge;

+    /// Newton iterations actually performed (≤ `max_iterations`).
    int  iterations    = 0;
+    /// Final infinity-norm of the gradient (Newton stopping criterion).
    FT   gradient_norm = FT(0);
+    /// `true` iff `gradient_norm < gradient_tolerance` at exit.
    bool converged     = false;
 };

--- a/code/include/CGAL/Discrete_inversive_distance.h
+++ b/code/include/CGAL/Discrete_inversive_distance.h
@@ -47,25 +47,49 @@ namespace CGAL {

 // ── Default traits for Inversive-Distance ────────────────────────────────────

+/*!
+\ingroup PkgConformalMapConcepts
+\brief Traits class for `discrete_inversive_distance_map()` — declares
+the kernel, mesh and property-map types used by Luo's 2004 vertex-based
+inversive-distance circle packing.
+
+Primary template; specialise it for non-`Surface_mesh` triangle meshes.
+*/
 template <typename TriangleMesh,
          typename Kernel_ = CGAL::Simple_cartesian<double>>
 struct Default_inversive_distance_traits;

+/*!
+\ingroup PkgConformalMapConcepts
+\brief Specialisation for `CGAL::Surface_mesh<P>`; the only one shipped
+in Phase 8b-Lite.
+*/
 template <typename K>
 struct Default_inversive_distance_traits<CGAL::Surface_mesh<typename K::Point_3>, K>
 {
+    /// CGAL kernel parameter (defaults to `Simple_cartesian<double>`).
    using Kernel        = K;
+    /// Scalar field type used for all inversive-distance DOFs.
    using FT            = typename K::FT;
+    /// 3-D point type (vertex coordinates).
    using Point_3       = typename K::Point_3;
+    /// Triangle-mesh type this specialisation targets.
    using Triangle_mesh = CGAL::Surface_mesh<Point_3>;

+    /// Boost-graph vertex descriptor for `Triangle_mesh`.
    using Vertex_descriptor   = typename boost::graph_traits<Triangle_mesh>::vertex_descriptor;
+    /// Boost-graph edge descriptor for `Triangle_mesh`.
    using Edge_descriptor     = typename boost::graph_traits<Triangle_mesh>::edge_descriptor;

    // Inversive-distance specific property maps.
+
+    /// Property map vertex → contiguous integer DOF index (legacy `iv:idx`).
    using Vertex_index_pmap = typename Triangle_mesh::template Property_map<Vertex_descriptor, int>;
+    /// Property map vertex → target cone angle Θᵥ in radians (legacy `iv:theta`).
    using Theta_v_pmap      = typename Triangle_mesh::template Property_map<Vertex_descriptor, FT>;
+    /// Property map vertex → initial radius r⁰ᵥ (legacy `iv:r0`).
    using R0_pmap           = typename Triangle_mesh::template Property_map<Vertex_descriptor, FT>;
+    /// Property map edge → inversive distance Iᵢⱼ (legacy `ie:I`).
    using I_e_pmap          = typename Triangle_mesh::template Property_map<Edge_descriptor,   FT>;
 };

--- a/code/include/clausen.hpp
+++ b/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)
--- a/code/include/conformal_mesh.hpp
+++ b/code/include/conformal_mesh.hpp
@@ -41,30 +41,42 @@ namespace conformallab {
 // ── Kernel ──────────────────────────────────────────────────────────────────
 // Simple double-precision Cartesian.  Conformal mapping algorithms never
 // need exact arithmetic — they operate on floating-point lengths and angles.
+
+/// CGAL kernel used by all conformallab algorithms (double precision).
 using Kernel = CGAL::Simple_cartesian<double>;
+/// 3-D point type (vertex coordinates).
 using Point3 = Kernel::Point_3;
+/// 2-D point type (UV-domain layout coordinates).
 using Point2 = Kernel::Point_2;

 // ── Mesh type ────────────────────────────────────────────────────────────────
+/// Triangle mesh carrying all conformal-map data as property maps.
 using ConformalMesh = CGAL::Surface_mesh<Point3>;

 // ── Index/descriptor aliases (CGAL 6.x naming) ───────────────────────────────
+/// Vertex descriptor of `ConformalMesh`.
 using Vertex_index   = ConformalMesh::Vertex_index;
+/// Half-edge descriptor of `ConformalMesh`.
 using Halfedge_index = ConformalMesh::Halfedge_index;
+/// Edge descriptor of `ConformalMesh`.
 using Edge_index     = ConformalMesh::Edge_index;
+/// Face descriptor of `ConformalMesh`.
 using Face_index     = ConformalMesh::Face_index;

 // ── Geometry type constant (replaces Java CoFace.type enum) ─────────────────
+/// Discrete geometry type that a face/mesh is interpreted in.
+/// Replaces the original Java `CoFace.type` enum.
 enum class GeometryType : int {
-    Euclidean  = 0,
-    Hyperbolic = 1,
-    Spherical  = 2
+    Euclidean  = 0, ///< Flat metric (ℝ²).
+    Hyperbolic = 1, ///< Hyperbolic metric (ℍ²).
+    Spherical  = 2  ///< Spherical metric (S²).
 };

 // ── Standard property-map bundles ────────────────────────────────────────────

-// Add the vertex properties used by all conformal-map functionals.
-// Returns {lambda, theta, idx}.
+/// Register and return the vertex-side property maps used by all five
+/// DCE functionals: `v:lambda` (log conformal factor), `v:theta` (target
+/// cone angle), `v:idx` (contiguous integer index).
 inline auto add_vertex_properties(ConformalMesh& mesh)
 {
    auto [lambda, ok1] = mesh.add_property_map<Vertex_index, double>("v:lambda", 0.0);
@@ -74,7 +86,8 @@ inline auto add_vertex_properties(ConformalMesh& mesh)
    return std::make_tuple(lambda, theta, idx);
 }

-// Add the edge intersection-angle property used by the hyperbolic functional.
+/// Register and return the edge intersection-angle property `e:alpha`
+/// (used by the hyper-ideal functional).
 inline auto add_edge_properties(ConformalMesh& mesh)
 {
    auto [alpha, ok] = mesh.add_property_map<Edge_index, double>("e:alpha", 0.0);
@@ -82,7 +95,7 @@ inline auto add_edge_properties(ConformalMesh& mesh)
    return alpha;
 }

-// Add the face geometry-type property.
+/// Register and return the per-face geometry-type property `f:type`.
 inline auto add_face_properties(ConformalMesh& mesh)
 {
    auto [ftype, ok] = mesh.add_property_map<Face_index, int>(
--- a/code/include/cp_euclidean_functional.hpp
+++ b/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,
--- a/code/include/cut_graph.hpp
+++ b/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();
--- a/code/include/discrete_elliptic_utility.hpp
+++ b/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) {
--- a/code/include/euclidean_functional.hpp
+++ b/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,
--- a/code/include/euclidean_geometry.hpp
+++ b/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)
 {
--- a/code/include/euclidean_hessian.hpp
+++ b/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,
--- a/code/include/fundamental_domain.hpp
+++ b/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,
--- a/code/include/gauss_bonnet.hpp
+++ b/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)
 {
--- a/code/include/hyper_ideal_functional.hpp
+++ b/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,
--- a/code/include/hyper_ideal_geometry.hpp
+++ b/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,
--- a/code/include/hyper_ideal_hessian.hpp
+++ b/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,
--- a/code/include/hyper_ideal_utility.hpp
+++ b/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)
--- a/code/include/hyper_ideal_visualization_utility.hpp
+++ b/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)
 {
--- a/code/include/inversive_distance_functional.hpp
+++ b/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,
--- a/code/include/layout.hpp
+++ b/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)
 {
--- a/code/include/matrix_utility.hpp
+++ b/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();
--- a/code/include/mesh_builder.hpp
+++ b/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;
--- a/code/include/mesh_io.hpp
+++ b/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))
--- a/code/include/newton_solver.hpp
+++ b/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,
--- a/code/include/p2_utility.hpp
+++ b/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,
--- a/code/include/period_matrix.hpp
+++ b/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;
--- a/code/include/projective_math.hpp
+++ b/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,
--- a/code/include/serialization.hpp
+++ b/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,
--- a/code/include/spherical_functional.hpp
+++ b/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,
--- a/code/include/spherical_geometry.hpp
+++ b/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;
--- a/code/include/spherical_hessian.hpp
+++ b/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,
--- a/code/include/viewer_utils.h
+++ b/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);

 } 
--- a/doc/api/headers.md
+++ b/doc/api/headers.md
@@ -1,111 +1,69 @@
+<!-- AUTO-GENERATED by scripts/gen-headers-md.py — do not edit by hand. -->
+<!-- Source of truth: the `\file` / leading-comment briefs in code/include/. -->
+
 # Public Headers (`code/include/`)

-All algorithms are header-only. Include the headers you need directly —
-there is no compiled library to link against (only GTest for tests and
-CGAL/Eigen for the CGAL-dependent headers).
+All algorithms are header-only.  Include the headers you need
+directly — there is no compiled library to link against (only
+GTest for tests and CGAL/Eigen for the CGAL-dependent headers).

-## Core mesh type
+This page is **regenerated** from Doxygen XML on every push to
+`main` that touches a public header (see
+`.gitea/workflows/doxygen-pages.yml` and `scripts/gen-headers-md.py`).
+To improve a description, edit the `\file` brief at the top of
+the corresponding header, then re-run `bash scripts/regen-docs.sh`.

-| Header | Description |
-|---|---|
-| `conformal_mesh.hpp` | `ConformalMesh` = `CGAL::Surface_mesh<Point3>`. Property-map naming convention. Index type aliases. `GeometryType` enum. |
-| `constants.hpp` | `conformallab::PI`, `TWO_PI` |
-| `mesh_builder.hpp` | `make_triangle()` / `make_tetrahedron()` / `make_quad_strip()` / `make_fan()` / `make_open_cylinder()` — test mesh factories |
-| `mesh_io.hpp` | `load_mesh()` / `save_mesh()` via `CGAL::IO` (OFF / OBJ / PLY) |
-| `mesh_utils.hpp` | `cgal_to_eigen()` — convert `ConformalMesh` vertex positions to `Eigen::MatrixXd` |
+## CGAL public API (`<CGAL/...>`)

-## Special functions
+| Header | Brief | Public symbols |
+|---|---|---|
+| `CGAL/Conformal_layout.h` | Thin CGAL-style wrapper around the legacy euclidean_layout(), spherical_layout() and hyper_ideal_layout() functions defined in code/include/layout.hpp. | `Layout2D`, `Layout3D`, `HolonomyData`, `CutGraph` |
+| `CGAL/Conformal_map_traits.h` | Defines the ConformalMapTraits concept and the default model Default_conformal_map_traits<TriangleMesh, K> for the package. | _(no public symbols)_ |
+| `CGAL/Discrete_circle_packing.h` | User-facing entry for the face-based circle-packing functional of Bobenko-Pinkall-Springborn 2010. | `Circle_packing_result` |
+| `CGAL/Discrete_conformal_map.h` | User-facing entry point for the Discrete_conformal_map package. | `Conformal_map_result`, `Hyper_ideal_map_result` |
+| `CGAL/Discrete_inversive_distance.h` | User-facing entry for the vertex-based inversive-distance circle- packing functional of Luo (2004), with the Bowers-Stephenson (2004) initialisation. | _(no public symbols)_ |

-| Header | Description |
-|---|---|
-| `clausen.hpp` | `clausen_cl2(θ)` (Clausen Cl₂), `lobachevsky(θ)` (Л), `im_li2(θ)` (ImLi₂ = imaginary part of dilogarithm) |
+## CGAL internals (`<CGAL/Conformal_map/...>`)

-## HyperIdeal geometry (H²)
+| Header | Brief | Public symbols |
+|---|---|---|
+| `CGAL/Conformal_map/doxygen_groups.h` | Copyright (c) 2024-2026 Tarik Moussa. | _(no public symbols)_ |
+| `CGAL/Conformal_map/doxygen_namespaces.h` | Copyright (c) 2024-2026 Tarik Moussa. | _(no public symbols)_ |
+| `CGAL/Conformal_map/internal/parameters.h` | Copyright (c) 2024-2026 Tarik Moussa. | _(no public symbols)_ |

-| Header | Description |
-|---|---|
-| `hyper_ideal_geometry.hpp` | `ζ₁₃`, `ζ₁₄`, `ζ₁₅` (Springborn 2020), `l_from_zeta()`, `alpha_ij()`, `beta_i()`, `sigma_i()`, `sigma_ij()` |
-| `hyper_ideal_utility.hpp` | Tetrahedron volumes (Meyerhoff formula, Kolpakov–Mednykh) |
-| `hyper_ideal_visualization_utility.hpp` | Poincaré disk projection, circumcircle helpers, Lorentz boost |
-| `hyper_ideal_functional.hpp` | `HyperIdealMaps`, `setup_hyper_ideal_maps()`, `compute_hyper_ideal_lambda0_from_mesh()`, `assign_all_dof_indices()`, `hyper_ideal_gradient()`, `hyper_ideal_energy()` |
-| `hyper_ideal_hessian.hpp` | `hyper_ideal_hessian()` — symmetric FD Hessian (Phase 9b: analytic) |
+## Core (`conformallab` namespace, `<...>`)

-## Spherical geometry (S²)
+| Header | Brief | Public symbols |
+|---|---|---|
+| `clausen.hpp` | Clausen integral, Lobachevsky function, and Im(Li2). | _(no public symbols)_ |
+| `conformal_mesh.hpp` | conformal_mesh.hpp Central mesh type for the discrete conformal mapping algorithms. | _(no public symbols)_ |
+| `constants.hpp` | constants.hpp Single source of truth for mathematical constants used throughout conformallab++. | _(no public symbols)_ |
+| `cp_euclidean_functional.hpp` | cp_euclidean_functional.hpp Phase 9a.1 — Circle-Packing Euclidean functional (CP-Euclidean). | `CPEuclideanMaps` |
+| `cut_graph.hpp` | cut_graph.hpp Phase 6 — Tree-cotree algorithm for computing a cut graph of a triangulated surface. | `CutGraph` |
+| `discrete_elliptic_utility.hpp` | Ported from de.varylab.discreteconformal.util.DiscreteEllipticUtility (Java). | _(no public symbols)_ |
+| `euclidean_functional.hpp` | euclidean_functional.hpp Energy and gradient of the Euclidean discrete conformal functional (EuclideanCyclicFunctional) evaluated on a ConformalMesh. | `EuclideanMaps`, `EuclideanResult` |
+| `euclidean_geometry.hpp` | euclidean_geometry.hpp Corner-angle formula for Euclidean triangles in the discrete conformal (log-length) parametrisation. | `EuclideanFaceAngles` |
+| `euclidean_hessian.hpp` | euclidean_hessian.hpp Analytical Hessian of the Euclidean discrete conformal energy — the cotangent-Laplace operator. | `EuclCotWeights` |
+| `fundamental_domain.hpp` | fundamental_domain.hpp Phase 7 — Fundamental domain polygon for closed surfaces. | `FundamentalDomain` |
+| `gauss_bonnet.hpp` | gauss_bonnet.hpp Phase 6 — Gauss–Bonnet consistency check for prescribed target angles. | _(no public symbols)_ |
+| `hyper_ideal_functional.hpp` | hyper_ideal_functional.hpp Energy and gradient of the hyper-ideal discrete conformal map functional evaluated on a ConformalMesh (CGAL::Surface_mesh). | `HyperIdealMaps`, `HyperIdealResult`, `FaceAngleOutputs`, `FaceAngles` |
+| `hyper_ideal_geometry.hpp` | hyper_ideal_geometry.hpp Pure-math building blocks for the hyper-ideal discrete conformal map. | _(no public symbols)_ |
+| `hyper_ideal_hessian.hpp` | hyper_ideal_hessian.hpp Phase 4a — Hessian of the hyper-ideal discrete conformal functional. | _(no public symbols)_ |
+| `hyper_ideal_utility.hpp` | Hyperbolic tetrahedron volume formulas. | _(no public symbols)_ |
+| `hyper_ideal_visualization_utility.hpp` | Port of the static helper HyperIdealVisualizationPlugin.getEuclideanCircleFromHyperbolic() from de.varylab.discreteconformal.plugin. | _(no public symbols)_ |
+| `inversive_distance_functional.hpp` | inversive_distance_functional.hpp Phase 9a.2 — Inversive-distance circle-packing functional (Luo 2004). | `InversiveDistanceMaps` |
+| `layout.hpp` | layout.hpp Phase 5/6/7 — Layout / embedding: DOF vector → vertex coordinates in the target geometry via BFS-trilateration. | `MobiusMap`, `Layout2D`, `Layout3D`, `HolonomyData` |
+| `matrix_utility.hpp` | 4x4 mapping matrix from corresponding point pairs. | _(no public symbols)_ |
+| `mesh_builder.hpp` | mesh_builder.hpp Factory functions that build simple reference meshes for testing and examples. | _(no public symbols)_ |
+| `mesh_io.hpp` | mesh_io.hpp Phase 4b — CGAL::IO wrappers for ConformalMesh. | _(no public symbols)_ |
+| `mesh_utils.hpp` | mesh_utils.hpp Conversions between CGAL::Surface_mesh and Eigen matrices. | _(no public symbols)_ |
+| `newton_solver.hpp` | newton_solver.hpp Phase 4a — Newton solver for all three discrete conformal functionals. | `NewtonResult` |
+| `p2_utility.hpp` | 2-D projective geometry utilities for the Euclidean signature. | _(no public symbols)_ |
+| `period_matrix.hpp` | period_matrix.hpp Phase 7 — Period matrix for closed surfaces with Euclidean (flat) metric. | `PeriodData` |
+| `projective_math.hpp` | Projective and hyperbolic geometry utilities. | _(no public symbols)_ |
+| `serialization.hpp` | serialization.hpp Phase 5 — Save and load conformal map results in JSON and XML formats. | _(no public symbols)_ |
+| `spherical_functional.hpp` | spherical_functional.hpp Energy and gradient of the spherical discrete conformal functional evaluated on a ConformalMesh (CGAL::Surface_mesh). | `SphericalMaps`, `SphericalResult` |
+| `spherical_geometry.hpp` | spherical_geometry.hpp Pure-math building blocks for the spherical discrete conformal map. | `SphericalFaceAngles` |
+| `spherical_hessian.hpp` | spherical_hessian.hpp Analytical Hessian of the spherical discrete conformal energy — the spherical cotangent-Laplace operator. | `SpherCotWeights` |
+| `viewer_utils.h` | _(undocumented — add a `\file` brief at the top of the header)_ | _(no public symbols)_ |

-| Header | Description |
-|---|---|
-| `spherical_geometry.hpp` | Spherical arc-length, half-angle formula, spherical law of cosines |
-| `spherical_functional.hpp` | `SphericalMaps`, `setup_spherical_maps()`, `compute_spherical_lambda0_from_mesh()`, `spherical_gradient()`, `spherical_energy()` |
-| `spherical_hessian.hpp` | `spherical_hessian()` — analytic Hessian via ∂α/∂u from law of cosines |
-
-## Euclidean geometry (ℝ²)
-
-| Header | Description |
-|---|---|
-| `euclidean_geometry.hpp` | Euclidean corner angle (t-value / atan2), edge lengths from DOF vector |
-| `euclidean_functional.hpp` | `EuclideanMaps`, `setup_euclidean_maps()`, `compute_euclidean_lambda0_from_mesh()`, `euclidean_gradient()`, `euclidean_energy()` |
-| `euclidean_hessian.hpp` | `euclidean_hessian()` — cotangent Laplacian (Pinkall–Polthier 1993) |
-
-## Circle-packing functionals (Phase 9a)
-
-| Header | Description |
-|---|---|
-| `cp_euclidean_functional.hpp` | **Face-based** CP-Euclidean (Bobenko–Pinkall–Springborn 2010), port of Java `CPEuclideanFunctional`.  `CPEuclideanMaps`, `setup_cp_euclidean_maps()`, `assign_cp_euclidean_face_dof_indices()`, `cp_euclidean_gradient()`, `cp_euclidean_energy()`, **analytic** `cp_euclidean_hessian()` (2×2-per-edge `h_jk = sin θ / (cosh Δρ − cos θ)`) |
-| `inversive_distance_functional.hpp` | **Vertex-based** inversive-distance (Luo 2004, no Java original).  `InversiveDistanceMaps`, `setup_inversive_distance_maps()`, `compute_inversive_distance_init_from_mesh()` (Bowers–Stephenson 2004 init), `inversive_distance_gradient()`, `inversive_distance_energy()` |
-
-## Solver
-
-| Header | Description |
-|---|---|
-| `newton_solver.hpp` | Five Newton solvers — `newton_euclidean()`, `newton_spherical()`, `newton_hyper_ideal()` (uses block-FD Hessian since Phase 9b), `newton_cp_euclidean()` (analytic Hessian, Phase 9a-Newton), `newton_inversive_distance()` (FD Hessian, Phase 9a-Newton). Plus `solve_linear_system()` (SimplicialLDLT + SparseQR fallback) and the `NewtonResult` struct. |
-
-## Preprocessing
-
-| Header | Description |
-|---|---|
-| `gauss_bonnet.hpp` | `euler_characteristic()`, `genus()`, `gauss_bonnet_sum()`, `gauss_bonnet_rhs()`, `gauss_bonnet_deficit()`, `check_gauss_bonnet()`, `enforce_gauss_bonnet()` |
-
-## Layout and holonomy
-
-| Header | Description |
-|---|---|
-| `cut_graph.hpp` | `CutGraph` struct, `compute_cut_graph()` — tree-cotree algorithm (Erickson–Whittlesey 2005), produces 2g seam edges |
-| `layout.hpp` | `euclidean_layout()`, `spherical_layout()`, `hyper_ideal_layout()`, `normalise_{euclidean,hyperbolic,spherical}()`. Structs: `Layout2D`, `Layout3D`, `HolonomyData`. `MobiusMap` (T(z)=(az+b)/(cz+d), `from_three`, `compose`, `inverse`, `apply`). Priority-BFS, `halfedge_uv`. |
-
-## Post-processing
-
-| Header | Description |
-|---|---|
-| `period_matrix.hpp` | `PeriodData`, `compute_period_matrix()`, `reduce_to_fundamental_domain()`, `is_in_fundamental_domain()` — period ratio τ = ω₂/ω₁ ∈ ℍ, SL(2,ℤ) reduction |
-| `fundamental_domain.hpp` | `FundamentalDomain`, `compute_fundamental_domain()` (genus 1: CCW parallelogram; genus > 1: empty, TODO Phase 9c), `tiling_copy()`, `tiling_neighbourhood()` |
-| `serialization.hpp` | `save_result_json()`, `load_result_json()`, `save_result_xml()`, `load_result_xml()`, `save_layout_off()` |
-
-## Math utilities
-
-These headers contain pure-math helpers used internally across the
-package.  All are also re-exported on the public API so downstream
-consumers (e.g. user code that wants to assemble a custom functional)
-can use them directly.
-
-| Header | Description |
-|---|---|
-| `matrix_utility.hpp` | Small linear-algebra helpers used by `period_matrix.hpp` (2×2 determinant, …) |
-| `projective_math.hpp` | Real projective 2-space helpers (point on line, lines through points, dual) |
-| `p2_utility.hpp` | Higher-level P² utilities used by hyper-ideal layout (perpendicular bisectors in Poincaré disk, etc.) |
-| `discrete_elliptic_utility.hpp` | Genus-1 / elliptic-curve helpers — half-period ratio τ from texture coords; bridge to genus-1 fundamental domain |
-
-## CGAL public API (Phase 8b-Lite)
-
-These headers live under `include/CGAL/` and are the user-facing entry
-points.  They are thin wrappers over the implementation in
-`code/include/*.hpp`, exposing the five DCE models through a
-CGAL-conventional interface.
-
-| Header | Description |
-|---|---|
-| `CGAL/Conformal_map_traits.h` | `ConformalMapTraits` concept + `Default_conformal_map_traits<TriangleMesh, K>` (Euclidean-flavoured base trait) |
-| `CGAL/Discrete_conformal_map.h` | `discrete_conformal_map_euclidean()`, `discrete_conformal_map_spherical()`, `discrete_conformal_map_hyper_ideal()`.  Result types `Conformal_map_result<FT>` and `Hyper_ideal_map_result<FT>` |
-| `CGAL/Discrete_circle_packing.h` | `Default_cp_euclidean_traits` + `discrete_circle_packing_euclidean()`.  Result: `Circle_packing_result<FT>` with `rho_per_face` |
-| `CGAL/Discrete_inversive_distance.h` | `Default_inversive_distance_traits` + `discrete_inversive_distance_map()`.  Result: `Conformal_map_result<FT>` reused |
-| `CGAL/Conformal_layout.h` | Thin re-export of `euclidean_layout`, `spherical_layout`, `hyper_ideal_layout` into the `CGAL::` namespace |
-| `CGAL/Conformal_map/internal/parameters.h` | (internal) Named-parameter tags for `CGAL::parameters::*` |
--- a/doc/api/tests.md
+++ b/doc/api/tests.md
@@ -45,7 +45,7 @@ All tests have CTest prefix `cgal.` (set via `TEST_PREFIX "cgal."` in CMakeLists
 | `Normalisation` | `test_phase6.cpp` | 4 | Euclidean centroid, length ratios, Möbius centring |
 | `MobiusMap` | `test_phase7.cpp` | 8 | Identity, inverse, compose, `from_three`, `apply(Vector2d)` |
 | `BestRootFace` | `test_phase7.cpp` | 2 | Valid root face selection, interior bonus |
-| `HalfedgeUV` | `test_phase7.cpp` | 4 | Size = #halfedges, seam consistency, boundary halfedges = 0 |
+| `HalfedgeUV` | `test_phase7.cpp` | 4 | Size = number of half-edges, seam consistency, boundary half-edges = 0 |
 | `PriorityBFS` | `test_phase7.cpp` | 3 | Success, no seam on open meshes, all vertices placed |
 | `NormaliseEuclidean` | `test_phase7.cpp` | 2 | UV centroid = 0, halfedge_uv centroid = 0 |
 | `PeriodMatrix` | `test_phase7.cpp` | 7 | τ ∈ ℍ, SL(2,ℤ) reduction, exception outside ℍ |
--- a/doc/architecture/locked-vs-flexible.md
+++ b/doc/architecture/locked-vs-flexible.md
@@ -167,7 +167,7 @@ the CGAL-canonical syntax AND we are willing to fork CGAL upstream.
 | Status         | 🟢 opportunistic |
 | Locked since   | Phase 3 + 9a (prefix `ev:`/`sv:`/`v:`/`cf:`/`ce:`/`iv:`/`ie:` set when each functional was introduced) |
 | Cost to change | One sed-replace + recompile.  No user-visible effect because the names are an *internal* convention; the CGAL public API never exposes them. |
-| When to revisit | If a future functional reuses an existing letter prefix.  Already discussed in [`locked-vs-flexible.md`](#4-five-dce-models-on-the-same-mesh). |
+| When to revisit | If a future functional reuses an existing letter prefix.  Already discussed in §4 "Five DCE models on the same mesh" above. |

 ---

--- a/doc/architecture/overall_pipeline.md
+++ b/doc/architecture/overall_pipeline.md
@@ -150,7 +150,7 @@ maps.theta_v[v] = M_PI / 3;          // 60° cone singularity

 Before solving, the prescribed angles must satisfy:

-$$\sum_{v} (2\pi - \Theta_v) = 2\pi \cdot \chi(M)$$
+> ∑ᵥ (2π − Θᵥ) = 2π · χ(M)

 ```cpp
 check_gauss_bonnet(mesh, maps);         // throws if violated
@@ -283,7 +283,7 @@ Both `uv` and `halfedge_uv` are transformed identically.
 From the two holonomy translations ω₁, ω₂ ∈ ℂ read off from the cut graph,
 the conformal type of a flat torus is the SL(2,ℤ)-orbit of:

-$$\tau = \omega_2 / \omega_1 \in \mathbb{H}$$
+> τ = ω₂ / ω₁ ∈ ℍ

 ```cpp
 PeriodData pd = compute_period_matrix(hol);
--- a/doc/roadmap/java-parity.md
+++ b/doc/roadmap/java-parity.md
@@ -67,14 +67,91 @@ They are candidates for Phase 9 or Phase 10.

 Note: items marked as *new research* (e.g. Inversive Distance, HyperIdeal Hessian variants)
 are tracked separately in `doc/roadmap/research-track.md`.
-| `HomotopyUtility` | Homotopy generators | 9c |
-| `SpanningTreeUtility` | Spanning tree algorithms | 8 / infrastructure |
-| `SurgeryUtility` | Mesh surgery (cut/glue) | — |
-| `StitchingUtility` | Seam stitching | — |
-| `CuttingUtility` | Advanced cutting (beyond tree-cotree) | 9c |
-| `HyperellipticUtility` | Hyperelliptic surfaces | 10 |
-| `LaplaceUtility` | Discrete Laplace operators | 9 / infrastructure |
-| `ConformalStructureUtility` | Conformal structure extraction | 10 |
+
+---
+
+## Java packages not yet in the roadmap — 2026 full-library scan
+
+A complete scan of the Java source tree (`de.varylab.discreteconformal`) revealed the
+following packages and classes not yet covered by Phases 1–9c or the Phase-10 plan.
+Organised by value / effort.
+
+### Functional package (`de.varylab.discreteconformal.functional`)
+
+| Java class | What it does | Proposed phase |
+|---|---|---|
+| `ConesUtility` (in `unwrapper/`) | Cone singularity detection, BFS-path cutting from cone to boundary, auto-placement via conjugate gradient, quantization to π/2 / π/3 / π/6 — fills the "⚠️ data structure only" gap in the parity table | **9d** |
+| `MobiusCenteringFunctional` | Variational Möbius centering via Lorentz geometry: E = Σ log(−⟨x,p⟩/√(−⟨x,x⟩)). Supplies full gradient + Hessian — more principled than the iterative Fréchet mean in `normalise_hyperbolic()` | **9d** |
+| `ElectrostaticSphereFunctional` | Repulsive electrostatic energy on S² (E = Σ 0.5/d² + sphere constraint). Useful as initialization heuristic for spherical uniformization | **9d** (optional) |
+| `EuclideanCyclicFunctional` | Euclidean DCE functional reduced to a cyclic-symmetry quotient — reduces DOFs for surfaces with cyclic symmetry group | **10g** |
+| `HyperbolicCyclicFunctional` | Hyperbolic analogue of above | **10g** |
+
+### Circle pattern package (`de.varylab.discreteconformal.unwrapper.circlepattern`)
+
+| Java class | What it does | Proposed phase |
+|---|---|---|
+| `CirclePatternUtility` | Computes circle pattern radii (ρ per face) via Newton trust-region on CPEuclidean energy — the solver side | **9e** |
+| `CirclePatternLayout` | Embeds a circle pattern in the plane from the ρ values — the layout side. Required complement to `cp_euclidean_functional.hpp` already ported in 9a.1 | **9e** |
+| `CPEuclideanRotation` | Rotation-invariant variant of the CP-Euclidean functional | **9e** |
+
+### Uniformization package (`de.varylab.discreteconformal.uniformization`)
+
+| Java class | What it does | Proposed phase |
+|---|---|---|
+| `CutAndGlueUtility` | Mesh surgery beyond tree-cotree: arbitrary cut paths, gluing cut surfaces back together — needed for genus g > 1 canonical polygon construction | **9c** (add to existing plan) |
+| `VisualizationUtility` | Java/jReality specific — do not port | — |
+| `Uniformizer` | High-level pipeline driver — already covered by our CLI + pipeline API | — |
+
+### Unwrapper package — additional classes
+
+| Java class | What it does | Proposed phase |
+|---|---|---|
+| `StereographicUnwrapper` | Stereographic projection S²→ℂ∪{∞} + Möbius centring — converts spherical DCE output to 2-D atlas for genus-0 | **9d** |
+| `SphereUtility` | Sphere-specific utilities (area centroid, antipodal, normalization helpers) | **9d** |
+| `CircleDomainUnwrapper` | Conformal map of multiply-connected planar region onto disk-with-holes (Koebe-Andreev-Thurston) | **10d** |
+
+### Quasi-isothermic package (`de.varylab.discreteconformal.unwrapper.quasiisothermic`)
+
+Quasi-isothermic maps generalize conformal maps to meshes where exact conformality is impossible.
+Entirely absent from the current roadmap.
+
+| Java class | What it does | Proposed phase |
+|---|---|---|
+| `QuasiisothermicDelaunay` | Delaunay-conformal triangulation as QI preprocessing | **10e** |
+| `QuasiisothermicLayout` | Embedding from quasi-isothermic DOFs | **10e** |
+| `QuasiisothermicUtility` | Lawson-correspondence parameterization (~800 lines) | **10e** |
+| `DBFSolution` | Discrete Beltrami field solution | **10e** |
+| `SinConditionApplication` | Sin-condition functional for QI maps | **10e** |
+| `ConformalStructureUtility` | Conformal structure extraction from QI solution | **10e** |
+
+### Koebe package (`de.varylab.discreteconformal.unwrapper.koebe`)
+
+| Java class | What it does | Proposed phase |
+|---|---|---|
+| `KoebePolyhedron` | Koebe–Andreev–Thurston theorem: realization of every 3-connected planar graph as a convex polyhedron with edges tangent to the unit sphere (321 lines) | **10f** |
+
+### Util package — additional classes not yet in roadmap
+
+| Java class | What it does | Proposed phase |
+|---|---|---|
+| `DualityUtility` | Hodge-star operator + dual cycles via cotangent weights — **prerequisite for Phase 10a** (discrete holomorphic forms need the dual mesh) | **10a prerequisite** (consider 9c) |
+| `HyperellipticUtility` | Hyperelliptic surfaces (genus g ≥ 2 with Z₂ symmetry) — period matrix has block-diagonal structure | **10b** |
+| `HyperIdealHyperellipticUtility` | HyperIdeal variant for hyperelliptic surfaces | **10b** |
+| `PathUtility` | Mesh path operations — supporting infrastructure for Phase 9c homology basis | **9c** |
+| `LaplaceUtility` | Discrete Laplace operators (cotangent + combinatorial) | **9 / infrastructure** |
+| `EdgeUtility` | Edge orientation and classification helpers | **infrastructure** |
+| `StitchingUtility` | Seam stitching after cut-and-glue | **9c** |
+
+### Do not port — Java-specific or superseded
+
+| Java class | Reason |
+|---|---|
+| `ColtIterationReporterImpl` | Colt sparse matrix library — replaced by Eigen |
+| `NodeIndexComparator` | Java Comparator — replaced by `std::less` |
+| `SimpleMatrixPrintUtility` | Debug print — Eigen `.format()` is sufficient |
+| `EuclideanUnwrapperPETSc` / `SphericalNormalizerPETSc` | PETSc solver binding — replaced by Eigen |
+| `Search` | CoHDS-specific graph search — replaced by CGAL halfedge iteration |
+| `SparseUtility` | Colt sparse matrix utils — replaced by Eigen |

 ---

--- a/doc/roadmap/phases.md
+++ b/doc/roadmap/phases.md
@@ -129,6 +129,89 @@ mesh type.
                 layer, +1 week integration.
 ```

+9d — Cone metrics + sphere utilities  (Java port — 2026 library scan)
+────────────────────────────────────────────────────────────────────
+
+```
+9d.1   ConesUtility  (Java port: unwrapper/ConesUtility.java)
+         → cone_singularities.hpp
+         Fills the "⚠️ data structure only" gap in java-parity.md:
+         - Detect interior cone vertices (angle deficit ≠ 0)
+         - BFS path from cone to mesh boundary → cut edge set
+         - Auto-placement: conjugate gradient on Θ-gradient magnitude
+         - Quantization: snap cone angles to π/2, π/3, π/6 for
+           quad / triangle / hexagonal atlas targets
+         Java reference: unwrapper/ConesUtility.java
+
+9d.2   StereographicUnwrapper + SphereUtility  (Java port)
+         → stereographic_layout.hpp
+         Stereographic projection S²→ℂ∪{∞} + Möbius centering for
+         genus-0 surfaces. Converts spherical DCE output to a flat 2-D atlas.
+         Java reference: unwrapper/StereographicUnwrapper.java (266 lines)
+
+9d.3   MobiusCenteringFunctional  (Java port, optional upgrade)
+         → integrate into layout.hpp normalise_hyperbolic()
+         Variational Möbius centering via Lorentz geometry:
+           E = Σ log(-⟨x,p⟩/√(-⟨x,x⟩))
+         Supplies gradient + Hessian — replaces iterative Fréchet mean.
+         Java reference: functional/MobiusCenteringFunctional.java
+```
+
+9e — Circle pattern layout  (Java port — complement to Phase 9a.1)
+────────────────────────────────────────────────────────────────────
+
+```
+9e     CirclePatternLayout + CirclePatternUtility  (Java port)
+         → circle_pattern_layout.hpp
+         Phase 9a.1 ported the CPEuclidean energy + solver; this phase
+         adds the embedding step (ρ values → actual vertex positions in ℝ²).
+         - CirclePatternUtility: compute per-face radii ρ via NTR solver
+         - CirclePatternLayout: embed from ρ values (intersection-angle model)
+         - CPEuclideanRotation: rotation-invariant CP functional variant
+         Java references: unwrapper/circlepattern/CirclePattern{Layout,Utility}.java
+                          unwrapper/circlepattern/CPEuclideanRotation.java
+```
+
+---
+
+## ◼ New research directions — Phases 10d–10g  (2026 library scan)
+
+These were identified by a full scan of the Java source tree in 2026.
+They extend significantly beyond the Java port into new mathematical territory.
+
+```
+10d   CircleDomainUnwrapper  (Koebe–Andreev–Thurston)
+         → circle_domain_unwrapper.hpp
+         Conformal map of a multiply-connected planar region onto a
+         canonical disk-with-holes (classical complex-analysis result).
+         Java reference: unwrapper/CircleDomainUnwrapper.java (570 lines)
+         Mathematical basis: Koebe–Andreev–Thurston + Beardon–Stephenson 1990
+
+10e   Quasi-isothermic maps
+         → quasiisothermic.hpp
+         Generalisation of conformal maps for meshes where exact conformality
+         is unachievable (high Gaussian curvature, coarse triangulation).
+         Includes: Delaunay pre-conditioning, discrete Beltrami field,
+                   sin-condition functional, Lawson-correspondence parameterization.
+         Java references: unwrapper/quasiisothermic/ (~1 200 lines total)
+         Mathematical basis: Lam 2015 + Bohle–Lam–Pinkall–Reitebuch 2015
+
+10f   Koebe polyhedra
+         → koebe_polyhedron.hpp
+         Koebe–Andreev–Thurston theorem: realize every 3-connected planar
+         graph as a convex polyhedron with edges tangent to the unit sphere.
+         Connects circle packing with 3-D convex geometry.
+         Java reference: unwrapper/koebe/KoebePolyhedron.java (321 lines)
+         Mathematical basis: Koebe 1936 + Thurston 1997 (lecture notes)
+
+10g   Cyclic-symmetry functionals
+         → cyclic_functional.hpp
+         Euclidean and hyperbolic DCE functionals reduced to a cyclic-symmetry
+         quotient — dramatically reduces DOFs for ornamental / symmetric surfaces.
+         Java references: functional/EuclideanCyclicFunctional.java
+                          functional/HyperbolicCyclicFunctional.java (~530 lines)
+```
+
 ---

 ## ◼ Optional / Hypothetical — geometry-central Cross-Comparison
--- a/scripts/doxygen-coverage.sh
+++ b/scripts/doxygen-coverage.sh
@@ -0,0 +1,159 @@
+#!/bin/bash
+# scripts/doxygen-coverage.sh
+#
+# Measure Doxygen documentation coverage of the public C++ API by parsing
+# the Doxygen XML output.  Reports:
+#   * total documentable members (functions, classes, structs, enums,
+#     typedefs, variables) in code/include/**
+#   * how many have a non-empty briefdescription/detaileddescription
+#   * coverage % and list of undocumented members
+#
+# Prerequisite: doxygen must have been run with GENERATE_XML=YES (which
+# the project's Doxyfile sets).  This script invokes it if XML is missing.
+#
+# Usage:
+#   bash scripts/doxygen-coverage.sh               # short summary
+#   bash scripts/doxygen-coverage.sh --list-undoc  # list undocumented members
+#   bash scripts/doxygen-coverage.sh --threshold 95  # fail if coverage < 95 %
+#
+# Exit codes:
+#   0  coverage ≥ threshold (default 0 — informational only)
+#   1  coverage <  threshold
+#   2  XML output missing / could not be parsed
+
+set -eu
+
+XML_DIR="doc/doxygen/xml"
+THRESHOLD=0
+LIST_UNDOC=0
+
+INCLUDE_DETAIL=0
+
+while [ $# -gt 0 ]; do
+  case "$1" in
+    --threshold) THRESHOLD="$2"; shift 2 ;;
+    --list-undoc) LIST_UNDOC=1; shift ;;
+    --include-detail) INCLUDE_DETAIL=1; shift ;;
+    *) echo "Unknown arg: $1" >&2; exit 2 ;;
+  esac
+done
+
+if [ ! -d "$XML_DIR" ]; then
+  echo "XML output missing — running doxygen..."
+  doxygen Doxyfile >/dev/null 2>&1
+fi
+
+if [ ! -d "$XML_DIR" ]; then
+  echo "ERROR: $XML_DIR still missing after doxygen run" >&2
+  exit 2
+fi
+
+python3 - "$XML_DIR" "$LIST_UNDOC" "$THRESHOLD" "$INCLUDE_DETAIL" <<'PYEOF'
+import sys, os, glob, xml.etree.ElementTree as ET
+
+xml_dir, list_undoc, threshold, include_detail = \
+    sys.argv[1], int(sys.argv[2]), float(sys.argv[3]), int(sys.argv[4])
+
+# Implementation-detail namespaces — not part of the public API surface.
+# Skipped by default; pass --include-detail to count them too.
+DETAIL_NAMES = ("::detail::", "::detail_xml::", "::cp_detail::", "::id_detail::",
+                "::detail$", "::detail_xml$", "::cp_detail$", "::id_detail$")
+
+def is_detail(qualified_name: str) -> bool:
+    if include_detail:
+        return False
+    return any(qualified_name.find(d.rstrip("$")) >= 0 for d in DETAIL_NAMES)
+
+# Restrict to compounds whose location is under code/include/ (the
+# public API).  XML output also includes README.md and CLAUDE.md as
+# "file" kind compounds, which we want to skip.
+PUBLIC_PREFIX = os.path.abspath("code/include") + os.sep
+
+KINDS = {"function", "class", "struct", "enum", "typedef", "variable", "namespace"}
+
+total = 0
+documented = 0
+undoc = []
+
+for path in sorted(glob.glob(os.path.join(xml_dir, "*.xml"))):
+    if os.path.basename(path) in {"index.xml", "Doxyfile.xml", "indexpage.xml"}:
+        continue
+    if os.path.basename(path).startswith(("namespacestd", "md_")):
+        continue
+    try:
+        tree = ET.parse(path)
+    except ET.ParseError:
+        continue
+    for cd in tree.iter("compounddef"):
+        kind = cd.attrib.get("kind", "")
+        # only count compounds living in our public include tree
+        loc = cd.find("location")
+        if loc is None:
+            continue
+        file_attr = loc.attrib.get("file", "")
+        if not file_attr.startswith(PUBLIC_PREFIX) and \
+           not file_attr.startswith("code/include/"):
+            continue
+
+        # The compound itself
+        if kind in {"class", "struct", "namespace"}:
+            cname = cd.findtext("compoundname", "?")
+            if not is_detail(cname):
+                total += 1
+                brief = cd.find("briefdescription")
+                detail = cd.find("detaileddescription")
+                has_doc = (brief is not None and len("".join(brief.itertext()).strip()) > 0) or \
+                          (detail is not None and len("".join(detail.itertext()).strip()) > 0)
+                if has_doc:
+                    documented += 1
+                else:
+                    undoc.append(f"{kind:9s}  {cname}  ({file_attr}:{loc.attrib.get('line','?')})")
+
+        # Members inside the compound
+        for memberdef in cd.iter("memberdef"):
+            mkind = memberdef.attrib.get("kind", "")
+            if mkind not in KINDS:
+                continue
+            prot = memberdef.attrib.get("prot", "public")
+            if prot != "public":
+                continue
+            name = memberdef.findtext("name", "?")
+            qual = memberdef.findtext("qualifiedname", name)
+            if is_detail(qual):
+                continue
+            total += 1
+            brief = memberdef.find("briefdescription")
+            detail = memberdef.find("detaileddescription")
+            has_doc = (brief is not None and len("".join(brief.itertext()).strip()) > 0) or \
+                      (detail is not None and len("".join(detail.itertext()).strip()) > 0)
+            if has_doc:
+                documented += 1
+            else:
+                mloc = memberdef.find("location")
+                fl = mloc.attrib.get("file", "?") if mloc is not None else "?"
+                ln = mloc.attrib.get("line", "?") if mloc is not None else "?"
+                undoc.append(f"{mkind:9s}  {qual}  ({fl}:{ln})")
+
+if total == 0:
+    print("ERROR: no public members found — check that GENERATE_XML=YES and EXTRACT_ALL=YES")
+    sys.exit(2)
+
+pct = 100.0 * documented / total
+print(f"Doxygen coverage (public symbols under code/include/):")
+print(f"  documented:   {documented}")
+print(f"  total:        {total}")
+print(f"  coverage:     {pct:.1f}%")
+print(f"  undocumented: {total - documented}")
+
+if list_undoc:
+    print()
+    print("Undocumented symbols:")
+    for s in undoc:
+        print(f"  {s}")
+
+if pct < threshold:
+    print(f"\nFAIL: coverage {pct:.1f}% < threshold {threshold}%", file=sys.stderr)
+    sys.exit(1)
+
+sys.exit(0)
+PYEOF
--- a/scripts/gen-headers-md.py
+++ b/scripts/gen-headers-md.py
@@ -0,0 +1,242 @@
+#!/usr/bin/env python3
+"""scripts/gen-headers-md.py
+
+Auto-generate `doc/api/headers.md` from the Doxygen XML index.  For
+each public header under `code/include/`, the script extracts:
+
+  * the header's `@file` brief description (if present), else falls
+    back to the first leading `//` comment block in the source file;
+  * the names of all public classes, structs, free functions, enums
+    and typedefs declared at file scope (one-line list).
+
+Headers are grouped by their directory:
+
+  * `code/include/CGAL/`        →  "CGAL public API"
+  * `code/include/CGAL/.../`    →  nested CGAL subgroups (internal_np, …)
+  * `code/include/*.hpp`        →  "Core (conformallab namespace)"
+
+The file is auto-generated; do NOT hand-edit.  Regenerate via:
+
+    bash scripts/regen-docs.sh          # convenience wrapper
+    # or:
+    doxygen Doxyfile && python3 scripts/gen-headers-md.py
+
+Prerequisite: `GENERATE_XML = YES` in the Doxyfile (already set).
+"""
+
+from __future__ import annotations
+import os, sys, glob, xml.etree.ElementTree as ET
+from collections import defaultdict
+
+REPO_ROOT = os.path.abspath(os.path.dirname(os.path.dirname(__file__)))
+XML_DIR   = os.path.join(REPO_ROOT, "doc", "doxygen", "xml")
+OUT_PATH  = os.path.join(REPO_ROOT, "doc", "api", "headers.md")
+PUBLIC_ABS = os.path.join(REPO_ROOT, "code", "include") + os.sep
+PUBLIC_REL = "code/include/"
+
+
+def is_public(path: str) -> bool:
+    return path.startswith(PUBLIC_ABS) or path.startswith(PUBLIC_REL)
+
+
+def relpath(path: str) -> str:
+    if path.startswith(PUBLIC_ABS):
+        return path[len(PUBLIC_ABS):]
+    if path.startswith(PUBLIC_REL):
+        return path[len(PUBLIC_REL):]
+    return path
+
+
+def text_of(elem: ET.Element | None) -> str:
+    """Concatenate all text under elem, trimming whitespace."""
+    if elem is None:
+        return ""
+    return " ".join("".join(elem.itertext()).split()).strip()
+
+
+def first_sentence(text: str, max_chars: int = 240) -> str:
+    """Return the first sentence (or first `max_chars` chars) of text.
+    Heuristic: split at the first ". " (period + space) that follows
+    at least 30 characters, else hard-truncate."""
+    text = text.strip()
+    if not text:
+        return ""
+    # Strip Doxygen tag-class artefacts that bleed into briefs
+    text = text.replace("CGAL::Discrete_conformal_map", "")  # spurious self-reference
+    for cut in [". ", ".\n", "; "]:
+        i = text.find(cut, 30)
+        if 0 < i < max_chars:
+            return text[: i + 1].strip()
+    if len(text) > max_chars:
+        return text[: max_chars - 1].rstrip() + "…"
+    return text
+
+
+def first_leading_comment(path: str) -> str:
+    """Fallback: extract the first contiguous //-comment block at top of file."""
+    try:
+        with open(path, encoding="utf-8") as f:
+            lines = []
+            started = False
+            for ln in f:
+                s = ln.strip()
+                if s.startswith("//"):
+                    lines.append(s.lstrip("/").strip())
+                    started = True
+                elif started:
+                    break
+                elif not s or s.startswith("#pragma"):
+                    continue
+                else:
+                    break
+            return " ".join(l for l in lines if l)
+    except OSError:
+        return ""
+
+
+def scan_xml() -> dict[str, dict]:
+    """Walk all compounddef kind='file' XMLs and return a dict
+    keyed by absolute path → {'brief': str, 'symbols': [list]}.
+    """
+    out: dict[str, dict] = {}
+    for xml in sorted(glob.glob(os.path.join(XML_DIR, "*.xml"))):
+        bn = os.path.basename(xml)
+        if bn in {"index.xml", "Doxyfile.xml", "indexpage.xml"}:
+            continue
+        try:
+            tree = ET.parse(xml)
+        except ET.ParseError:
+            continue
+        for cd in tree.iter("compounddef"):
+            if cd.attrib.get("kind") != "file":
+                continue
+            loc = cd.find("location")
+            if loc is None:
+                continue
+            path = loc.attrib.get("file", "")
+            if not is_public(path):
+                continue
+            # Normalise to absolute for downstream file IO
+            if path.startswith(PUBLIC_REL):
+                path = os.path.join(REPO_ROOT, path)
+            ext = os.path.splitext(path)[1]
+            if ext not in {".h", ".hpp"}:
+                continue
+            brief = text_of(cd.find("briefdescription"))
+            if not brief:
+                brief = text_of(cd.find("detaileddescription"))
+            if not brief:
+                brief = first_leading_comment(path)
+            brief = first_sentence(brief)
+
+            symbols: list[str] = []
+            # Free symbols defined directly in the file
+            for sect in cd.iter("sectiondef"):
+                for m in sect.iter("memberdef"):
+                    mk = m.attrib.get("kind", "")
+                    if mk not in {"function", "typedef", "enum", "variable"}:
+                        continue
+                    if m.attrib.get("prot", "public") != "public":
+                        continue
+                    name = m.findtext("name", "?")
+                    qual = m.findtext("qualifiedname", name)
+                    if "::detail" in qual or "::cp_detail" in qual or \
+                       "::id_detail" in qual or "::detail_xml" in qual:
+                        continue
+                    if mk == "function":
+                        symbols.append(f"`{name}()`")
+                    elif mk == "typedef":
+                        symbols.append(f"`{name}`")
+                    elif mk == "enum":
+                        symbols.append(f"`enum {name}`")
+                    else:
+                        symbols.append(f"`{name}`")
+            # Classes/structs declared in the file (skip template
+            # specialisations whose name contains angle brackets — those
+            # appear as duplicates of the primary template).
+            for cref in cd.iter("innerclass"):
+                cname = cref.text or ""
+                if "detail" in cname or "<" in cname:
+                    continue
+                short = cname.rsplit("::", 1)[-1]
+                symbols.append(f"`{short}`")
+            # De-duplicate, preserve order
+            seen = set()
+            uniq = []
+            for s in symbols:
+                if s in seen:
+                    continue
+                seen.add(s); uniq.append(s)
+            out[path] = {"brief": brief, "symbols": uniq}
+    return out
+
+
+def group_key(path: str) -> tuple[str, int]:
+    rel = relpath(path)
+    parts = rel.split(os.sep)
+    if parts[0] == "CGAL":
+        if len(parts) == 2:
+            return ("CGAL public API (`<CGAL/...>`)", 1)
+        return (f"CGAL internals (`<CGAL/{parts[1]}/...>`)", 2)
+    return ("Core (`conformallab` namespace, `<...>`)", 3)
+
+
+def render(headers: dict[str, dict]) -> str:
+    groups: dict[tuple[str, int], list[tuple[str, dict]]] = defaultdict(list)
+    for path, info in headers.items():
+        groups[group_key(path)].append((path, info))
+
+    lines: list[str] = []
+    lines.append("<!-- AUTO-GENERATED by scripts/gen-headers-md.py — do not edit by hand. -->")
+    lines.append("<!-- Source of truth: the `\\file` / leading-comment briefs in code/include/. -->")
+    lines.append("")
+    lines.append("# Public Headers (`code/include/`)")
+    lines.append("")
+    lines.append("All algorithms are header-only.  Include the headers you need")
+    lines.append("directly — there is no compiled library to link against (only")
+    lines.append("GTest for tests and CGAL/Eigen for the CGAL-dependent headers).")
+    lines.append("")
+    lines.append("This page is **regenerated** from Doxygen XML on every push to")
+    lines.append("`main` that touches a public header (see")
+    lines.append("`.gitea/workflows/doxygen-pages.yml` and `scripts/gen-headers-md.py`).")
+    lines.append("To improve a description, edit the `\\file` brief at the top of")
+    lines.append("the corresponding header, then re-run `bash scripts/regen-docs.sh`.")
+    lines.append("")
+
+    for key in sorted(groups.keys(), key=lambda k: (k[1], k[0])):
+        title, _ = key
+        lines.append(f"## {title}")
+        lines.append("")
+        lines.append("| Header | Brief | Public symbols |")
+        lines.append("|---|---|---|")
+        for path, info in sorted(groups[key], key=lambda x: x[0]):
+            rel = relpath(path).replace(os.sep, "/")
+            brief = info["brief"] or "_(undocumented — add a `\\file` brief at the top of the header)_"
+            # Defensive: collapse pipe/newline chars that would break the table
+            brief = brief.replace("|", "\\|").replace("\n", " ").strip()
+            syms  = ", ".join(info["symbols"][:10]) or "_(no public symbols)_"
+            if len(info["symbols"]) > 10:
+                syms += f", … (+{len(info['symbols']) - 10} more)"
+            lines.append(f"| `{rel}` | {brief} | {syms} |")
+        lines.append("")
+
+    return "\n".join(lines) + "\n"
+
+
+def main() -> int:
+    if not os.path.isdir(XML_DIR):
+        print(f"ERROR: {XML_DIR} missing — run `doxygen Doxyfile` first", file=sys.stderr)
+        return 2
+    headers = scan_xml()
+    if not headers:
+        print("ERROR: no public headers found in XML output", file=sys.stderr)
+        return 2
+    text = render(headers)
+    with open(OUT_PATH, "w", encoding="utf-8") as f:
+        f.write(text)
+    print(f"Wrote {OUT_PATH} ({len(headers)} headers, {sum(len(h['symbols']) for h in headers.values())} symbols)")
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())
--- a/scripts/regen-docs.sh
+++ b/scripts/regen-docs.sh
@@ -0,0 +1,15 @@
+#!/bin/sh
+# scripts/regen-docs.sh — convenience wrapper:
+#   1. run doxygen to refresh XML + HTML output
+#   2. regenerate doc/api/headers.md from the XML
+#   3. report coverage so the operator sees the impact
+#
+# Intended for local use ("I just added a \file brief — regenerate
+# the markdown landing page"); the same steps run automatically in
+# .gitea/workflows/doxygen-pages.yml on every push to main that touches
+# the public headers.
+set -eu
+cd "$(dirname "$0")/.."
+doxygen Doxyfile >/dev/null
+python3 scripts/gen-headers-md.py
+bash scripts/doxygen-coverage.sh