docs: full audit — fix 4 wrong port/research labels + consolidated research-track

A full audit of `doc/` plus root-level markdown files (27 files) against
the actual ground truth in the C++ code and the local Java repository at
`/Users/tarikmoussa/Desktop/conformallab/` revealed four pre-existing
mis-labels and a stale test count.  All are corrected here.

Audit findings — corrected
─────────────────────────

1. **`InversiveDistanceFunctional` mis-labelled as Java port** (4 doc sites)
   Empirical verification:
       find /Users/tarikmoussa/Desktop/conformallab -iname "*nversive*"
       (zero matches)
   The class does NOT exist in `de.varylab.discreteconformal`.  The C++
   implementation is built from Luo 2004 + Glickenstein 2011 + Bowers-
   Stephenson 2004 — new research, not a port.
   Fixed in: java-parity.md, references.md, add-inversive-distance.md.

2. **HyperIdeal Hessian mis-labelled as "Java has analytic Hessian"**
   Empirical verification: `HyperIdealFunctional.java:295-298`:
       public boolean hasHessian() { return false; }
   Java has NO Hessian at all.  Both the FD (Phase 4a) and the block-FD
   (Phase 9b) Hessians in C++ are research beyond the Java port.  The
   chain rule (b,a) → ℓ → ζ → α/β is the *mathematical formulation*
   from Springborn 2020, not something Java implements.
   Fixed in: java-parity.md.

3. **Stale test count** README:87 said "28 suites, 170 tests" — current
   actual is 35 suites, 176 CGAL + 36 non-CGAL.  Fixed.

4. **Tutorial framing** — `add-inversive-distance.md` was framed as
   "porting an InversiveDistanceFunctional.java" that does not exist.
   Rewritten as "Implementing the Inversive-Distance functional from
   Luo 2004" with prominent verification block at top.

New document: `doc/roadmap/research-track.md`
─────────────────────────────────────────────

Consolidates everything in conformallab++ that goes beyond a Java port:

* Items already on `main`: HyperIdeal FD Hessian, period matrix τ
  partial-research components, Möbius holonomy storage.
* Items on open PRs: CP-Euclidean (PR #8, port), Inversive-Distance
  (PR #8, research), block-FD Hessian (PR #9, research).
* Planned research with full citations:
  - **Phase 9b-analytic** — full analytic HyperIdeal Hessian via
    Schläfli identity (Schläfli 1858/60) and chain rule through
    ζ₁₃/ζ₁₄/ζ₁₅, citing Springborn 2020 §4, Cho-Kim 1999,
    Glickenstein 2011 §4.  Includes acceptance-criteria checklist
    (per-case derivative cross-checks, gauge null space, PSD,
    measured ≥ 3× speed-up, LaTeX correctness note).
  - **Phase 9a.2-analytic** — analytic inversive-distance Hessian
    via Glickenstein 2011 eq. (4.6).
  - **Phase 10c** — full uniformization for genus g ≥ 2 (Fuchsian
    group representation) — fully new research, no Java reference.
  - **geometry-central** GC-1/2/3 exploratory track.

* Java backlog summary: 11 worth-porting Java classes identified by
  the parallel survey (FundamentalPolygonUtility, DiscreteHarmonicForm-
  Utility, DiscreteHolomorphicFormUtility, CanonicalBasisUtility,
  HyperbolicCyclicFunctional, QuasiisothermicUtility, KoebePolyhedron, …).
  ~6 500 Java lines, ~5 months of porting work, organised by phase.

Updated documents
─────────────────

* CLAUDE.md
  - New "Port-vs-research maintenance rule" with empirical verification
    command and the four corrected mis-labels.
  - Doc map: 23 → 24 documents (research-track.md added).

* README.md
  - Test count corrected (170 → 176+36).

* doc/math/references.md
  - Luo 2004 entry corrected ("new research" instead of "not yet ported").
  - New entries for Bowers-Stephenson 2004, Glickenstein 2011,
    Bobenko-Pinkall-Springborn 2010, Schläfli 1858/60.

* doc/roadmap/phases.md
  - Phase 9 reorganised: 9a split into 9a.1 (port) / 9a.2 (research),
    9b clarified as research (Java has no Hessian), 9c expanded with
    Java line counts and effort estimates.
  - Phase 10 reorganised: 10a/10b/10c with their Java prerequisites
    explicitly listed; 10c flagged as "fully new research".
  - Phase 10b' added: parallel research track (hyperbolic functional,
    quasi-isothermic, Möbius centering).
  - Phase 10c' added: optional Java-port additions (Koebe, circle
    patterns, electrostatic sphere).

