docs: Doxygen-API + Validierungsprotokoll + Porting-Tutorial

Für einen Mathematiker der unabhängig validieren und eigene Forschung
einbringen möchte.

Doxygen-Kommentare (code/include/):
  newton_solver.hpp — newton_euclidean(), newton_spherical(), newton_hyper_ideal()
    je mit \param, \return, \note, \see inkl. mathematischer Begründung
    (Konvexität, Vorzeichenkonvention, SparseQR-Fallback-Erklärung)
  layout.hpp — euclidean_layout(), spherical_layout(), hyper_ideal_layout()
    mit vollständiger Parameter-Doku, halfedge_uv-Semantik, Poincaré-Disk-Note

Neues Dokument:
  doc/math/validation-protocol.md
    7 reproduzierbare Checks mit konkreten Befehlen und erwartetem Output:
    0. 170 Tests, 1 Skip
    1. Gauss–Bonnet exakt (1e-10)
    2. FD-Gradientencheck < 1e-6 für alle 3 Geometrien
    3. Newton-Konvergenz < 50 Iterationen
    4. τ ∈ SL(2,ℤ)-Fundamentaldomäne (3 Invarianten)
    5. Möbius-Arithmetik (Inverse, Compose, from_three)
    6. End-to-End-Pipeline
    7. Manueller τ-Check für torus_4x4.off (Codebeispiel)

Neues Tutorial:
  doc/tutorials/add-inversive-distance.md
    Vollständiger Step-by-Step-Port von Phase 9a (Luo 2004):
    Header anlegen, Energie/Gradient implementieren, FD-Check,
    Newton-Wrapper, CMakeLists, Java-Referenzvergleich, Checkliste.

doc/getting-started.md:
  Abschnitt "Known issues": macOS-Finder-Duplikate (rm-Befehl),
  Warnung "First build 30–90s" (Tarball-Extraktion)

README.md:
  Zwei neue Links in der Dokumentationstabelle (validation-protocol,
  tutorial)

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

README.md

@@ -85,6 +85,8 @@ Layout2D layout = euclidean_layout(mesh, res.x, maps);
 | **Project structure** — directory tree + build targets | [doc/architecture/project-structure.md](doc/architecture/project-structure.md) |
 | **Discrete conformal theory** — mathematical background for collaborators | [doc/math/discrete-conformal-theory.md](doc/math/discrete-conformal-theory.md) |
 | **Validation** — known analytic results + how to verify them | [doc/math/validation.md](doc/math/validation.md) |
+| **Validation protocol** — concrete commands with expected outputs | [doc/math/validation-protocol.md](doc/math/validation-protocol.md) |
+| **Tutorial: add a new functional** — step-by-step Inversive-Distance port | [doc/tutorials/add-inversive-distance.md](doc/tutorials/add-inversive-distance.md) |
 | **Geometry modes** — Euclidean / Spherical / HyperIdeal comparison | [doc/math/geometry-modes.md](doc/math/geometry-modes.md) |
 | **References** — all papers by module | [doc/math/references.md](doc/math/references.md) |
 | **Roadmap** — Phases 1–10 | [doc/roadmap/phases.md](doc/roadmap/phases.md) |

code/include/layout.hpp

@@ -452,6 +452,29 @@ inline void set_root_huv_2d(
 } // namespace detail
 
 // ── Euclidean layout ──────────────────────────────────────────────────────────