* doc/roadmap/java-parity.md
  - Inversive-distance row: ❌ Java, ❌ C++ (Phase 9a.2) — new research.
  - CP-Euclidean row added: ✅ Java, ❌ C++ (Phase 9a.1) — port.
  - HyperIdeal Hessian row: ❌ Java, ⚠️ FD + block-FD in C++.
  - Worth-porting table replaced with the survey results (12 classes,
    Java line counts, suggested phases).
  - "HyperIdeal Hessian: FD vs analytic" section rewritten with the
    correction notice.

* doc/tutorials/add-inversive-distance.md
  - Rewritten end-to-end with prominent verification block at top.
  - Now correctly framed as "Implementing the Inversive-Distance
    functional from Luo 2004" — research, not port.
  - Includes the four required cross-validations:
    limit cases, Bowers-Stephenson round-trip, FD-vs-analytic,
    cross-validation against euclidean_functional at u=0.
  - New "How to know if it's a port or research" closing section
    with the empirical verification command.

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

doc/tutorials/add-inversive-distance.md

@@ -1,33 +1,84 @@
-# Tutorial: Porting the Inversive-Distance Functional (Phase 9a)
+# Tutorial: Implementing the Inversive-Distance Functional (Phase 9a.2)
 
-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).
+This tutorial walks through adding a **new** discrete-conformal functional
+to conformallab++.  The running example is the **vertex-based inversive-
+distance functional** of Luo (2004), used as Phase 9a.2 of the roadmap.
+
+> ## ⚠️ This is research, not a port
+>
+> An earlier draft of this document claimed this functional was a port of
+> `de.varylab.discreteconformal.functional.InversiveDistanceFunctional`.
+> **That Java class does not exist.**  Verified empirically:
+>
+> ```bash
+> $ find /Users/tarikmoussa/Desktop/conformallab -iname "*nversive*"
+> (zero results)
+> $ grep -r "InversiveDistance" /Users/tarikmoussa/Desktop/conformallab/src
+> (zero matches)
+> ```
+>
+> The closest Java cousin is `CPEuclideanFunctional.java`, which implements
+> the **face-based** circle-packing variant (Phase 9a.1).  The vertex-based
+> inversive-distance functional (this tutorial) is built **from the
+> literature**, not from a Java reference, and the correctness validation
+> is cross-checked against three sources:
+>
+> 1. **Luo, F.** (2004). *Combinatorial Yamabe Flow on Surfaces.* Comm. Contemp. Math. 6(5), 765–780.
+> 2. **Bowers, P. L. & Stephenson, K.** (2004). *Uniformizing dessins and Belyĭ maps via circle packing.* Mem. AMS 170(805).
+> 3. **Glickenstein, D.** (2011). *Discrete conformal variations and scalar curvature on piecewise flat manifolds.* J. Differential Geometry 87(2), 201–238.
+>
+> The tutorial below has been re-written to match this reality.
 
 **Prerequisite:** Read [doc/api/extending.md](../api/extending.md) first
-for the general pattern.  This tutorial fills in every detail for one
-specific case.
+for the general functional-porting pattern.  This tutorial fills in the
+mathematical and code details for one specific case.
 
 ---
 
 ## Mathematical background
 
-Inversive distance (Luo 2004) uses a different edge-length update formula:
+Inversive-distance circle packing parametrises **each vertex** by a circle
+of radius `r_i = exp(u_i)`.  Two adjacent circles have an *inversive
+distance* `I_ij` that is a fixed constant of the edge, derived once from
+the initial geometry via the Bowers–Stephenson identity:
 
 ```
-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) · ...
+I_ij = ( ℓ_ij² − r_i² − r_j² ) / ( 2 r_i r_j )           (Bowers-Stephenson 2004)
 ```
 
-For the Euclidean version the angle formula is the same law of cosines,
-but the log-lengths Λ̃ are computed differently from the u-vector.
+Geometric interpretation of `I_ij`:
 
-**Reference:** Luo, F. (2004). *Combinatorial Yamabe Flow on Surfaces*.
-Communications in Contemporary Mathematics, 6(5), 765–780.
+| Range | Configuration |
+|---|---|
+| `I_ij = +1`            | tangent circles (Koebe-style) |
+| `I_ij ∈ (0, 1)`        | overlapping with intersection angle `φ`, `I = cos φ` |
+| `I_ij = 0`             | orthogonal circles |
+| `I_ij ∈ (−1, 0)`       | disjoint circles, inversive distance > 1 |
+| `I_ij ≤ −1`            | impossible packing |
 
-**Java source:** `de.varylab.discreteconformal.functional.InversiveDistanceFunctional`
+The edge length under a state `u` is then determined by Luo's formula:
+
+```
+ℓ_ij(u)²  =  exp(2 u_i) + exp(2 u_j) + 2 I_ij exp(u_i + u_j)
+          =  r_i² + r_j² + 2 I_ij r_i r_j                       (Luo 2004 §3)
+```
+
+The angle formula is the same numerically-stable half-tangent law of
+cosines used by `euclidean_functional.hpp`; only the way `ℓ_ij` is
+computed from `u` is different.
+
+The gradient is the standard Yamabe-flow gradient:
+
+```
+∂E/∂u_v  =  Θ_v  −  Σ_{T ∋ v}  α_v(T)                          (Luo 2004 Lemma 3.1)
+```
+
+The energy is a path integral of the gradient (Luo's 1-form is closed on
+the domain where every triangle is valid); we use the same 10-point
+Gauss-Legendre quadrature as `euclidean_functional.hpp`.
+
+The Hessian is finite-difference for the MVP; an analytic form
+(Glickenstein 2011 eq. 4.6) is documented in the research-track roadmap.
 
 ---
 
@@ -38,165 +89,212 @@ cp code/include/euclidean_functional.hpp \
    code/include/inversive_distance_functional.hpp
 ```
 
-Edit the new file:
+Modify the maps struct: replace `lambda0` (Euclidean log-length) with the
+inversive-distance constant `I_e` and the initial radius `r0` (used for
+the Bowers-Stephenson init).
 
 ```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)
+    ConformalMesh::Property_map<Vertex_index, int>     v_idx;    // DOF index (−1 = pinned)
+    ConformalMesh::Property_map<Vertex_index, double>  theta_v;  // target cone angle
+    ConformalMesh::Property_map<Vertex_index, double>  r0;       // initial radius r_i^(0)
+    ConformalMesh::Property_map<Edge_index,   double>  I_e;      // inversive distance per edge
 };
-
-// ... (follow euclidean_functional.hpp exactly, replacing lambda0 with inv_dist0
-//      and the angle formula with the Luo (2004) version)
 ```
 
+Then implement the four entry points that any functional needs in
+conformallab++:
+
+- `setup_inversive_distance_maps(mesh)` — create maps with defaults.
+- `compute_inversive_distance_init_from_mesh(mesh, m)` — choose `r_i^(0)`
+  from the input geometry, then compute `I_ij` via Bowers-Stephenson.
+- `inversive_distance_gradient(mesh, x, m)` — Luo's `Θ − Σ α`.
+- `inversive_distance_energy(mesh, x, m)` — 10-point Gauss-Legendre path integral.
+
+For the full implementation, see `code/include/inversive_distance_functional.hpp`
+(part of PR #8).
+
 ---
 
-## Step 2 — Implement the energy, gradient, and angle formula
+## Step 2 — Edge-length kernel
 
-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:
+The single new pure-math primitive is the Luo edge-length formula.  Wrap
+it in a small detail helper so the gradient function reads cleanly:
 
 ```cpp