+
+/// Embed the mesh in ℝ² using the Euclidean metric encoded in x.
+///
+/// Runs priority-BFS trilateration: places vertices in order of BFS depth
+/// from the root face (largest 3D area), so errors accumulate last.
+/// For closed genus-g surfaces a CutGraph must be supplied — otherwise
+/// the layout will have a seam discontinuity (Layout2D::has_seam = true).
+///
+/// \param mesh      Input surface mesh.  Must have lambda0 and v_idx set in maps.
+/// \param x         DOF vector returned by newton_euclidean().
+/// \param maps      EuclideanMaps (lambda0, v_idx, e_idx).
+/// \param cut       Optional cut graph (compute_cut_graph()).  Pass nullptr for
+///                  open meshes or if seams are acceptable.
+/// \param holonomy  If non-null and cut != nullptr, receives the lattice
+///                  translations ω_i ∈ ℂ per cut edge.
+///                  Pass to compute_period_matrix() for the conformal modulus τ.
+/// \param normalise If true, calls normalise_euclidean() on the result:
+///                  centroid → origin, major axis → x-axis (PCA).
+/// \return          Layout2D with .uv[v] (per-vertex UV) and
+///                  .halfedge_uv[h] (per-halfedge UV for texture atlasing).
+///
+/// \note halfedge_uv differs from uv at seam edges: the two sides of a cut
+///       carry different UV coordinates for proper GPU texture atlasing.
 inline Layout2D euclidean_layout(
     ConformalMesh&             mesh,
     const std::vector<double>& x,
@@ -571,6 +594,20 @@ inline Layout2D euclidean_layout(
 }
 
 // ── Spherical layout ──────────────────────────────────────────────────────────
+
+/// Embed the mesh on the unit sphere S² using the spherical metric encoded in x.
+///
+/// Runs priority-BFS trilateration using the spherical law of cosines.
+/// Typical use: genus-0 (sphere-like) surfaces after newton_spherical().
+///
+/// \param mesh      Input genus-0 surface mesh.
+/// \param x         DOF vector returned by newton_spherical().
+/// \param maps      SphericalMaps.
+/// \param cut       Optional cut graph (rarely needed for genus-0).
+/// \param holonomy  If non-null, receives rotational holonomies (spherical).
+/// \param normalise If true, calls normalise_spherical(): rotates the centroid
+///                  to the north pole (Rodrigues rotation formula).
+/// \return          Layout3D with .xyz[v] ∈ S² ⊂ ℝ³ for each vertex.
 inline Layout3D spherical_layout(
     ConformalMesh&             mesh,
     const std::vector<double>& x,
@@ -670,6 +707,29 @@ inline Layout3D spherical_layout(
 }
 
 // ── HyperIdeal layout (Poincaré disk) — exact trilateration ──────────────────
+
+/// Embed the mesh in the Poincaré disk (H²) using the hyperbolic metric encoded in x.
+///
+/// Runs priority-BFS trilateration using exact Möbius-isometric placement:
+/// each new vertex is located by solving the hyperbolic law of cosines and
+/// applying a Möbius map to position it in the disk.
+/// For closed genus-g surfaces (g ≥ 1) a CutGraph is required.
+///
+/// \param mesh      Input genus-g surface mesh (g ≥ 1).
+/// \param x         DOF vector returned by newton_hyper_ideal()
+///                  (vertex b_v and edge a_e variables).
+/// \param maps      HyperIdealMaps (lambda0, v_idx, e_idx).
+/// \param cut       CutGraph from compute_cut_graph().  Required for closed surfaces.
+/// \param holonomy  If non-null, receives the Möbius maps T_i ∈ SU(1,1) per cut edge.
+///                  Pass to compute_period_matrix() for holonomy analysis.
+/// \param normalise If true, calls normalise_hyperbolic(): iterative face-area-weighted
+///                  Möbius centring (Fréchet mean, 30 iterations) → disk origin.
+/// \return          Layout2D with .uv[v] ∈ Poincaré disk (|uv| < 1) for each vertex,
+///                  and .halfedge_uv[h] for seam-aware texture atlasing.
+///
+/// \note All vertex positions satisfy |uv[v]| < 1 (inside the Poincaré disk)
+///       if the metric is hyperbolic.  Points on or outside the boundary indicate
+///       a non-hyperbolic metric (Gauss–Bonnet violation).
 inline Layout2D hyper_ideal_layout(
     ConformalMesh&             mesh,
     const std::vector<double>& x,

code/include/newton_solver.hpp

@@ -139,9 +139,29 @@ inline std::vector<double> line_search(
 } // namespace detail
 
 // ── Euclidean Newton solver ────────────────────────────────────────────────────
-//
+
-// Minimises the Euclidean discrete conformal energy by solving G(x) = 0.
+/// Solve the Euclidean discrete conformal problem: find u ∈ ℝ^V such that
-// The Hessian H is PSD; Eigen::SimplicialLDLT is used directly.
+/// Σ_{faces adj v} α_v(u) = Θ_v for all vertices v.
+///
+/// Starting from x0, Newton's method minimises E(u) (the Euclidean DCE energy,
+/// which is convex) by iterating  u ← u − H⁻¹·G  with backtracking line search.
+/// The Hessian H is the cotangent Laplacian — PSD with one zero eigenvalue on
+/// closed surfaces (gauge mode).  A SparseQR fallback handles this automatically.
+///
+/// \param mesh      Input triangulated surface (edges must carry lambda0 + theta_v).
+/// \param x0        Initial DOF vector (length = number of free vertices).
+///                  Pass all-zeros for a flat start (typical).
+/// \param m         EuclideanMaps: lambda0[e], theta_v[v], v_idx[v] must be set.
+///                  Call setup_euclidean_maps() + compute_euclidean_lambda0_from_mesh()
+///                  + enforce_gauss_bonnet() before passing here.
+/// \param tol       Convergence threshold on max |G_i|.  Default: 1e-8.
+/// \param max_iter  Maximum Newton iterations.  Default: 200.
+/// \return          NewtonResult{x*, iterations, grad_inf_norm, converged}.
+///
+/// \note On closed meshes without a pinned vertex, SimplicialLDLT detects the
+///       gauge singularity and falls back to SparseQR automatically.
+///
+/// \see doc/math/discrete-conformal-theory.md §3 for the mathematical background.
 inline NewtonResult newton_euclidean(
     ConformalMesh&             mesh,
     std::vector<double>        x0,
@@ -199,10 +219,28 @@ inline NewtonResult newton_euclidean(
 }
 
 // ── Spherical Newton solver ───────────────────────────────────────────────────
-//
+
-// Solves G(x) = 0 for the spherical discrete conformal functional.
+/// Solve the spherical discrete conformal problem: find u ∈ ℝ^V such that
-// The Hessian H is NSD at the solution; −H is PSD, so we factorise −H and
+/// Σ_{faces adj v} α_v(u) = Θ_v for all vertices v (genus-0 / sphere-like surfaces).
-// solve (−H)·Δx = G  ⟺  H·Δx = −G.
+///
+/// The spherical DCE energy is *concave*, so the Hessian H is NSD at the solution.
+/// The solver factorises −H (which is PSD) and solves (−H)·Δx = G.
+/// A gauge vertex must be pinned (set v_idx = -1) to remove the rotational mode.
+///
+/// \param mesh      Input triangulated surface, genus 0.
+/// \param x0        Initial DOF vector (length = free vertices, excluding gauge_vertex).
+///                  All-zeros is a good start.
+/// \param m         SphericalMaps: lambda0[e], theta_v[v], v_idx[v], gauge_vertex set.
+///                  Call setup_spherical_maps() + compute_spherical_lambda0_from_mesh()
+///                  + enforce_gauss_bonnet() (checks Σ(2π-Θ) > 0) before passing here.
+/// \param tol       Convergence threshold on max |G_i|.  Default: 1e-8.
+/// \param max_iter  Maximum Newton iterations.  Default: 200.
+/// \return          NewtonResult{x*, iterations, grad_inf_norm, converged}.
+///
+/// \note Unlike the Euclidean solver, the spherical solver does NOT need a SparseQR
+///       fallback — the gauge vertex pins the null mode directly.
+///
+/// \see doc/math/geometry-modes.md §Spherical for sign-convention details.
 inline NewtonResult newton_spherical(
     ConformalMesh&             mesh,
     std::vector<double>        x0,
@@ -258,17 +296,31 @@ inline NewtonResult newton_spherical(
 }
 
 // ── HyperIdeal Newton solver ──────────────────────────────────────────────────
-//
+
-// Solves G(x) = 0 for the hyper-ideal discrete conformal functional.
+/// Solve the hyper-ideal discrete conformal problem: find (b, a) ∈ ℝ^{V+E} such that
-//
+/// Σ β_v(b,a) = Θ_v  and  Σ α_e(b,a) = θ_e  for all vertices v and edges e.
-// Gradient sign convention (opposite to Euclidean/Spherical):
+///
-//   G_v = Σ β_v − Θ_v,   G_e = Σ α_e − θ_e   (actual − target)
+/// Used for genus-g surfaces (g ≥ 1) under hyperbolic cone metrics.
-//
+/// The energy is *strictly convex* (Springborn 2020, Theorem 1.3), so Newton
-// The hyper-ideal energy is strictly convex (Springborn 2020), so H is PSD
+/// converges globally from any starting point.
-// and SimplicialLDLT (with SparseQR fallback) applies directly.
+///
-//
+/// DOF layout: first V_free entries are vertex variables b_v (hyper-ideal radii),
-// The Hessian is computed by symmetric finite differences of G (see
+/// followed by E entries for edge variables a_e (intersection angles).
-// hyper_ideal_hessian.hpp); replace with an analytical Hessian in Phase 5.
+/// Use assign_all_dof_indices(mesh, maps) to set v_idx and e_idx automatically —
+/// no vertex needs to be pinned.
+///
+/// \param mesh      Input triangulated surface, genus g ≥ 1.
+/// \param x0        Initial DOF vector (length = V + E).  All-zeros typical.
+/// \param m         HyperIdealMaps: lambda0[e], theta_v[v], v_idx[v], e_idx[e] set.
+///                  Call setup_hyper_ideal_maps() + compute_hyper_ideal_lambda0_from_mesh().
+/// \param tol       Convergence threshold on max |G_i|.  Default: 1e-8.
+/// \param max_iter  Maximum Newton iterations.  Default: 200.
+/// \param hess_eps  Finite-difference step for Hessian approximation.  Default: 1e-5.
+///                  (Phase 9b will replace this with an analytic Hessian.)
+/// \return          NewtonResult{x*, iterations, grad_inf_norm, converged}.
+///
+/// \see Springborn (2020), Theorem 1.3 for the strict convexity proof.
+/// \see doc/math/geometry-modes.md §Hyper-ideal for DOF layout details.
 inline NewtonResult newton_hyper_ideal(
     ConformalMesh&             mesh,
     std::vector<double>        x0,

doc/getting-started.md

@@ -127,6 +127,30 @@ After a full build (`-DWITH_CGAL=ON`):
 
 ---
 
+## Known issues
+
+### macOS Finder duplicates
+
+macOS Finder sometimes creates duplicate files named `foo 2.hpp` when copying
+the repository.  These cause confusing compile errors ("redefinition of …").
+
+**Fix (run once from the repo root):**
+```bash
+find code/include -name "* 2.*" -delete
+find code/include -name "*\ 2.*" -delete
+```
+
+Files without a ` 2` suffix are always canonical — the duplicates are safe to delete.
+
+### First build is slow
+
+The first CMake configure extracts four tarballs (Eigen 3.4, CGAL 6.1.1,
+libigl 2.6, GLFW 3.4) and downloads GTest via `FetchContent`.
+Allow **30–90 seconds** for the first configure.  Subsequent builds are fast
+(< 10 s incremental).
+
+---
+
 ## Rebuilding the CI Docker image
 
 The CI runner is a self-hosted Raspberry Pi (ARM64). After changes to

doc/math/validation-protocol.md

@@ -0,0 +1,213 @@
+# Validation Protocol
+
+Concrete, reproducible steps to verify the mathematical correctness of
+conformallab++.  Every check below has a deterministic expected outcome.
+
+---
+
+## Prerequisites
+
+```bash
+cmake -S code -B build -DWITH_CGAL=ON -DCMAKE_BUILD_TYPE=Release
+cmake --build build --target conformallab_cgal_tests -j$(nproc)
+```
+
+---
+
+## Check 0 — All tests pass
+
+```bash
+ctest --test-dir build -R "^cgal\." --output-on-failure
+```
+
+**Expected output (last lines):**
+```
+100% tests passed, 0 tests failed out of 170
+
+The following tests did not run:
+    206 - cgal.HomologyGenerators.Genus2_FourGeneratorPaths_BLOCKED (Skipped)
+```
+
+If any test fails, stop — the implementation is broken.
+
+---
+
+## Check 1 — Gauss–Bonnet (topological identity, error < 1e-10)
+
+```bash
+./build/conformallab_cgal_tests --gtest_filter="GaussBonnet.*" -v
+```
+
+**Expected:** all 8 tests `[  PASSED  ]`
+
+What is verified: for each test mesh (tetrahedron χ=2, torus χ=0, open mesh χ=1):
+
+```
+Σᵥ (2π − Θᵥ)   =   2π · χ(M)   ± 1e-10
+```
+
+This is a **pure topology check** — it fails only if vertex/face counts or
+property-map assignments are wrong.  It does not depend on the Newton solver.
+
+---
+
+## Check 2 — Euclidean gradient consistency (FD vs. analytic, error < 1e-6)
+
+```bash
+./build/conformallab_cgal_tests --gtest_filter="EuclideanFunctional.GradientCheck*" -v
+```
+
+**Expected:** all `GradientCheck_*` tests `[  PASSED  ]`
+
+What is verified: for ε = 1e-5,
+
+```
+|G(u)ᵢ  −  (E(u + εeᵢ) − E(u − εeᵢ)) / (2ε)| < 1e-6
+```
+
+This check **proves that the energy and its gradient are mathematically consistent**.
+A failing FD check means Newton will converge to the wrong point — it is the most
+important correctness check for any new functional.
+
+Also run for Spherical and HyperIdeal:
+```bash
+./build/conformallab_cgal_tests --gtest_filter="SphericalFunctional.GradientCheck*" -v
+./build/conformallab_cgal_tests --gtest_filter="HyperIdealFunctional.GradientCheck*" -v
+```
+
+---
+
+## Check 3 — Newton convergence on canonical test meshes
+
+```bash
+./build/conformallab_cgal_tests --gtest_filter="NewtonSolver.*" -v
+```
+
+**Expected:** all 11 tests `[  PASSED  ]`
+
+Each test verifies:
+- `res.converged == true`
+- `res.grad_inf_norm < 1e-8`
+- `res.iterations < 50`  (typically 5–20 for the small test meshes)
+
+---
+
+## Check 4 — Period matrix: SL(2,ℤ)-reduction invariants
+
+```bash
+./build/conformallab_cgal_tests --gtest_filter="PeriodMatrix.*" -v
+```
+
+**Expected:** all 7 tests `[  PASSED  ]`
+
+The three mathematical invariants checked for **any** genus-1 output τ:
+
+| Property | Condition | Why |
+|---|---|---|
+| Upper half-plane | `Im(τ) > 0` | τ encodes a positive-area lattice |
+| Outside unit disk | `|τ| ≥ 1 − 1e-10` | SL(2,ℤ) reduction step |
+| Vertical strip | `|Re(τ)| ≤ 0.5 + 1e-10` | SL(2,ℤ) reduction step |
+
+These hold for **any** well-formed genus-1 mesh — they are topology, not geometry.
+
+---
+
+## Check 5 — Möbius arithmetic (complex analysis correctness)
+
+```bash
+./build/conformallab_cgal_tests --gtest_filter="MobiusMap.*" -v
+```
+
+**Expected:** all 8 tests `[  PASSED  ]`
+
+What is verified:
+- `T ∘ T⁻¹ = Id` (inverse is correct)
+- `(T₁ ∘ T₂)(z) = T₁(T₂(z))` (composition is associative)
+- `from_three(z₁, z₂, z₃)` maps z₁→0, z₂→1, z₃→∞ (unique Möbius transformation)
+
+A bug here would corrupt **all** hyperbolic holonomy computation.
+
+---
+
+## Check 6 — End-to-end pipeline (build + solve + layout)
+
+```bash
+./build/conformallab_cgal_tests --gtest_filter="Pipeline.*" -v
+```
+
+**Expected:** all 5 tests `[  PASSED  ]`
+
+What is verified: starting from a mesh file, the full pipeline
+(setup → Gauss-Bonnet → Newton → layout → serialise → reload)
+produces a consistent result.
+
+---
+
+## Check 7 — Manual torus τ verification
+
+This check requires adding a small program (or modifying an existing test).
+It validates that `torus_4x4.off` produces τ in the fundamental domain:
+
+```cpp
+#include "conformal_mesh.hpp"
+#include "euclidean_functional.hpp"
+#include "newton_solver.hpp"
+#include "cut_graph.hpp"
+#include "layout.hpp"
+#include "period_matrix.hpp"
+#include "mesh_io.hpp"
+#include <iostream>
+
+int main() {
+    conformallab::ConformalMesh mesh;
+    conformallab::load_mesh(mesh, "code/data/off/torus_4x4.off");
+
+    auto maps = conformallab::setup_euclidean_maps(mesh);
+    conformallab::compute_euclidean_lambda0_from_mesh(mesh, maps);
+    conformallab::enforce_gauss_bonnet(mesh, maps);
+
+    auto res = conformallab::newton_euclidean(mesh, std::vector<double>(maps.n_dof, 0.0), maps);
+    std::cout << "Converged: " << res.converged
+              << "  iterations: " << res.iterations
+              << "  |G|∞: " << res.grad_inf_norm << "\n";
+
+    auto cg  = conformallab::compute_cut_graph(mesh);
+    conformallab::HolonomyData hol;
+    conformallab::euclidean_layout(mesh, res.x, maps, &cg, &hol, true);
+    auto pd  = conformallab::compute_period_matrix(hol);
+
+    std::cout << "τ = " << pd.tau_reduced.real()
+              << " + " << pd.tau_reduced.imag() << "i\n";
+    std::cout << "|τ| = " << std::abs(pd.tau_reduced) << "\n";
+    std::cout << "|Re(τ)| = " << std::abs(pd.tau_reduced.real()) << "\n";
+}
+```
+
+**Expected output (torus_4x4.off, R=2, r=1 torus of revolution):**
+```
+Converged: 1   iterations: <30   |G|∞: <1e-8
+τ = [small]  + [positive]i        (Re close to 0 by 4-fold symmetry)
+|τ| ≥ 1.0                          (fundamental domain)
+|Re(τ)| ≤ 0.5                      (fundamental domain)
+```
+
+The exact value of Im(τ) depends on the 3D embedding (R=2, r=1 gives unequal
+inner/outer edge lengths).  Use `torus_8x8.off` for a finer approximation.
+
+---
+
+## Summary checklist
+
+```
+[ ] Check 0: 170 tests pass, 1 skip
+[ ] Check 1: Gauss–Bonnet exact (1e-10)
+[ ] Check 2: FD gradient < 1e-6 for all 3 geometries
+[ ] Check 3: Newton convergence < 50 iterations
+[ ] Check 4: τ in SL(2,ℤ) fundamental domain
+[ ] Check 5: Möbius arithmetic (inverse, compose, from_three)
+[ ] Check 6: End-to-end pipeline
+[ ] Check 7: Torus τ in upper half-plane with correct symmetry
+```
+
+All checks are deterministic and do not depend on random initialization or
+floating-point non-determinism beyond standard IEEE-754.

doc/tutorials/add-inversive-distance.md

@@ -0,0 +1,202 @@
+# Tutorial: Porting the Inversive-Distance Functional (Phase 9a)
+
+This is a complete, step-by-step example of how to add a new functional
+to conformallab++.  It ports `InversiveDistanceFunctional.java` from the
+Java reference implementation (Luo 2004).
+
+**Prerequisite:** Read [doc/api/extending.md](../api/extending.md) first
+for the general pattern.  This tutorial fills in every detail for one
+specific case.
+
+---
+
+## Mathematical background
+
+Inversive distance (Luo 2004) uses a different edge-length update formula:
+
+```
+Given inversive distances I_{ij} ∈ ℝ for each edge {i,j}:
+
+  cosh(l̃_{ij}) = I_{ij} · cosh((u_i + u_j) / 2)
+               + (cosh²((u_i - u_j) / 2) - 1) · ...
+```
+
+For the Euclidean version the angle formula is the same law of cosines,
+but the log-lengths Λ̃ are computed differently from the u-vector.
+
+**Reference:** Luo, F. (2004). *Combinatorial Yamabe Flow on Surfaces*.
+Communications in Contemporary Mathematics, 6(5), 765–780.
+
+**Java source:** `de.varylab.discreteconformal.functional.InversiveDistanceFunctional`
+
+---
+
+## Step 1 — Create the header
+
+```bash
+cp code/include/euclidean_functional.hpp \
+   code/include/inversive_distance_functional.hpp
+```
+
+Edit the new file:
+
+```cpp
+#pragma once
+// inversive_distance_functional.hpp
+//
+// Phase 9a — Inversive-distance discrete conformal functional (Luo 2004).
+// Ported from de.varylab.discreteconformal.functional.InversiveDistanceFunctional.
+//
+// Usage: identical to euclidean_functional.hpp — replace lambda0 with
+// inversive_distance0 (the initial inversive distances per edge).
+
+#include "conformal_mesh.hpp"
+#include <Eigen/Dense>
+#include <cmath>
+
+namespace conformallab {
+
+struct InversiveDistanceMaps {
+    CGAL::Surface_mesh<CGAL::Simple_cartesian<double>::Point_3>
+        ::Property_map<Edge_index, double>   inv_dist0; ///< I_{ij} per edge
+    CGAL::Surface_mesh<CGAL::Simple_cartesian<double>::Point_3>
+        ::Property_map<Vertex_index, double> theta_v;   ///< target corner angles Θ_v
+    CGAL::Surface_mesh<CGAL::Simple_cartesian<double>::Point_3>
+        ::Property_map<Vertex_index, int>    v_idx;     ///< DOF index (-1 = pinned)
+};
+
+// ... (follow euclidean_functional.hpp exactly, replacing lambda0 with inv_dist0
+//      and the angle formula with the Luo (2004) version)
+```
+
+---
+
+## Step 2 — Implement the energy, gradient, and angle formula
+
+The Euclidean angle formula uses the law of cosines on edge lengths
+derived from log-lengths.  For inversive distance, the edge lengths
+are derived from inversive distances I_{ij} and the u-vector:
+
+```cpp
+// Euclidean (reference):
+//   lambda_ij = lambda0_ij + u_i + u_j
+//   l_ij      = exp(lambda_ij / 2)
+
+// Inversive distance (Luo 2004):
+//   cosh(l̃_ij) = I_ij * cosh((u_i + u_j) / 2)   [simplified form]
+//   l̃_ij       = acosh(I_ij * cosh((u_i + u_j) / 2))
+```
+
+Then the **corner angle** at vertex k opposite edge {i,j} in triangle {i,j,k}
+is computed by the standard law of cosines from l̃_{ij}, l̃_{jk}, l̃_{ki}.
+
+The **gradient** is the same as Euclidean:
+
+```cpp
+G_v = Σ_{faces containing v} alpha_v(face, u)  −  Theta_v
+```
+
+---
+
+## Step 3 — Write the gradient-check test
+
+Create `code/tests/cgal/test_inversive_distance.cpp`:
+
+```cpp
+#include "inversive_distance_functional.hpp"
+#include "mesh_factory.hpp"
+#include <gtest/gtest.h>
+
+// Finite-difference gradient check: copy the pattern from
+// test_euclidean_functional.cpp :: EuclideanFunctional.GradientCheck_Triangle
+TEST(InversiveDistance, GradientCheck_Triangle)
+{
+    // Build a single equilateral triangle (open mesh, no boundary issues).
+    ConformalMesh mesh = MeshFactory::make_open_mesh(MeshFactory::Kind::triangle);
+    InversiveDistanceMaps maps = setup_inversive_distance_maps(mesh);
+    compute_inversive_distance0_from_mesh(mesh, maps);  // I_ij from 3D edge lengths
+
+    const int n = /* count free DOFs */;
+    std::vector<double> x0(n, 0.0);
+    constexpr double eps = 1e-5;
+
+    auto G  = inversive_distance_gradient(mesh, x0, maps);
+
+    for (int i = 0; i < n; ++i) {
+        auto xp = x0; xp[i] += eps;
+        auto xm = x0; xm[i] -= eps;
+        double Ep = inversive_distance_energy(mesh, xp, maps);
+        double Em = inversive_distance_energy(mesh, xm, maps);
+        double fd = (Ep - Em) / (2.0 * eps);
+        EXPECT_NEAR(G[i], fd, 1e-6) << "gradient mismatch at DOF " << i;
+    }
+}
+```
+
+Run with:
+
+```bash
+./build/conformallab_cgal_tests --gtest_filter="InversiveDistance.*"
+```
+
+---
+
+## Step 4 — Register in CMakeLists
+
+Add to `code/tests/cgal/CMakeLists.txt` inside the `add_executable` block:
+
+```cmake
+    # ── Phase 9a: Inversive-distance functional ────────────────────────────────
+    test_inversive_distance.cpp
+```
+
+---
+
+## Step 5 — Add a Newton wrapper
+
+In `code/include/newton_solver.hpp` add:
+
+```cpp
+/// Solve the inversive-distance conformal problem.
+/// \see newton_euclidean() — identical structure.
+inline NewtonResult newton_inversive_distance(
+    ConformalMesh&                  mesh,
+    std::vector<double>             x0,
+    const InversiveDistanceMaps&    m,
+    double                          tol      = 1e-8,
+    int                             max_iter = 200)
+{
+    // Copy newton_euclidean() exactly, replacing euclidean_* with
+    // inversive_distance_*.  The Hessian is PSD (Luo 2004, Thm. 1),
+    // so SimplicialLDLT with SparseQR fallback applies directly.
+}
+```
+
+---
+
+## Step 6 — Verify against Java output
+
+The Java library outputs text results for a triangulated torus.  To compare:
+
+1. Run Java ConformalLab on the same OFF mesh with Inversive-Distance mode.
+2. Record the final gradient norm and angle sums.
+3. Run the C++ equivalent and check:
+
+```cpp
+EXPECT_LT(res.grad_inf_norm, 1e-8);
+EXPECT_LT(res.iterations, 50);  // inversive-distance converges faster than hyper-ideal
+```
+
+**Java reference:** `InversiveDistanceFunctionalTest.java` in `de.varylab.discreteconformal.test`
+
+---
+
+## Checklist
+
+- [ ] `code/include/inversive_distance_functional.hpp` compiles
+- [ ] `GradientCheck_Triangle` passes
+- [ ] `GradientCheck_OpenMesh` passes (copy from Euclidean tests)
+- [ ] Newton converges on `torus_4x4.off`
+- [ ] Result matches Java output (gradient norm < 1e-8 on same mesh)
+- [ ] Registered in `CMakeLists.txt`
+- [ ] `doc/roadmap/java-parity.md` updated: Phase 9a → ✅