-// Euclidean (reference):
-//   lambda_ij = lambda0_ij + u_i + u_j
-//   l_ij      = exp(lambda_ij / 2)
+namespace id_detail {
 
-// 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))
+// ℓ² = exp(2u_i) + exp(2u_j) + 2 I exp(u_i + u_j)
+// Returns -1 on degenerate input (no valid packing).
+inline double edge_length_squared(double u_i, double u_j, double I_ij) {
+    double ri = std::exp(u_i);
+    double rj = std::exp(u_j);
+    double l2 = ri*ri + rj*rj + 2.0 * I_ij * ri * rj;
+    return l2 > 0.0 ? l2 : -1.0;
+}
+
+} // namespace id_detail
 ```
 
-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
-```
+This is the only place where the inversive-distance model differs from
+the Euclidean one.  All downstream code (angle computation, gradient
+accumulation, energy integration) is structurally identical.
 
 ---
 
-## Step 3 — Write the gradient-check test
+## Step 3 — Reuse `euclidean_angles()`
 
-Create `code/tests/cgal/test_inversive_distance.cpp`:
+The half-tangent law of cosines is independent of how lengths were
+obtained.  Feed `log(ℓ²)` to the existing helper to compute the three
+corner angles per face:
 
 ```cpp
-#include "inversive_distance_functional.hpp"
-#include "mesh_factory.hpp"
-#include <gtest/gtest.h>
+auto fa = euclidean_angles(std::log(l12sq), std::log(l23sq), std::log(l31sq));
+```
 
-// 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
+This is the **non-trivial reuse** that justifies the structural similarity
+to the Euclidean functional — we get the law-of-cosines numerics for free,
+and only the edge-length input changes.
 
-    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);
+## Step 4 — Validation tests
 
-    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;
+The acceptance criteria for this functional are stricter than for a Java
+port because there is no reference implementation to compare against.  We
+need **three independent validations**:
+
+### 4.1  Limit-case edge lengths
+
+Each of Luo's special cases (`I = 1` tangent, `I = 0` orthogonal,
+`I = −1` inside-tangent) gives a closed-form `ℓ` that must be reproduced
+to machine precision:
+
+```cpp
+TEST(InversiveDistanceFunctional, EdgeLengthFormula_TangentialLimit) {
+    // r_i=1, r_j=2, I=1:  ℓ² = 1 + 4 + 2·1·1·2 = 9 ⇒ ℓ = 3 = r_i + r_j
+    double l2 = id_detail::edge_length_squared(0.0, std::log(2.0), 1.0);
+    EXPECT_NEAR(std::sqrt(l2), 3.0, 1e-12);
+}
+```
+
+Three such tests cover the diagnostic special cases.
+
+### 4.2  Bowers-Stephenson round-trip
+
+The initialisation `compute_inversive_distance_init_from_mesh` must be
+self-consistent: starting from `(ℓ_3d, r_i, r_j)` and computing `I_ij`,
+the round-trip back through Luo's formula must give the original `ℓ`.
+
+```cpp
+TEST(InversiveDistanceFunctional, BowersStephensonRoundTrip) {
+    auto mesh = make_triangle();
+    auto m = setup_inversive_distance_maps(mesh);
+    compute_inversive_distance_init_from_mesh(mesh, m);
+
+    for (auto e : mesh.edges()) {
+        double l_3d  = /* mesh 3-D edge length */;
+        double ri    = m.r0[mesh.source(mesh.halfedge(e))];
+        double rj    = m.r0[mesh.target(mesh.halfedge(e))];
+        double l_rec = std::sqrt(ri*ri + rj*rj + 2.0 * m.I_e[e] * ri * rj);
+        EXPECT_NEAR(l_rec, l_3d, 1e-12);
     }
 }
 ```
 
-Run with:
+### 4.3  FD-vs-analytic gradient check
 
-```bash
-./build/conformallab_cgal_tests --gtest_filter="InversiveDistance.*"
-```
+Standard pattern from every functional in conformallab++ — see
+`test_euclidean_functional.cpp`.  Compare the analytic gradient
+to a symmetric finite difference of the energy.
 
----
+### 4.4  Cross-validation against `euclidean_functional.hpp`
 
-## 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:
+At `u = 0`, both functionals reconstruct the input 3-D edge length
+exactly (Euclidean via `compute_lambda0`, inversive distance via
+Bowers-Stephenson).  Therefore the actual angle sums are identical,
+and the two gradients (with default `Θ_v = 2π`) must match component-wise:
+
+```cpp
+TEST(InversiveDistanceFunctional, AngleDefectAtU0_AgreesWithEuclideanAtU0) {
+    auto G_id = inversive_distance_gradient(mesh, /*x=0*/, m_id);
+    auto G_eu = euclidean_gradient            (mesh, /*x=0*/, m_eu);
+    for (size_t i = 0; i < G_id.size(); ++i)
+        EXPECT_NEAR(G_id[i], G_eu[i], 1e-10);
+}
+```
+
+This is the empirical statement of Glickenstein 2011 §5: different
+parametrisations of the same initial discrete metric produce the same
+Newton-time-zero gradient.
+
+---
+
+## Step 5 — Register the tests
+
+In `code/tests/cgal/CMakeLists.txt`:
+
+```cmake
+    # ── Phase 9a.2: InversiveDistance (Luo 2004 + Glickenstein 2011) ─────────
+    #    Vertex-based inversive-distance circle-packing functional.  No Java
+    #    reference; implemented from the literature.  Cross-validated against
+    #    EuclideanCyclicFunctional at the natural initial geometry (u = 0).
+    test_inversive_distance_functional.cpp
+```
+
+Run:
+
+```bash
+ctest --test-dir build -R "InversiveDistance" --output-on-failure
+```
+
+---
+
+## Step 6 — Newton solver
+
+Once the functional passes its tests, wire a Newton wrapper into
+`newton_solver.hpp`:
 
 ```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.
-}
+    int                             max_iter = 200);
 ```
 
+The body is structurally identical to `newton_euclidean()` — same
+SimplicialLDLT + SparseQR fallback, same termination test.  Only the
+inner gradient / Hessian calls differ.
+
 ---
 
-## Step 6 — Verify against Java output
+## Checklist for a new functional
 
-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`
+- [ ] `code/include/<name>_functional.hpp` compiles
+- [ ] Limit-case edge-length tests pass at machine precision
+- [ ] Round-trip identity (init ⇄ length formula) verified
+- [ ] FD-vs-analytic gradient check passes on triangle, quad strip, tetra
+- [ ] Cross-validation test against an existing functional at `u = 0`
+- [ ] Newton wrapper added to `newton_solver.hpp`
+- [ ] Registered in `code/tests/cgal/CMakeLists.txt`
+- [ ] `doc/roadmap/java-parity.md` updated (port status or research note)
+- [ ] `doc/math/references.md` extended with the primary paper(s)
+- [ ] If this is *new research* beyond Java: add an entry in
+      `doc/roadmap/research-track.md` with citations and acceptance criteria
 
 ---
 
-## Checklist
+## How to know if it's a port or research
 
-- [ ] `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 → ✅
+Run the local Java-repo check **first** before writing any tutorial doc:
+
+```bash
+find /Users/tarikmoussa/Desktop/conformallab -iname "*<feature>*"
+grep -r "<ClassName>" /Users/tarikmoussa/Desktop/conformallab/src
+```
+
+If both return zero matches, the feature is **not** in Java and any C++
+implementation is **new research**, not a port.  The tutorial framing and
+the `doc/roadmap/research-track.md` entry should reflect this from day one.