let domains author their own pass-2 estimator as formulas, not Python

Pass 2 previously only knew one estimator: a hardcoded physics model that
matches dimensions literally named platform/actuator/energy_storage. Any
domain outside that shape (e.g. archery) got all-zero estimates and failed
every combo. Domains can now declare free variables and per-metric formulas
as data instead; a safe AST-based evaluator (engine/formula.py, no eval())
resolves declared entity properties via dep(key, constraint_type) and
generalizes the existing hand-nested mass-budget search into an N-variable
recursive optimizer. Fully additive -- the legacy platform/actuator/
energy_storage path is untouched and still runs unchanged for every domain
that declares no formulas.

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

CLAUDE.md

@@ -60,7 +60,7 @@ tests/                     # pytest, uses seeded_repo fixture from conftest.py
 ## Data flow (pipeline passes)
 
 1. **Pass 1 — Constraints**: `ConstraintResolver.resolve()` → blocked/conditional/valid. Blocked combos get a result row and `continue`.
-2. **Pass 2 — Estimation**: LLM or `_stub_estimate()` → raw metric values. Saved immediately via `save_raw_estimates()` (normalized_score=NULL).
+2. **Pass 2 — Estimation**: `_estimate_physics()` (deterministic physics engine; estimator-only, no LLM) → raw metric values. Saved immediately via `save_raw_estimates()` (normalized_score=NULL).
 3. **Pass 3 — Scoring**: `Scorer.score_combination()` → log-normalized scores + weighted geometric mean composite. Saves via `save_scores()` + `save_result()`.
 4. **Pass 4 — LLM Review**: Only for above-threshold combos with an LLM provider. No real provider yet (only `MockLLMProvider`).
 5. **Pass 5 — Human Review**: Manual via web UI results page.

LOGIC DOCS/003-gpu-batching-for-scaled-optimizer.md

@@ -0,0 +1,34 @@
+# GPU batching for the mass-allocation optimizer — not yet, here's the threshold
+
+## Context
+
+`Pipeline._decide_masses`'s joint platform/actuator/storage optimizer (coarse-to-fine grid search, see `_search_best_allocation`) calls its objective function roughly 11,700 times per combo. Profiling confirmed this dominates pipeline runtime: for 50 combos, 582,920 objective-function calls, each doing scalar arithmetic (power_density, the drag cubic solve, normalize, composite_score) on one `(platform, actuator, storage)` triple. The cost is Python's per-call overhead (bytecode dispatch, refcounting, attribute lookups), not the arithmetic itself — the individual formulas are cheap.
+
+## Why GPU doesn't fit today
+
+A single combo's grid is only ~169 points per round (13×13). GPUs pay off when there's enough independent parallel work to amortize kernel-launch and host↔device transfer overhead (each typically tens of microseconds to low milliseconds); 169 elements doesn't come close, and that overhead would be paid repeatedly — once per grid round, ~6-10 rounds per combo.
+
+The parallelism that actually exists is **across combos**, not within one combo's grid — every combo's optimization is fully independent of every other's. At current scale (~180 combos reach the optimizer per domain after Pass 1 filtering), batching every combo's grid into one array gives ~180×169 ≈ 30K elements per round — borderline, probably a wash against plain CPU numpy.
+
+## The actual threshold
+
+Combo count scales **multiplicatively** with added dimensions or entities per dimension (today: 11 platforms × 15 actuators × 18 storages ≈ 2,970 combos, ~180 of which reach the optimizer). Add a 4th dimension with even 10 options and total combos scale to ~30,000, with optimizer-eligible combos likely growing roughly proportionally to ~1,800/domain — batched grid size ≈ 300K elements/round. A 5th dimension does it again, into the low millions. That's the regime where a GPU's thousands of cores start meaningfully outrunning a CPU's 4-16-wide SIMD lanes.
+
+So: not "more dimensions" directly, but the combo×grid batch size those dimensions produce. Rough rule of thumb from this discussion:
+
+- **Tens of thousands of elements/round** (current scale, or a modest one-dimension addition): plain CPU numpy vectorization is enough, no GPU.
+- **Hundreds of thousands to low millions**: GPU batching across combos starts being worth evaluating.
+
+## What GPU batching would actually require
+
+Not just "swap numpy for cupy." It means restructuring `_process_pass2` from combo-first (one combo through the optimizer at a time) to batch-first (a chunk of N combos' grids evaluated together as one array with a "combo" axis, broadcasting each combo's own constants — `k_act`, `k_med`, `e_dens`, drag coefficients, mass bounds — across that axis). That's a real architectural change, not a drop-in acceleration:
+
+- **CLAUDE.md documents the pipeline as deliberately combo-first**: "each combo goes through all requested passes before the next combo starts... Progress is persisted per-combo (crash-safe, resumable)." Batching means checkpointing per-*batch*, not per-combo — a real (if manageable) tradeoff against that resumability guarantee, not something that comes free alongside the speedup.
+- Every branch in `_raw_physics_from_masses` / `_solve_achievable_speed_mps` (biological floors, ambient energy forms, degenerate fallbacks, the cubic's edge cases) needs to become `np.where(condition, a, b)` instead of `if/else` — careful, error-prone translation work, not mechanical.
+
+## Decision
+
+Don't build this now — no current need at ~180 combos/domain. If dimension count grows enough to matter, do it in two steps:
+
+1. **CPU numpy vectorization first** (batch one combo's grid into arrays, evaluate with vectorized ops instead of a Python double-loop). This is needed regardless of GPU or not, since it's the same rewrite either way, and profiling suggests it could plausibly give 10-50x on its own by replacing ~11,700 Python calls/combo with a couple dozen numpy batch calls.
+2. **Re-profile at the new scale.** Only reach for GPU batching-across-combos if CPU numpy is still the dominant cost after step 1, and only once the batched element count is actually in GPU-favorable territory (see thresholds above) — this is a "measure, then decide" call, not something to build ahead of need.

src/physcom/db/repository.py

@@ -9,7 +9,7 @@ from datetime import datetime, timezone
 from typing import Sequence
 
 from physcom.models.entity import Dependency, Entity
-from physcom.models.domain import Domain, DomainConstraint, MetricBound
+from physcom.models.domain import Domain, DomainConstraint, FreeVariable, MetricBound, MetricFormula
 from physcom.models.combination import Combination
 
 
@@ -286,6 +286,10 @@ class Repository:
                     "INSERT OR IGNORE INTO domain_constraints (domain_id, key, value) VALUES (?, ?, ?)",
                     (domain.id, dc.key, val),
                 )
+        for fv in domain.free_variables:
+            self.add_free_variable(domain.id, fv, commit=False)
+        for mf in domain.metric_formulas:
+            self.add_metric_formula(domain.id, mf, commit=False)
         self.conn.commit()
         return domain
 
@@ -299,6 +303,31 @@ class Repository:
             by_key.setdefault(r["key"], []).append(r["value"])
         return [DomainConstraint(key=k, allowed_values=v) for k, v in by_key.items()]
 
+    def _load_free_variables(self, domain_id: int) -> list[FreeVariable]:
+        rows = self.conn.execute(
+            """SELECT id, name, sort_order, floor_formula, ceiling_formula
+               FROM domain_free_variables WHERE domain_id = ? ORDER BY sort_order""",
+            (domain_id,),
+        ).fetchall()
+        return [
+            FreeVariable(
+                id=r["id"], name=r["name"], sort_order=r["sort_order"],
+                floor_formula=r["floor_formula"], ceiling_formula=r["ceiling_formula"],
+            )
+            for r in rows
+        ]
+
+    def _load_metric_formulas(self, domain_id: int) -> list[MetricFormula]:
+        rows = self.conn.execute(
+            """SELECT id, metric_name, formula
+               FROM domain_metric_formulas WHERE domain_id = ? ORDER BY metric_name""",
+            (domain_id,),
+        ).fetchall()
+        return [
+            MetricFormula(id=r["id"], metric_name=r["metric_name"], formula=r["formula"])
+            for r in rows
+        ]
+
     def _load_domain(self, where: str, param: str | int) -> Domain | None:
         row = self.conn.execute(f"SELECT * FROM domains WHERE {where} = ?", (param,)).fetchone()
         if not row:
@@ -325,6 +354,8 @@ class Repository:
                 for w in weights
             ],
             constraints=self._load_domain_constraints(row["id"]),
+            free_variables=self._load_free_variables(row["id"]),
+            metric_formulas=self._load_metric_formulas(row["id"]),
         )
 
     def get_domain(self, name: str) -> Domain | None:
@@ -376,12 +407,63 @@ class Repository:
         )
         self.conn.commit()
 
+    # ── Free variables & metric formulas ──────────────────────────
+
+    def add_free_variable(self, domain_id: int, fv: FreeVariable, commit: bool = True) -> FreeVariable:
+        cur = self.conn.execute(
+            """INSERT INTO domain_free_variables
+               (domain_id, name, sort_order, floor_formula, ceiling_formula)
+               VALUES (?, ?, ?, ?, ?)""",
+            (domain_id, fv.name, fv.sort_order, fv.floor_formula, fv.ceiling_formula),
+        )
+        fv.id = cur.lastrowid
+        if commit:
+            self.conn.commit()
+        return fv
+
+    def update_free_variable(self, fv_id: int, fv: FreeVariable) -> None:
+        self.conn.execute(
+            """UPDATE domain_free_variables
+               SET name = ?, sort_order = ?, floor_formula = ?, ceiling_formula = ?
+               WHERE id = ?""",
+            (fv.name, fv.sort_order, fv.floor_formula, fv.ceiling_formula, fv_id),
+        )
+        self.conn.commit()
+
+    def delete_free_variable(self, fv_id: int) -> None:
+        self.conn.execute("DELETE FROM domain_free_variables WHERE id = ?", (fv_id,))
+        self.conn.commit()
+
+    def add_metric_formula(self, domain_id: int, mf: MetricFormula, commit: bool = True) -> MetricFormula:
+        cur = self.conn.execute(
+            """INSERT OR REPLACE INTO domain_metric_formulas (domain_id, metric_name, formula)
+               VALUES (?, ?, ?)""",
+            (domain_id, mf.metric_name, mf.formula),
+        )
+        mf.id = cur.lastrowid
+        if commit:
+            self.conn.commit()
+        return mf
+
+    def update_metric_formula(self, mf_id: int, mf: MetricFormula) -> None:
+        self.conn.execute(
+            "UPDATE domain_metric_formulas SET metric_name = ?, formula = ? WHERE id = ?",
+            (mf.metric_name, mf.formula, mf_id),
+        )
+        self.conn.commit()
+
+    def delete_metric_formula(self, mf_id: int) -> None:
+        self.conn.execute("DELETE FROM domain_metric_formulas WHERE id = ?", (mf_id,))
+        self.conn.commit()
+
     def delete_domain(self, domain_id: int) -> None:
         self.conn.execute("DELETE FROM pipeline_runs WHERE domain_id = ?", (domain_id,))
         self.conn.execute("DELETE FROM combination_results WHERE domain_id = ?", (domain_id,))
         self.conn.execute("DELETE FROM combination_scores WHERE domain_id = ?", (domain_id,))
         self.conn.execute("DELETE FROM domain_metric_weights WHERE domain_id = ?", (domain_id,))
         self.conn.execute("DELETE FROM domain_constraints WHERE domain_id = ?", (domain_id,))
+        self.conn.execute("DELETE FROM domain_free_variables WHERE domain_id = ?", (domain_id,))
+        self.conn.execute("DELETE FROM domain_metric_formulas WHERE domain_id = ?", (domain_id,))
         self.conn.execute("DELETE FROM domains WHERE id = ?", (domain_id,))
         self.conn.commit()
 
@@ -475,7 +557,7 @@ class Repository:
         self, combo_id: int, status: str, block_reason: str | None = None, commit: bool = True
     ) -> None:
         # Don't downgrade from higher pass states — preserves human/LLM review data
-        if status in ("scored", "llm_reviewed") or status.endswith("_fail"):
+        if status in ("scored", "llm_reviewed", "valid") or status.endswith("_fail"):
             row = self.conn.execute(
                 "SELECT status FROM combinations WHERE id = ?", (combo_id,)
             ).fetchone()
@@ -488,6 +570,13 @@ class Repository:
                     return
                 if status == "llm_reviewed" and cur == "reviewed":
                     return
+                # "valid" is pass 1's domain-agnostic result -- a combo
+                # already at any later pass state (or a fail state) has
+                # progressed past pass 1 already, in this domain or
+                # another one sharing the same combo. Pass 1 re-running
+                # for a different domain must not silently revert that.
+                if status == "valid" and cur not in (None, "valid"):
+                    return
         self.conn.execute(
             "UPDATE combinations SET status = ?, block_reason = ? WHERE id = ?",
             (status, block_reason, combo_id),
@@ -899,6 +988,8 @@ class Repository:
         self.conn.execute("DELETE FROM entities")
         self.conn.execute("DELETE FROM domain_metric_weights")
         self.conn.execute("DELETE FROM domain_constraints")
+        self.conn.execute("DELETE FROM domain_free_variables")
+        self.conn.execute("DELETE FROM domain_metric_formulas")
         self.conn.execute("DELETE FROM domains")
         self.conn.execute("DELETE FROM metrics")
         self.conn.execute("DELETE FROM dimensions")

src/physcom/db/schema.py

@@ -120,6 +120,24 @@ CREATE TABLE IF NOT EXISTS domain_constraints (
     UNIQUE(domain_id, key, value)
 );
 
+CREATE TABLE IF NOT EXISTS domain_free_variables (
+    id              INTEGER PRIMARY KEY AUTOINCREMENT,
+    domain_id       INTEGER NOT NULL REFERENCES domains(id),
+    name            TEXT NOT NULL,
+    sort_order      INTEGER NOT NULL,
+    floor_formula   TEXT NOT NULL,
+    ceiling_formula TEXT NOT NULL,
+    UNIQUE(domain_id, name)
+);
+
+CREATE TABLE IF NOT EXISTS domain_metric_formulas (
+    id          INTEGER PRIMARY KEY AUTOINCREMENT,
+    domain_id   INTEGER NOT NULL REFERENCES domains(id),
+    metric_name TEXT NOT NULL,
+    formula     TEXT NOT NULL,
+    UNIQUE(domain_id, metric_name)
+);
+
 CREATE INDEX IF NOT EXISTS idx_deps_entity ON dependencies(entity_id);
 CREATE INDEX IF NOT EXISTS idx_deps_category_key ON dependencies(category, key);
 CREATE INDEX IF NOT EXISTS idx_combo_status ON combinations(status);

src/physcom/engine/constraint_resolver.py

@@ -32,6 +32,13 @@ CATEGORY_SEVERITY: dict[str, str] = {
     "energy": "block",
     "environment": "block",
     "infrastructure": "skip",
+    # Safety-critical physical necessities (radiation shielding, containment,
+    # etc.) -- same severity as energy/environment, not the softer default
+    # "warn" every other category falls through to. Missing this entry meant
+    # Nuclear Thermal Drive/Nuclear Fuel's "material" requires (radiation_
+    # shielding) defaulted to a non-blocking warning nothing in the catalog
+    # ever satisfies -- see LOGIC DOCS/002's "silent guardrail hole" pattern.
+    "material": "block",
 }
 
 # For provides-vs-range_min: deficit > this ratio = hard block, else warning
@@ -53,6 +60,40 @@ KEY_AGGREGATION: dict[str, str] = {
 OVERRUN_TOLERANCE: float = 0.10
 
 
+def aggregate_dependency_value(
+    combination: Combination,
+    key: str,
+    constraint_type: str,
+    key_aggregation: dict[str, str] | None = None,
+) -> float | None:
+    """Collapse every numeric dependency matching (key, constraint_type)
+    across a combination's entities into one system-level number: summed for
+    extensive keys (KEY_AGGREGATION says "sum", e.g. mass/footprint --
+    independent components sharing one physical vehicle), otherwise the
+    strongest single value wins (today's default pairwise behavior). Returns
+    None if no entity declares a matching numeric dependency. Shared by
+    ConstraintResolver._check_provides_vs_range and the formula evaluator's
+    injected dep() builtin (see engine/formula.py, engine/pipeline.py) so
+    there's one implementation of "how do these entities' declared numbers
+    combine," not two.
+    """
+    key_aggregation = KEY_AGGREGATION if key_aggregation is None else key_aggregation
+    values: list[float] = []
+    for entity in combination.entities:
+        for dep in entity.dependencies:
+            if dep.key != key or dep.constraint_type != constraint_type:
+                continue
+            try:
+                values.append(float(dep.value))
+            except (ValueError, TypeError):
+                continue
+    if not values:
+        return None
+    if key_aggregation.get(key) == "sum":
+        return sum(values)
+    return max(values)
+
+
 @dataclass
 class ConstraintResult:
     """Outcome of constraint resolution for a combination."""
@@ -243,11 +284,13 @@ class ConstraintResolver:
                     required.setdefault(dep.key, []).append((entity.name, val))
 
         for key in set(provided) & set(required):
+            prov_val = aggregate_dependency_value(
+                combination, key, "provides", self.key_aggregation
+            )
             if self.key_aggregation.get(key) == "sum":
                 prov_name = " + ".join(name for name, _ in provided[key])
-                prov_val = sum(val for _, val in provided[key])
             else:
-                prov_name, prov_val = max(provided[key], key=lambda t: t[1])
+                prov_name, _ = max(provided[key], key=lambda t: t[1])
 
             for req_name, req_val in required[key]:
                 if prov_val < req_val * self.deficit_threshold:

src/physcom/engine/formula.py

@@ -0,0 +1,139 @@
+"""Safe arithmetic expression language for domain-authored estimator formulas.
+
+No eval()/exec() anywhere -- compile_formula validates every AST node against
+a fixed whitelist (arithmetic, numeric/string constants, name lookups, calls
+to an explicitly supplied function table) before evaluate_formula ever walks
+it, so a formula can express "how much drawback_force does this combo
+produce" but never anything with attribute/subscript/import access.
+"""
+
+from __future__ import annotations
+
+import ast
+import math
+import operator
+from dataclasses import dataclass
+from typing import Callable
+
+
+class FormulaError(Exception):
+    """Raised for invalid formula syntax/structure or a failed evaluation."""
+
+
+_ALLOWED_NODES = (
+    ast.Expression, ast.BinOp, ast.UnaryOp, ast.Constant, ast.Name, ast.Load,
+    ast.Call, ast.keyword,
+    ast.Add, ast.Sub, ast.Mult, ast.Div, ast.Pow, ast.USub, ast.UAdd,
+)
+
+_BINOPS: dict[type, Callable[[float, float], float]] = {
+    ast.Add: operator.add,
+    ast.Sub: operator.sub,
+    ast.Mult: operator.mul,
+    ast.Div: operator.truediv,
+    ast.Pow: operator.pow,
+}
+
+_UNARYOPS: dict[type, Callable[[float], float]] = {
+    ast.USub: operator.neg,
+    ast.UAdd: operator.pos,
+}
+
+DEFAULT_FUNCTIONS: dict[str, Callable[..., float]] = {
+    "min": min,
+    "max": max,
+    "abs": abs,
+    "sqrt": math.sqrt,
+    "log": math.log,
+    "log1p": math.log1p,
+    "exp": math.exp,
+}
+
+
+@dataclass(frozen=True)
+class CompiledFormula:
+    source: str
+    _tree: ast.Expression
+
+
+def _validate(tree: ast.AST) -> None:
+    for node in ast.walk(tree):
+        if not isinstance(node, _ALLOWED_NODES):
+            raise FormulaError(
+                f"disallowed expression element: {type(node).__name__}"
+            )
+        if isinstance(node, ast.Constant):
+            if isinstance(node.value, bool) or not isinstance(node.value, (int, float, str)):
+                raise FormulaError(
+                    f"disallowed constant type: {type(node.value).__name__}"
+                )
+        if isinstance(node, ast.Name) and node.id.startswith("__"):
+            raise FormulaError(f"disallowed name: {node.id}")
+        if isinstance(node, ast.Call) and not isinstance(node.func, ast.Name):
+            raise FormulaError("only direct function calls are allowed")
+
+
+def compile_formula(source: str) -> CompiledFormula:
+    try:
+        tree = ast.parse(source, mode="eval")
+    except SyntaxError as exc:
+        raise FormulaError(f"invalid syntax in '{source}': {exc}") from exc
+    _validate(tree)
+    return CompiledFormula(source=source, _tree=tree)
+
+
+def _eval(node: ast.AST, variables: dict[str, float], functions: dict[str, Callable]):
+    if isinstance(node, ast.Expression):
+        return _eval(node.body, variables, functions)
+    if isinstance(node, ast.Constant):
+        # Numeric constants are cast to float (not left as int) so a formula
+        # like a**b**c can't build an arbitrary-precision giant int before
+        # ever raising -- float exponentiation overflows to inf/OverflowError
+        # quickly instead. String constants (dep() key/constraint_type args)
+        # pass through unchanged.
+        return node.value if isinstance(node.value, str) else float(node.value)
+    if isinstance(node, ast.Name):
+        if node.id not in variables:
+            raise FormulaError(f"unknown variable '{node.id}'")
+        return variables[node.id]
+    if isinstance(node, ast.BinOp):
+        op = _BINOPS.get(type(node.op))
+        if op is None:
+            raise FormulaError(f"unsupported operator: {type(node.op).__name__}")
+        return op(
+            _eval(node.left, variables, functions),
+            _eval(node.right, variables, functions),
+        )
+    if isinstance(node, ast.UnaryOp):
+        op = _UNARYOPS.get(type(node.op))
+        if op is None:
+            raise FormulaError(f"unsupported operator: {type(node.op).__name__}")
+        return op(_eval(node.operand, variables, functions))
+    if isinstance(node, ast.Call):
+        fname = node.func.id  # validated as ast.Name by _validate
+        func = functions.get(fname)
+        if func is None:
+            raise FormulaError(f"unknown function '{fname}'")
+        args = [_eval(a, variables, functions) for a in node.args]
+        kwargs = {kw.arg: _eval(kw.value, variables, functions) for kw in node.keywords}
+        return func(*args, **kwargs)
+    raise FormulaError(f"unsupported expression element: {type(node).__name__}")
+
+
+def evaluate_formula(
+    compiled: CompiledFormula,
+    variables: dict[str, float],
+    functions: dict[str, Callable] | None = None,
+) -> float:
+    effective_functions = {**DEFAULT_FUNCTIONS, **(functions or {})}
+    try:
+        result = _eval(compiled._tree, variables, effective_functions)
+    except FormulaError:
+        raise
+    except (TypeError, ValueError, ArithmeticError) as exc:
+        raise FormulaError(f"error evaluating '{compiled.source}': {exc}") from exc
+    if not isinstance(result, (int, float)) or isinstance(result, bool):
+        raise FormulaError(
+            f"formula '{compiled.source}' did not evaluate to a number"
+        )
+    return float(result)

src/physcom/engine/pipeline.py

@@ -106,13 +106,28 @@ ENERGY_FORM_RELIABILITY: dict[str, float] = {
 # validated it against the platform's ceiling), so it's used as the point
 # estimate rather than an invented one.
 
-# Human/animal actuators correctly declare mass_min=0 (a rider's body isn't
+# Radiation-pressure actuators (solar sails) don't declare a "mass" at
-# purchasable vehicle-borne mass and must not compete for the platform's
+# all -- thrust scales with sail area, not carried mass -- so their
-# mass budget), but that same 0 breaks power = power_density * mass. Fix:
+# effective mass is derived from declared footprint via a thin deployable
-# a fixed physiological reference mass used only in the power formula,
+# sail film's areal density. Used the same way as BIOLOGICAL_OPERATOR_MASS_KG
-# added to -- never substituted into -- the vehicle's own mass budget.
+# below: converts the entity's declared footprint FLOOR into a mass floor,
+# not a fixed value -- above it, effective mass is a free, budget-competing
+# variable like any other actuator (bigger sail = more collected power),
+# sized by the same joint optimizer, not a one-off product spec.
+SAIL_AREAL_DENSITY_KG_PER_M2: float = 0.05
+
+# Human/animal actuators declare mass_min=0 (there's no minimum purchase
+# quantity for a rider the way there is for an engine), but treated as a
+# literal floor that lets the optimizer size a payload down toward 0kg of
+# operator -- nonsensical, and it also breaks power = power_density * mass.
+# Used as a FLOOR (not a fixed value) on top of the declared mass_min: at
+# least one real operator must be present. Above that floor, actuator mass
+# is a free, budget-competing, structurally-carried variable exactly like
+# any mechanical actuator -- the "size" slider means more or bigger
+# operators (a loaded cargo trike, a two-horse team), sized by the same
+# joint optimizer everything else uses, not a fixed physiological constant.
 BIOLOGICAL_OPERATOR_MASS_KG: dict[str, float] = {
-    "biological": 70.0,  # human rider; Animal Traction shares this form too
+    "biological": 70.0,  # one average human rider; Animal Traction shares this form too
 }
 
 # A platform's declared mass range often spans a whole real-world class, not
@@ -182,13 +197,48 @@ def _solve_two_requirement_masses(
         return a_min, s_min
     return max(a, a_min), max(s, s_min)
 
-# Ambient energy forms (sun, wind, gravity, food) aren't a depletable
+
-# onboard store the way a fuel tank is -- "distance before running out"
+def _solve_achievable_speed_mps(
-# doesn't apply (a sailboat doesn't run out of wind). Rather than
+    power_density: float, floor_total: float, k_med: float, drag_coeff: float,
-# degenerate to 0 (mass_min=0, energy_density often undeclared entirely),
+) -> float:
-# range_fuel reports the domain's own declared ceiling for these: full
+    """Invert specific_power = k_med*v + (drag_coeff/floor_total)*v^3 for v
-# marks is the physically honest answer, not an error.
+    -- the steady-state speed at which a build's actual power output exactly
-AMBIENT_ENERGY_FORMS: set[str] = {"biological", "wind", "radiation_pressure", "gravitational"}
+    balances mass-proportional resistance plus mass-independent aerodynamic
+    drag. A depressed cubic (no v^2 term) with drag_coeff/floor_total > 0
+    and k_med >= 0: A*v^3 + B*v - C = 0 is strictly increasing for v >= 0
+    (derivative 3*A*v^2 + B > 0 everywhere), so it has exactly one
+    non-negative real root -- solved directly via Cardano's formula, no
+    iteration needed. Falls back to the plain linear model (v = power/k_med)
+    when there's no drag coefficient for this medium, so ungraded media
+    behave exactly as before."""
+    if power_density <= 0 or k_med is None:
+        return 0.0
+    if drag_coeff <= 0 or floor_total <= 0:
+        return power_density / k_med if k_med else 0.0
+    A = drag_coeff / floor_total
+    B = k_med
+    C = power_density
+    p, q = B / A, -C / A
+
+    def cbrt(x: float) -> float:
+        return math.copysign(abs(x) ** (1 / 3), x) if x else 0.0
+
+    discriminant = (q / 2) ** 2 + (p / 3) ** 3  # always >= 0 given p, C >= 0
+    sqrt_disc = math.sqrt(discriminant)
+    v = cbrt(-q / 2 + sqrt_disc) + cbrt(-q / 2 - sqrt_disc)
+    return max(v, 0.0)
+
+
+# Ambient energy forms (sun, wind, gravity) aren't a depletable onboard
+# store the way a fuel tank is -- "distance before running out" doesn't
+# apply (a sailboat doesn't run out of wind). Rather than degenerate to 0
+# (mass_min=0, energy_density often undeclared entirely), range_fuel
+# reports the domain's own declared ceiling for these: full marks is the
+# physically honest answer, not an error. Food is deliberately NOT here:
+# stopping to eat is a resupply, the same category as refuelling a tank,
+# not a genuinely external/inexhaustible power source -- Biological Feed
+# uses the normal storage-mass-limited range_fuel formula.
+AMBIENT_ENERGY_FORMS: set[str] = {"wind", "radiation_pressure", "gravitational"}
 
 # Resistive energy cost of travel, J per kg of vehicle per meter --
 # rolling resistance for ground vehicles, cruise-flight lift/drag for
@@ -210,24 +260,52 @@ SPECIFIC_ENERGY_CONSUMPTION_J_PER_KG_M: dict[str, float] = {
 # so it falls through to the old placeholder formula in the code below
 # rather than silently claiming a resistance-based number that isn't real.
 #
-# KNOWN GAP: this whole table is mass-proportional resistance only (rolling
+# FORMERLY A KNOWN GAP, now fixed below: the table above is mass-proportional
-# resistance, effectively) -- there's no aerodynamic drag term (force ~
+# resistance only (rolling resistance, effectively) -- no aerodynamic drag
-# frontal_area * velocity^2, independent of mass). That's a reasonable
+# term (force ~ frontal_area * velocity^2, independent of mass). That's a
-# approximation for something car-scale, where rolling resistance genuinely
+# reasonable approximation for something car-scale, where rolling resistance
-# dominates at typical speeds and this was validated against real car range.
+# genuinely dominates at typical speeds and this was validated against real
-# It badly overestimates range for light/human-scale vehicles, where drag
+# car range. It badly overestimated range for light/human-scale vehicles,
-# is the dominant resistance term and doesn't scale down with mass the way
+# where drag is the dominant resistance term and doesn't scale down with
-# this formula assumes -- confirmed on a real combo (Light Personal Vehicle +
+# mass the way this formula assumes -- confirmed on a real combo (Light
-# Electric Motor + Rechargeable Battery, #876): a sane 9kg battery on a
+# Personal Vehicle + Electric Motor + Rechargeable Battery, #876): a sane
-# realistic 31kg vehicle came out to ~1,977km, a 6-9x overestimate against
+# 9kg battery on a realistic 31kg vehicle came out to ~1,977km, a 6-9x
-# real e-bikes on comparable battery energy (~50-80km on ~500Wh). The mass
+# overestimate against real e-bikes on comparable battery energy (~50-80km
-# allocation itself was fine (correctly floor-clamped, nothing oversized) --
+# on ~500Wh). It also meant an achieved-speed metric derived from power
-# this is a missing term in the resistance formula, not an allocation bug,
+# alone (see DRAG_POWER_COEFF_BY_MEDIUM / _solve_achievable_speed_mps below)
-# so a mass-allocation optimizer wouldn't fix it either. Real fix needs a
+# had no ceiling at all -- without a v^2-scaling force to push back, more
-# genuine drag term (frontal-area-ish figure -- `footprint` exists but is a
+# power always bought proportionally more speed, forever.
-# ground-footprint number, not obviously the right proxy for cross-sectional
+#
-# area facing the wind -- and a drag coefficient assumption), scoped
+# DRAG_POWER_COEFF_BY_MEDIUM below adds that missing term: a mass-INDEPENDENT
-# separately from the resistance-constant tuning already done here.
+# drag power coefficient (0.5 * air_density * drag_coefficient * frontal_area,
+# W per (m/s)^3) added on top of the existing mass-proportional term.
+#
+# "air" was originally left out of this table on the reasoning that its
+# L/D-based cruise model was already a reasonable velocity-roughly-linear
+# approximation -- true for computing energy per meter during a normal
+# cruise, but WRONG for the exact same reason ground was wrong: L/D-based
+# drag force is also roughly velocity-independent within a design cruise
+# band, so it's still just "resistance = mass-proportional constant" with
+# no v^2 term, and inverting power/resistance for achieved speed still had
+# no ceiling. Confirmed live: a Rotorcraft + Gas Turbine combo showed a
+# "speed" of 3,127 m/s (Mach 9) with phi4's pass-4 review flagging it
+# directly ("unrealistic for urban commuting, likely indicating an
+# error"). Added below with an aircraft-like reference cross-section.
+# Water hull drag would need its own (different) treatment rather than
+# reusing either reference, so it's left as a known remaining gap.
+DRAG_POWER_COEFF_BY_MEDIUM: dict[str, float] = {
+    # 0.5 * rho_air(1.225 kg/m^3) * Cd(~0.3) * frontal_area(~2.2 m^2, small
+    # car reference) -- sanity check: at 30 m/s (108 km/h) this alone costs
+    # ~11kW, in the right ballpark for real highway cruise power.
+    "ground": 0.5 * 1.225 * 0.3 * 2.2,
+    # 0.5 * rho_air(1.225) * Cd(~0.2, streamlined fuselage) * frontal_area
+    # (~1.74 m^2, small aircraft/rotorcraft reference) -- sanity check: at
+    # 60 m/s (a fast urban rotorcraft cruise) this alone costs ~46kW, a
+    # plausible fraction of a light helicopter's real cruise power (most
+    # of the rest goes to induced/rotor drag, not modeled here -- this
+    # coefficient only covers fuselage parasite drag).
+    "air": 0.5 * 1.225 * 0.2 * 1.74,
+}
 
 # Structural manufacturing cost, $ per kg of platform mass -- certification
 # and materials overhead scale hugely by medium (aerospace-grade vs.
@@ -596,7 +674,7 @@ class Pipeline:
             result.pass2_estimated += 1
             return
 
-        raw_metrics, feasible = self._stub_estimate(combo, domain.metric_bounds)
+        raw_metrics, feasible = self._estimate_physics(combo, domain)
 
         if not feasible:
             # No platform mass within its own declared ceiling could
@@ -629,7 +707,7 @@ class Pipeline:
                 estimate_dicts.append({
                     "metric_id": mb.metric_id,
                     "raw_value": rval,
-                    "estimation_method": "stub",
+                    "estimation_method": "physics_calc",
                     "confidence": 1.0,
                 })
         if estimate_dicts:
@@ -767,7 +845,7 @@ class Pipeline:
             self._wait_for_rate_limit(run_id, exc.retry_after)
             try:
                 review_result = self.llm.review_plausibility(
-                    description, raw_dict, score_dict, domain.metric_bounds
+                    description, raw_dict, score_dict, domain
                 )
             except LLMRateLimitError:
                 return  # still limited; skip, retry next run
@@ -812,9 +890,9 @@ class Pipeline:
         self, combo: Combination, bounds_by_name: dict[str, MetricBound]
     ) -> "_PhysicsContext | None":
         """Derive the entity-level physics inputs that don't depend on a
-        mass allocation choice -- shared by _stub_estimate (which picks the
+        mass allocation choice -- shared by _estimate_physics (which picks
-        allocation via solve or a special case) and _optimize_allocation
+        the allocation via _decide_masses) and evaluate_allocation (the
-        (which searches over candidate allocations). Returns None if the
+        explore-panel's direct evaluation). Returns None if the
         combo doesn't have the platform/actuator/storage shape this whole
         formula assumes (shouldn't happen for real combos, but a domain
         without all three dimensions requested would hit this)."""
@@ -859,28 +937,87 @@ class Pipeline:
         ctx: "_PhysicsContext",
         actuator_mass: float,
         storage_mass: float,
-        power_mass: float,
-        denom_offset: float,
         bounds_by_name: dict[str, MetricBound],
         units_by_name: dict[str, str],
-        cargo_capacity_kg: float,
         platform_mass: float | None = None,
     ) -> dict[str, float]:
-        """power_density/range_fuel/cost_efficiency for an EXPLICIT mass
+        """power_density/range_fuel/cost_efficiency/cargo_capacity for an
-        allocation. `power_mass` is separate from `actuator_mass` for the
+        EXPLICIT mass allocation. `platform_mass` defaults to the
-        biological/radiation-pressure special cases (see _stub_estimate),
+        platform's representative mass (ctx.p_rep) -- pass an explicit
-        where the numerator mass isn't the same as the build-budget mass;
+        value to explore a specific weight class instead (see
-        for the normal (solved, optimized, or manually-explored) case
+        evaluate_allocation). Cargo capacity is derived from THIS build's
-        they're the same value. `platform_mass` defaults to the platform's
+        actual platform+actuator mass (not a separate declared-floor
-        representative mass (ctx.p_rep) -- pass an explicit value to
+        constant), with storage_mass subtracted out of that allowance --
-        explore a specific weight class instead (see evaluate_allocation)."""
+        see the deadweight/lightship comment at its computation below for
+        why fuel/battery competes with cargo instead of padding it. It
+        responds to the same optimizer/explore-slider choices every other
+        metric here does."""
         p_mass = ctx.p_rep if platform_mass is None else platform_mass
         out: dict[str, float] = {}
         floor_total = p_mass + actuator_mass + storage_mass
-        physics_denom = floor_total + denom_offset
+        power_density_value = (ctx.k_act * actuator_mass) / floor_total if floor_total else 0.0
 
         if "power_density" in bounds_by_name:
-            out["power_density"] = (ctx.k_act * power_mass) / physics_denom if physics_denom else 0.0
+            out["power_density"] = power_density_value
+
+        # Deadweight/lightship cargo capacity. Real deadweight tonnage is a
+        # FIXED allowance sized off the vessel's own empty (lightship) mass
+        # -- hull + machinery, NOT fuel or cargo -- and fuel and cargo then
+        # SHARE that one allowance: a ship that bunkers more fuel has that
+        # much less room left for cargo, and vice versa. platform+actuator
+        # is the lightship analog here (the vehicle's own hardware);
+        # storage_mass is the fuel/battery competing with cargo for the
+        # same pool, not part of the base the pool is sized from -- get
+        # that backwards (basing the pool on platform+actuator+storage,
+        # as an earlier version of this did) and more battery looks like it
+        # BUYS more cargo room instead of using it up. Two ratio
+        # conventions coexist because heavy freight/maritime vehicles
+        # genuinely carry a much larger multiple of their own mass in
+        # cargo than light personal/delivery vehicles do (see
+        # CARGO_KG_PER_STRUCTURAL_KG's module comment); which one a domain
+        # scores is just which metric_name it declares. Floored at 0: a
+        # storage mass bigger than the whole allowance leaves no cargo
+        # room, not negative room.
+        # Skip the arithmetic entirely for domains that don't score either
+        # cargo convention and don't need it as cost_efficiency's $/(kg·m)
+        # denominator either -- this runs on every one of the ~11,000 grid
+        # points the search below tries per combo, so a domain like
+        # interplanetary_travel (scores neither) shouldn't pay for it.
+        needs_cargo = (
+            "cargo_capacity" in bounds_by_name
+            or "cargo_capacity_kg" in bounds_by_name
+            or units_by_name.get("cost_efficiency") == "$/(kg·m)"
+        )
+        cargo_capacity_2_5x = 0.0
+        if needs_cargo:
+            lightship_mass = p_mass + actuator_mass
+            cargo_capacity_2_5x = max(0.0, lightship_mass * CARGO_KG_PER_STRUCTURAL_KG - storage_mass)
+            cargo_capacity_0_3x = max(0.0, lightship_mass * 0.3 - storage_mass)
+            if "cargo_capacity" in bounds_by_name:
+                out["cargo_capacity"] = cargo_capacity_2_5x
+            if "cargo_capacity_kg" in bounds_by_name:
+                out["cargo_capacity_kg"] = cargo_capacity_0_3x
+
+        # Achieved steady-state cruise speed, DERIVED from this specific
+        # build's actual power_density, the medium's mass-proportional
+        # resistance, and (ground/air, see DRAG_POWER_COEFF_BY_MEDIUM) a
+        # mass-independent aerodynamic drag term -- not a platform-declared
+        # constant. A build with more power than the platform's bare
+        # target_velocity requires achieves a genuinely higher speed here;
+        # an underbuilt one achieves less -- speed is an output of the
+        # build, not an input to it. Computed unconditionally (not just
+        # when "speed" is a scored metric) because range_fuel/cost_efficiency
+        # below both need it too: the energy actually spent per meter
+        # depends on how fast this build is actually going, drag included.
+        drag_coeff = DRAG_POWER_COEFF_BY_MEDIUM.get(ctx.medium, 0.0)
+        achieved_speed = _solve_achievable_speed_mps(power_density_value, floor_total, ctx.k_med, drag_coeff)
+        effective_k_med = (
+            (ctx.k_med + drag_coeff * achieved_speed ** 2 / floor_total)
+            if ctx.k_med is not None and floor_total > 0 else ctx.k_med
+        )
+
+        if "speed" in bounds_by_name:
+            out["speed"] = achieved_speed
 
         if "range_fuel" in bounds_by_name:
             if ctx.storage_energy_form in AMBIENT_ENERGY_FORMS or ctx.k_med is None:
@@ -888,7 +1025,7 @@ class Pipeline:
                 out["range_fuel"] = mb.norm_max if mb else 0.0
             elif floor_total > 0:
                 out["range_fuel"] = min(
-                    (ctx.e_dens * storage_mass) / (ctx.k_med * floor_total), 1e13
+                    (ctx.e_dens * storage_mass) / (effective_k_med * floor_total), 1e13
                 )
 
         if "cost_efficiency" in bounds_by_name:
@@ -907,15 +1044,30 @@ class Pipeline:
             )
             amortized_per_m = upfront_cost / lifetime_m
 
+            if ctx.k_med is not None:
                 fuel_price_per_mj = FUEL_PRICE_PER_MJ.get(ctx.storage_energy_form, 0.04)
-            energy_per_m_mj = (
+                energy_per_m_mj = (effective_k_med * floor_total) / 1e6
-                (ctx.k_med or SPECIFIC_ENERGY_CONSUMPTION_J_PER_KG_M["ground"]) * floor_total
-            ) / 1e6
                 operating_per_m = energy_per_m_mj * fuel_price_per_mj
+            else:
+                # No resistance model for this medium (space -- real range
+                # is governed by the rocket equation, not implemented
+                # here, see the comment above
+                # SPECIFIC_ENERGY_CONSUMPTION_J_PER_KG_M). Don't fabricate
+                # an operating cost from ground physics the way an earlier
+                # version of this did (`effective_k_med or ...["ground"]`)
+                # -- report upfront/amortized cost only, an honest partial
+                # answer, rather than a wrong number for an unmodeled term.
+                operating_per_m = 0.0
 
             cost_per_m = amortized_per_m + operating_per_m
             if units_by_name.get("cost_efficiency") == "$/(kg·m)":
-                out["cost_efficiency"] = cost_per_m / max(cargo_capacity_kg, 1.0)
+                # Divide by whichever cargo convention this domain actually
+                # scores, so cost-per-cargo-kg and the cargo_capacity number
+                # shown alongside it always agree; default to the heavy-
+                # vehicle ratio if a domain scores $/(kg·m) without scoring
+                # either cargo metric explicitly (matches prior behavior).
+                cargo_basis = out.get("cargo_capacity_kg", out.get("cargo_capacity", cargo_capacity_2_5x))
+                out["cost_efficiency"] = cost_per_m / max(cargo_basis, 1.0)
             else:
                 out["cost_efficiency"] = cost_per_m
 
@@ -926,35 +1078,42 @@ class Pipeline:
         ctx: "_PhysicsContext",
         bounds_by_name: dict[str, MetricBound],
         units_by_name: dict[str, str],
-        cargo_capacity_kg: float,
+    ) -> tuple[float, float, float, bool]:
-    ) -> tuple[float, float, float, float, float, bool]:
         """Pick the platform/actuator/storage mass for the build this domain
         actually scores. First, the platform's declared physical
         performance target (accel/thrust, or target_velocity/resistance)
         sets a FLOOR -- a rotorcraft that can't produce enough thrust to
         hover isn't a rotorcraft, regardless of how a smaller/cheaper
-        engine might score. That floor also sets the smallest platform
+        engine might score. Biological actuators (a rider's own body) and
-        mass that could structurally carry it (CARGO_KG_PER_STRUCTURAL_KG
+        radiation-pressure actuators (a solar sail) get the same treatment
-        again, applied to the platform carrying its own actuator+storage
+        with one addition: BIOLOGICAL_OPERATOR_MASS_KG / a footprint-derived
-        instead of cargo) -- below that, no actuator/storage choice is
+        floor (see SAIL_AREAL_DENSITY_KG_PER_M2) sets a floor under the
-        physically possible. Above that lower bound, platform mass is a
+        floor -- at least one real operator, or the sail's own declared
-        real THIRD search variable, not fixed at p_rep: a bigger platform
+        minimum footprint, even if the performance-derived requirement
-        also raises the structural cap on how much actuator+storage it can
+        would otherwise ask for less -- but above that, mass is a free
-        carry, so growing all three together can score higher than
+        variable exactly like a mechanical actuator's; "bigger" means more
-        minimizing platform down to what's merely required. Searched
+        or bigger operators, or a bigger sail, not a fixed constant. That
-        jointly (outer coarse-to-fine scan over platform mass, inner
+        floor also
-        coarse-to-fine scan over actuator/storage at each candidate) for
+        sets the smallest platform mass that could structurally carry it
-        whatever allocation maximizes this domain's own weighted composite
+        (CARGO_KG_PER_STRUCTURAL_KG again, applied to the platform
-        score, using the same normalize()/composite_score() the real
+        carrying its own actuator+storage instead of cargo) -- below that,
-        scoring pass uses. Not "just enough to function" and not "best
+        no actuator/storage choice is physically possible. Above that
-        score regardless of function" -- both, floor then optimize jointly.
+        lower bound, platform mass is a real THIRD search variable, not
-        Returns (actuator_mass, storage_mass, power_mass, denom_offset,
+        fixed at p_rep: a bigger platform also raises the structural cap
-        platform_mass, feasible); see _raw_physics_from_masses for what
+        on how much actuator+storage it can carry, so growing all three
-        power_mass and denom_offset mean. `feasible` is False only when no
+        together can score higher than minimizing platform down to what's
-        platform mass within its own declared ceiling could structurally
+        merely required. Searched jointly (outer coarse-to-fine scan over
-        carry the required floor -- power_density/range_fuel/cost_efficiency
+        platform mass, inner coarse-to-fine scan over actuator/storage at
-        are all per-kg ratios, so they don't naturally penalize a build
+        each candidate) for whatever allocation maximizes this domain's
-        whose absolute mass tramples its own platform's declared ceiling;
+        own weighted composite score, using the same normalize()/
+        composite_score() the real scoring pass uses. Not "just enough to
+        function" and not "best score regardless of function" -- both,
+        floor then optimize jointly. Returns (actuator_mass, storage_mass,
+        platform_mass, feasible). `feasible` is False only when no platform
+        mass within its own declared ceiling could structurally carry the
+        required floor -- power_density/range_fuel/cost_efficiency are all
+        per-kg ratios, so they don't naturally penalize a build whose
+        absolute mass tramples its own platform's declared ceiling;
         callers must treat an infeasible build as a hard fail rather than
         trusting the (still-computable, still ratio-plausible) score. Also
         used by evaluate_allocation to compute the slider's starting
@@ -967,17 +1126,6 @@ class Pipeline:
                     return float(dep.value)
             return None
 
-        if ctx.actuator_energy_form in BIOLOGICAL_OPERATOR_MASS_KG:
-            power_mass = BIOLOGICAL_OPERATOR_MASS_KG[ctx.actuator_energy_form]
-            return ctx.a_min, ctx.s_min, power_mass, power_mass, ctx.p_rep, True
-        if ctx.actuator_energy_form == "radiation_pressure":
-            # thrust scales with sail area, not carried mass -- derive an
-            # effective mass from declared footprint and a thin-film areal
-            # density estimate rather than the (undeclared) mass attribute.
-            footprint = dep_value(ctx.actuator, "footprint", "range_min") or 0.0
-            actuator_mass = footprint * 0.05  # kg/m^2, thin deployable sail film
-            return actuator_mass, ctx.s_min, actuator_mass, 0.0, ctx.p_rep, True
-
         # Step 1: the required floor (same solve as before -- now a floor
         # for the search below, not the final answer).
         min_accel = dep_value(ctx.platform, "min_effective_accel", "range_min")
@@ -986,7 +1134,12 @@ class Pipeline:
         range_bounds = bounds_by_name.get("range_fuel")
         target_range = range_bounds.norm_max if range_bounds else None
 
-        if min_accel and specific_thrust:
+        # `is not None`, not truthy -- dep_value() legitimately returns 0.0
+        # for a declared floor of zero (Spaceship declares
+        # min_effective_accel=0, a real "no acceleration floor" value, not
+        # "undeclared"). A truthy check would silently treat that the same
+        # as an absent requirement and fall through to the wrong branch.
+        if min_accel is not None and specific_thrust is not None:
             c1, r1 = specific_thrust, min_accel
         elif target_velocity and ctx.k_med:
             # Resistance alone (k_med) only covers steady-state cruise --
@@ -1019,10 +1172,23 @@ class Pipeline:
             required_actuator = ctx.a_min if ctx.a_min > 0.0 else 10.0
             required_storage = ctx.s_min
 
+        if ctx.actuator_energy_form in BIOLOGICAL_OPERATOR_MASS_KG:
+            # At least one real operator, regardless of what the bare
+            # performance solve above would have asked for -- see the
+            # BIOLOGICAL_OPERATOR_MASS_KG module comment.
+            required_actuator = max(required_actuator, BIOLOGICAL_OPERATOR_MASS_KG[ctx.actuator_energy_form])
+        elif ctx.actuator_energy_form == "radiation_pressure":
+            # At least the entity's own declared minimum sail footprint,
+            # regardless of what the bare performance solve above would
+            # have asked for -- see the SAIL_AREAL_DENSITY_KG_PER_M2
+            # module comment.
+            footprint_floor = dep_value(ctx.actuator, "footprint", "range_min") or 0.0
+            required_actuator = max(required_actuator, footprint_floor * SAIL_AREAL_DENSITY_KG_PER_M2)
+
         if ctx.p_max is None:
             # No declared mass ceiling (e.g. Spaceship) -- no bounded
             # budget to search within, use the requirement floor as-is.
-            return required_actuator, required_storage, required_actuator, 0.0, ctx.p_rep, True
+            return required_actuator, required_storage, ctx.p_rep, True
 
         a_floor = max(ctx.a_min, required_actuator)
         s_floor = max(ctx.s_min, required_storage)
@@ -1048,12 +1214,12 @@ class Pipeline:
             # per-kg ratios, so they don't naturally penalize a build whose
             # ABSOLUTE mass tramples its own platform's declared ceiling --
             # something else has to catch that).
-            return a_floor, s_floor, a_floor, 0.0, p_lo, False
+            return a_floor, s_floor, p_lo, False
 
         def objective(platform_mass: float, actuator_mass: float, storage_mass: float) -> float:
             raw = self._raw_physics_from_masses(
-                ctx, actuator_mass, storage_mass, actuator_mass, 0.0,
+                ctx, actuator_mass, storage_mass,
-                bounds_by_name, units_by_name, cargo_capacity_kg,
+                bounds_by_name, units_by_name,
                 platform_mass=platform_mass,
             )
             scores, weights = [], []
@@ -1115,7 +1281,7 @@ class Pipeline:
                 best_score, best_p = sc, p
 
         actuator_mass, storage_mass, _score = best_at_platform(best_p, grid=12, rounds=6)
-        return actuator_mass, storage_mass, actuator_mass, 0.0, best_p, True
+        return actuator_mass, storage_mass, best_p, True
 
     @staticmethod
     def _search_best_allocation(
@@ -1157,21 +1323,32 @@ class Pipeline:
 
         return best_a, best_s, best_score
 
-    def _stub_estimate(
+    def _estimate_physics(
-        self, combo: Combination, metric_bounds: list[MetricBound]
+        self, combo: Combination, domain: Domain
     ) -> tuple[dict[str, float], bool]:
-        """Deterministic estimation from declared entity attributes (no LLM).
+        """Deterministic physics-based estimation from declared entity
+        attributes (no LLM) -- pass 2's estimator.
 
-        power_density, range_fuel, and cost_efficiency are computed from the
+        Domains that declare their own metric_formulas (see
-        platform's declared mass envelope treated as a combo-wide budget —
+        _estimate_via_formulas) are estimated entirely by those
-        see the module-level comment above BIOLOGICAL_OPERATOR_MASS_KG for
+        domain-authored formulas instead -- this method's hardcoded
-        the full formula rationale.
+        platform/actuator/energy_storage physics model below is additive,
+        not the only path: it's untouched and still runs unchanged for
+        every domain that declares no formulas (every existing transport
+        domain today).
 
-        safety/availability/reliability/cargo_capacity/environmental_impact
+        power_density, speed, range_fuel, cost_efficiency, and
-        are untouched — these are judgment calls (regulatory, economic,
+        cargo_capacity/cargo_capacity_kg are all computed from the
-        qualitative), not physics, and stay on the categorical lookup-table
+        platform's declared mass envelope treated as a combo-wide budget,
-        heuristics below (actuator's thrust_profile and energy_form and the
+        jointly optimized by _decide_masses -- see the module-level
-        combo's infrastructure requirements).
+        comment above BIOLOGICAL_OPERATOR_MASS_KG for the full formula
+        rationale.
+
+        safety/availability/reliability/environmental_impact are untouched
+        — these are judgment calls (regulatory, economic, qualitative),
+        not physics, and stay on the categorical lookup-table heuristics
+        below (actuator's thrust_profile and energy_form and the combo's
+        infrastructure requirements).
 
         cost_efficiency additionally checks the domain's declared unit:
         "$/(kg·m)" (freight-style domains) isn't a rescaling of "$/m" — it's
@@ -1187,6 +1364,10 @@ class Pipeline:
         whose absolute mass tramples its own platform's declared ceiling
         can still produce perfectly plausible-looking per-kg ratios.
         """
+        if domain.metric_formulas:
+            return self._estimate_via_formulas(combo, domain)
+
+        metric_bounds = domain.metric_bounds
         metric_names = [mb.metric_name for mb in metric_bounds]
         units_by_name = {mb.metric_name: mb.unit for mb in metric_bounds}
         bounds_by_name = {mb.metric_name: mb for mb in metric_bounds}
@@ -1196,7 +1377,6 @@ class Pipeline:
         # drives the untouched blocks below).
         power_density = 0.0  # W/kg
         energy_density = 0.0  # J/kg
-        mass_total = 0.0  # kg, extensive — components share one vehicle
         thrust_profile: str | None = None
         energy_form: str | None = None
         infra_matches: list[float] = []
@@ -1206,8 +1386,6 @@ class Pipeline:
                     power_density = max(power_density, float(dep.value))
                 if dep.key == "energy_density" and dep.constraint_type == "provides":
                     energy_density = max(energy_density, float(dep.value))
-                if dep.key == "mass" and dep.constraint_type == "range_min":
-                    mass_total += float(dep.value)
                 if dep.key == "thrust_profile" and dep.constraint_type == "provides":
                     thrust_profile = dep.value
                 if dep.key == "energy_form" and dep.constraint_type == "requires":
@@ -1216,20 +1394,19 @@ class Pipeline:
                     match = INFRASTRUCTURE_AVAILABILITY.get((dep.key, dep.value))
                     if match is not None:
                         infra_matches.append(match)
-        mass = mass_total if mass_total > 0 else 100.0  # kg, default if undeclared
-        cargo_capacity_kg = mass * CARGO_KG_PER_STRUCTURAL_KG
 
         # ── platform/actuator/storage-specific extraction, for
-        # power_density / range_fuel / cost_efficiency only ──────────────
+        # power_density / range_fuel / cost_efficiency / cargo_capacity
+        # only ─────────────────────────────────────────────────────────
         ctx = self._physics_context(combo, bounds_by_name)
         feasible = True
         if ctx is not None:
-            actuator_mass, storage_mass, power_mass, denom_offset, platform_mass, feasible = self._decide_masses(
+            actuator_mass, storage_mass, platform_mass, feasible = self._decide_masses(
-                ctx, bounds_by_name, units_by_name, cargo_capacity_kg
+                ctx, bounds_by_name, units_by_name
             )
             raw.update(self._raw_physics_from_masses(
-                ctx, actuator_mass, storage_mass, power_mass, denom_offset,
+                ctx, actuator_mass, storage_mass,
-                bounds_by_name, units_by_name, cargo_capacity_kg,
+                bounds_by_name, units_by_name,
                 platform_mass=platform_mass,
             ))
 
@@ -1251,11 +1428,8 @@ class Pipeline:
         if "range_degradation" in raw:
             raw["range_degradation"] = 365 * 86400
 
-        if "cargo_capacity" in raw:
+        # cargo_capacity / cargo_capacity_kg are set by _raw_physics_from_masses
-            raw["cargo_capacity"] = cargo_capacity_kg
+        # above, from the actual build mass -- not recomputed here.
-
-        if "cargo_capacity_kg" in raw:
-            raw["cargo_capacity_kg"] = mass * 0.3
 
         if "environmental_impact" in raw:
             raw["environmental_impact"] = max(0.0, power_density * 2e-7)
@@ -1265,6 +1439,169 @@ class Pipeline:
 
         return raw, feasible
 
+    def _estimate_via_formulas(
+        self, combo: Combination, domain: Domain
+    ) -> tuple[dict[str, float], bool]:
+        """Generic estimator for domains that declare their own
+        metric_formulas/free_variables instead of matching the
+        platform/actuator/energy_storage shape _estimate_physics assumes.
+
+        Every formula (free-variable floor/ceiling, and each metric) is
+        evaluated through the same safe engine.formula evaluator, given two
+        kinds of names to resolve: dep(key, constraint_type="provides",
+        agg=None) -- a declared entity property, aggregated across the
+        combo's entities via constraint_resolver.aggregate_dependency_value
+        (the exact sum-vs-max rule ConstraintResolver itself already
+        applies to mass/footprint) -- and any free variable already chosen
+        earlier in its declared sort_order, as a bare name. Free variables
+        are resolved by _search_free_variables, a generalization of
+        _search_best_allocation/_decide_masses's hand-nested 2-3-variable
+        search to an arbitrary domain-declared list; with zero free
+        variables (most new domains -- nothing needs sizing) it degenerates
+        to a single direct formula evaluation, no search at all.
+
+        An invalid formula (bad syntax, unknown name, division by zero for
+        this combo's values) degrades that one metric to 0.0 / that one free
+        variable to being skipped, rather than failing the whole combo --
+        symmetric with how a metric_bounds entry with no matching raw value
+        already defaults to 0.0 elsewhere in this pipeline.
+        """
+        from physcom.engine.constraint_resolver import aggregate_dependency_value
+        from physcom.engine.formula import FormulaError, compile_formula, evaluate_formula
+
+        def dep(key: str, constraint_type: str = "provides", agg: str | None = None) -> float:
+            key_aggregation = None if agg is None else {key: agg}
+            value = aggregate_dependency_value(combo, key, constraint_type, key_aggregation)
+            return value if value is not None else 0.0
+
+        functions = {"dep": dep}
+
+        compiled_metrics: dict[str, object] = {}
+        for mf in domain.metric_formulas:
+            try:
+                compiled_metrics[mf.metric_name] = compile_formula(mf.formula)
+            except FormulaError:
+                compiled_metrics[mf.metric_name] = None
+
+        compiled_vars: list[tuple[str, object, object]] = []
+        for fv in sorted(domain.free_variables, key=lambda v: v.sort_order):
+            try:
+                compiled_vars.append((
+                    fv.name, compile_formula(fv.floor_formula), compile_formula(fv.ceiling_formula),
+                ))
+            except FormulaError:
+                continue
+
+        bounds_by_name = {mb.metric_name: mb for mb in domain.metric_bounds}
+
+        def evaluate_metrics(resolved: dict[str, float]) -> dict[str, float]:
+            raw: dict[str, float] = {}
+            for name, compiled in compiled_metrics.items():
+                if compiled is None:
+                    raw[name] = 0.0
+                    continue
+                try:
+                    raw[name] = evaluate_formula(compiled, resolved, functions)
+                except FormulaError:
+                    raw[name] = 0.0
+            return raw
+
+        def objective(resolved: dict[str, float]) -> float:
+            raw = evaluate_metrics(resolved)
+            scores, weights = [], []
+            for mb in bounds_by_name.values():
+                val = raw.get(mb.metric_name)
+                if val is None:
+                    continue
+                n = normalize(val, mb.norm_min, mb.norm_max)
+                if mb.lower_is_better:
+                    n = 1.0 - n
+                scores.append(n)
+                weights.append(mb.weight)
+            return composite_score(scores, weights)
+
+        resolved, _score, feasible = self._search_free_variables(
+            compiled_vars, {}, functions, objective
+        )
+        return evaluate_metrics(resolved), feasible
+
+    def _search_free_variables(
+        self,
+        var_specs: list[tuple[str, object, object]],
+        resolved: dict[str, float],
+        functions: dict,
+        objective,
+        grid: int = 10,
+        rounds: int = 5,
+    ) -> tuple[dict[str, float], float, bool]:
+        """Recursive generalization of _search_best_allocation/_decide_masses's
+        hand-nested outer-platform/inner-(actuator,storage) search to an
+        arbitrary ordered list of free variables: resolve var_specs[0]'s
+        floor/ceiling against what's already been resolved plus dep(...),
+        coarse-to-fine grid-search it, and recurse into the rest for each
+        candidate -- exactly the same nesting the hardcoded 3-variable
+        search already does by hand, just generalized to N. grid/rounds
+        shrink with remaining depth (cost is grid**depth * rounds**depth)
+        for the same reason the hardcoded search's outer loop already uses
+        a cheaper inner search per candidate; tuned for a handful of free
+        variables (transport's own 3 is the reference point), not domains
+        declaring dozens. Returns (resolved, score, feasible); feasible is
+        False only when no candidate at some depth has floor <= ceiling.
+        """
+        from physcom.engine.formula import FormulaError, evaluate_formula
+
+        if not var_specs:
+            return dict(resolved), objective(resolved), True
+
+        name, floor_f, ceiling_f = var_specs[0]
+        rest = var_specs[1:]
+        try:
+            floor = evaluate_formula(floor_f, resolved, functions)
+            ceiling = evaluate_formula(ceiling_f, resolved, functions)
+        except FormulaError:
+            result = dict(resolved)
+            result[name] = 0.0
+            return result, -1.0, False
+
+        if ceiling < floor:
+            result = dict(resolved)
+            result[name] = floor
+            return result, -1.0, False
+
+        inner_grid = max(4, grid - 2 * len(rest))
+        inner_rounds = max(2, rounds - len(rest))
+
+        if ceiling == floor:
+            candidate = dict(resolved)
+            candidate[name] = floor
+            return self._search_free_variables(
+                rest, candidate, functions, objective, inner_grid, inner_rounds
+            )
+
+        win_lo, win_hi = floor, ceiling
+        best_resolved, best_score, best_feasible = None, -1.0, False
+        for _round in range(rounds):
+            for i in range(grid + 1):
+                v = win_lo + (win_hi - win_lo) * i / grid
+                candidate = dict(resolved)
+                candidate[name] = v
+                r_resolved, r_score, r_feasible = self._search_free_variables(
+                    rest, candidate, functions, objective, inner_grid, inner_rounds
+                )
+                if r_score > best_score:
+                    best_resolved, best_score, best_feasible = r_resolved, r_score, r_feasible
+            if best_resolved is None:
+                break
+            span = max((win_hi - win_lo) / grid * 2, 1e-9)
+            center = best_resolved[name]
+            win_lo, win_hi = max(floor, center - span), min(ceiling, center + span)
+
+        if best_resolved is None:
+            result = dict(resolved)
+            result[name] = floor
+            return result, -1.0, False
+        return best_resolved, best_score, best_feasible
+
     def evaluate_allocation(
         self,
         combo: Combination,
@@ -1282,31 +1619,28 @@ class Pipeline:
         is ever persisted.
 
         Any mass left as None defaults to what the real requirement-based
-        solve already picked (see _decide_masses / _stub_estimate), so a
+        solve already picked (see _decide_masses / _estimate_physics), so a
         slider opens on today's actual build, not an arbitrary point.
         Explicit values are floor-clamped to each component's own declared
         minimum (platform is also ceiling-clamped to its declared max) --
         never silently allowed below what pass 1 would have rejected.
 
-        Returns None for combos with no free actuator mass to explore
+        Returns None only for combos with no declared platform mass
-        (biological actuators, radiation-pressure sails -- see
+        ceiling to bound a weight-class slider (e.g. Spaceship). Every
-        _stub_estimate's module note) or with no declared platform mass
+        actuator type gets sliders, including biological (rider/operator
-        ceiling to bound a weight-class slider.
+        mass, see BIOLOGICAL_OPERATOR_MASS_KG) and radiation-pressure
+        (effective sail mass derived from footprint, see
+        SAIL_AREAL_DENSITY_KG_PER_M2) -- both are real, budget-competing
+        variables like any mechanical actuator's mass.
         """
         bounds_by_name = {mb.metric_name: mb for mb in domain.metric_bounds}
         units_by_name = {mb.metric_name: mb.unit for mb in domain.metric_bounds}
         ctx = self._physics_context(combo, bounds_by_name)
         if ctx is None or ctx.p_max is None:
             return None
-        if (
-            ctx.actuator_energy_form in BIOLOGICAL_OPERATOR_MASS_KG
-            or ctx.actuator_energy_form == "radiation_pressure"
-        ):
-            return None
 
-        cargo_capacity_kg = (ctx.p_min + ctx.a_min + ctx.s_min) * CARGO_KG_PER_STRUCTURAL_KG
+        default_actuator, default_storage, default_platform, _feasible = self._decide_masses(
-        default_actuator, default_storage, _power_mass, _denom_offset, default_platform, _feasible = self._decide_masses(
+            ctx, bounds_by_name, units_by_name
-            ctx, bounds_by_name, units_by_name, cargo_capacity_kg
         )
         p_mass = default_platform if platform_mass is None else platform_mass
         a_mass = default_actuator if actuator_mass is None else actuator_mass
@@ -1317,8 +1651,8 @@ class Pipeline:
         s_mass = max(ctx.s_min, s_mass)
 
         raw = self._raw_physics_from_masses(
-            ctx, a_mass, s_mass, a_mass, 0.0,
+            ctx, a_mass, s_mass,
-            bounds_by_name, units_by_name, cargo_capacity_kg,
+            bounds_by_name, units_by_name,
             platform_mass=p_mass,
         )
 

src/physcom/models/domain.py

@@ -26,6 +26,32 @@ class DomainConstraint:
     allowed_values: list[str] = field(default_factory=list)  # e.g. ["ground", "air"]
 
 
+@dataclass
+class FreeVariable:
+    """A domain-declared quantity pass 2's estimator searches to maximize
+    the composite score (see Pipeline._estimate_via_formulas), e.g.
+    "actuator_mass". floor_formula/ceiling_formula are evaluated per-combo
+    and may reference dep(...) and any free variable declared at a lower
+    sort_order (mirrors the outer/inner nesting the built-in transport
+    physics model already does by hand)."""
+
+    name: str
+    floor_formula: str
+    ceiling_formula: str
+    sort_order: int = 0
+    id: int | None = None
+
+
+@dataclass
+class MetricFormula:
+    """A domain-declared formula computing one metric's raw value, evaluated
+    against dep(...) lookups and the domain's resolved free variables."""
+
+    metric_name: str
+    formula: str
+    id: int | None = None
+
+
 @dataclass
 class Domain:
     """A context frame that defines what 'good' means (e.g., urban_commuting)."""
@@ -34,4 +60,6 @@ class Domain:
     description: str = ""
     metric_bounds: list[MetricBound] = field(default_factory=list)
     constraints: list[DomainConstraint] = field(default_factory=list)
+    free_variables: list[FreeVariable] = field(default_factory=list)
+    metric_formulas: list[MetricFormula] = field(default_factory=list)
     id: int | None = None

src/physcom/seed/transport_example.py

@@ -95,6 +95,7 @@ WATER_PLATFORMS: list[Entity] = [
             Dependency("environment", "gravity", "true", None, "provides"),
             Dependency("physical", "footprint", "200", "m²", "range_max"),
             Dependency("physical", "footprint", "20", "m²", "range_min"),
+            Dependency("physical", "mass", "20000000", "kg", "range_max"),
             Dependency("physical", "mass", "10000", "kg", "range_min"),
             Dependency("environment", "medium", "water", None, "requires"),
             Dependency("physical", "energy_density", "720000", "J/kg", "range_min"),
@@ -198,6 +199,20 @@ MULTI_PLATFORMS: list[Entity] = [
             Dependency("physical", "footprint", "5", "m²", "range_min"),
             Dependency("physical", "mass", "10000", "kg", "range_max"),
             Dependency("physical", "mass", "1500", "kg", "range_min"),
+            # No requires here previously -- vacuously satisfied every
+            # domain's medium DomainConstraint (check_domain_constraints
+            # only flags a violation when an entity DECLARES a requires
+            # for the constrained key), including space-only
+            # interplanetary_travel. The current requires/domain-constraint
+            # model only supports one value per key -- there's no OR
+            # mechanism for "ground or water" -- so this picks ground
+            # (its primary, most-common domain) rather than leaving it
+            # undeclared. Real tradeoff: it can no longer participate in
+            # maritime_shipping (water-only) either, losing the water half
+            # of "amphibious." Closes the vacuous-pass hole; genuine
+            # multi-medium support would need OR semantics added to
+            # check_domain_constraints, a separate, bigger change.
+            Dependency("environment", "medium", "ground", None, "requires"),
         ],
     ),
 ]
@@ -741,9 +756,16 @@ URBAN_COMMUTING = Domain(
         # this project hasn't done -- not scored anywhere for now rather than
         # pretend a quick formula or an equally uninformed LLM guess settles it.
         # Weights renormalized to sum to 1.0 across the remaining metrics.
-        MetricBound("power_density", weight=0.4167, norm_min=1, norm_max=2000, unit="W/kg"),
+        # speed and cargo_capacity_kg added -- a commute's actual travel
-        MetricBound("cost_efficiency", weight=0.4167, norm_min=1e-5, norm_max=2e-3, unit="$/m", lower_is_better=True),
+        # time and whether the vehicle can carry groceries/passengers/gear
-        MetricBound("range_fuel", weight=0.1666, norm_min=5000, norm_max=500000, unit="m"),
+        # both matter as much as raw power_density did on their own; speed
+        # is a genuine build OUTPUT (see _raw_physics_from_masses), not a
+        # platform-declared constant.
+        MetricBound("power_density", weight=0.25, norm_min=1, norm_max=2000, unit="W/kg"),
+        MetricBound("cost_efficiency", weight=0.25, norm_min=1e-5, norm_max=2e-3, unit="$/m", lower_is_better=True),
+        MetricBound("speed", weight=0.25, norm_min=2, norm_max=30, unit="m/s"),
+        MetricBound("range_fuel", weight=0.10, norm_min=5000, norm_max=500000, unit="m"),
+        MetricBound("cargo_capacity_kg", weight=0.15, norm_min=1, norm_max=500, unit="kg"),
     ],
     constraints=[DomainConstraint("medium", ["ground", "air"])],
 )

src/physcom/snapshot.py

@@ -6,7 +6,7 @@ from datetime import datetime, timezone
 
 from physcom.db.repository import Repository
 from physcom.models.entity import Entity, Dependency
-from physcom.models.domain import Domain, DomainConstraint, MetricBound
+from physcom.models.domain import Domain, DomainConstraint, FreeVariable, MetricBound, MetricFormula
 from physcom.models.combination import Combination
 
 
@@ -70,11 +70,27 @@ def export_snapshot(repo: Repository) -> dict:
                 "key": dc.key,
                 "allowed_values": dc.allowed_values,
             })
+        fvs = []
+        for fv in d.free_variables:
+            fvs.append({
+                "name": fv.name,
+                "sort_order": fv.sort_order,
+                "floor_formula": fv.floor_formula,
+                "ceiling_formula": fv.ceiling_formula,
+            })
+        mfs = []
+        for mf in d.metric_formulas:
+            mfs.append({
+                "metric_name": mf.metric_name,
+                "formula": mf.formula,
+            })
         domain_list.append({
             "name": d.name,
             "description": d.description,
             "metric_bounds": mbs,
             "constraints": dcs,
+            "free_variables": fvs,
+            "metric_formulas": mfs,
         })
 
     # Export combinations
@@ -207,11 +223,26 @@ def import_snapshot(repo: Repository, data: dict, *, clear: bool = False) -> dic
                 )
                 for dc in d_data.get("constraints", [])
             ]
+            fvs = [
+                FreeVariable(
+                    name=fv["name"],
+                    floor_formula=fv["floor_formula"],
+                    ceiling_formula=fv["ceiling_formula"],
+                    sort_order=fv.get("sort_order", 0),
+                )
+                for fv in d_data.get("free_variables", [])
+            ]
+            mfs = [
+                MetricFormula(metric_name=mf["metric_name"], formula=mf["formula"])
+                for mf in d_data.get("metric_formulas", [])
+            ]
             domain = Domain(
                 name=d_data["name"],
                 description=d_data.get("description", ""),
                 metric_bounds=mbs,
                 constraints=dcs,
+                free_variables=fvs,
+                metric_formulas=mfs,
             )
             repo.add_domain(domain)
         counts["domains"] += 1

src/physcom_web/routes/domains.py

@@ -4,7 +4,7 @@ from __future__ import annotations
 
 from flask import Blueprint, flash, redirect, render_template, request, url_for
 
-from physcom.models.domain import Domain, MetricBound
+from physcom.models.domain import Domain, FreeVariable, MetricBound, MetricFormula
 from physcom_web.app import get_repo
 
 bp = Blueprint("domains", __name__, url_prefix="/domains")
@@ -117,3 +117,96 @@ def metric_delete(domain_id: int, metric_id: int):
     flash("Metric removed.", "success")
     domain = repo.get_domain_by_id(domain_id)
     return render_template("domains/_metrics_table.html", domain=domain)
+
+
+# ── Free variable CRUD (HTMX partials) ───────────────────────
+
+
+@bp.route("/<int:domain_id>/free-vars/add", methods=["POST"])
+def free_var_add(domain_id: int):
+    repo = get_repo()
+    name = request.form["name"].strip()
+    floor_formula = request.form.get("floor_formula", "").strip()
+    ceiling_formula = request.form.get("ceiling_formula", "").strip()
+    try:
+        sort_order = int(request.form.get("sort_order", "0"))
+    except ValueError:
+        sort_order = 0
+    if not name or not floor_formula or not ceiling_formula:
+        flash("Name, floor formula, and ceiling formula are required.", "error")
+    else:
+        fv = FreeVariable(
+            name=name, floor_formula=floor_formula, ceiling_formula=ceiling_formula,
+            sort_order=sort_order,
+        )
+        repo.add_free_variable(domain_id, fv)
+        flash(f"Free variable '{name}' added.", "success")
+    domain = repo.get_domain_by_id(domain_id)
+    return render_template("domains/_free_vars_table.html", domain=domain)
+
+
+@bp.route("/<int:domain_id>/free-vars/<int:fv_id>/edit", methods=["POST"])
+def free_var_edit(domain_id: int, fv_id: int):
+    repo = get_repo()
+    try:
+        sort_order = int(request.form.get("sort_order", "0"))
+    except ValueError:
+        sort_order = 0
+    fv = FreeVariable(
+        name=request.form["name"].strip(),
+        floor_formula=request.form.get("floor_formula", "").strip(),
+        ceiling_formula=request.form.get("ceiling_formula", "").strip(),
+        sort_order=sort_order,
+    )
+    repo.update_free_variable(fv_id, fv)
+    flash("Free variable updated.", "success")
+    domain = repo.get_domain_by_id(domain_id)
+    return render_template("domains/_free_vars_table.html", domain=domain)
+
+
+@bp.route("/<int:domain_id>/free-vars/<int:fv_id>/delete", methods=["POST"])
+def free_var_delete(domain_id: int, fv_id: int):
+    repo = get_repo()
+    repo.delete_free_variable(fv_id)
+    flash("Free variable removed.", "success")
+    domain = repo.get_domain_by_id(domain_id)
+    return render_template("domains/_free_vars_table.html", domain=domain)
+
+
+# ── Metric formula CRUD (HTMX partials) ──────────────────────
+
+
+@bp.route("/<int:domain_id>/formulas/add", methods=["POST"])
+def formula_add(domain_id: int):
+    repo = get_repo()
+    metric_name = request.form["metric_name"].strip()
+    formula = request.form.get("formula", "").strip()
+    if not metric_name or not formula:
+        flash("Metric name and formula are required.", "error")
+    else:
+        repo.add_metric_formula(domain_id, MetricFormula(metric_name=metric_name, formula=formula))
+        flash(f"Formula for '{metric_name}' added.", "success")
+    domain = repo.get_domain_by_id(domain_id)
+    return render_template("domains/_formulas_table.html", domain=domain)
+
+
+@bp.route("/<int:domain_id>/formulas/<int:formula_id>/edit", methods=["POST"])
+def formula_edit(domain_id: int, formula_id: int):
+    repo = get_repo()
+    mf = MetricFormula(
+        metric_name=request.form["metric_name"].strip(),
+        formula=request.form.get("formula", "").strip(),
+    )
+    repo.update_metric_formula(formula_id, mf)
+    flash("Formula updated.", "success")
+    domain = repo.get_domain_by_id(domain_id)
+    return render_template("domains/_formulas_table.html", domain=domain)
+
+
+@bp.route("/<int:domain_id>/formulas/<int:formula_id>/delete", methods=["POST"])
+def formula_delete(domain_id: int, formula_id: int):
+    repo = get_repo()
+    repo.delete_metric_formula(formula_id)
+    flash("Formula removed.", "success")
+    domain = repo.get_domain_by_id(domain_id)
+    return render_template("domains/_formulas_table.html", domain=domain)

src/physcom_web/routes/pipeline.py

@@ -32,6 +32,8 @@ def _run_pipeline_in_background(
     from physcom.engine.scorer import Scorer
     from physcom.engine.pipeline import Pipeline
 
+    conn = None
+    repo = None
     try:
         conn = init_db(db_path)
         repo = Repository(conn)
@@ -58,6 +60,7 @@ def _run_pipeline_in_background(
             run_id=run_id,
         )
     except Exception as exc:
+        if repo is not None:
             try:
                 repo.update_pipeline_run(
                     run_id, status="failed",
@@ -65,7 +68,15 @@ def _run_pipeline_in_background(
                 )
             except Exception:
                 pass
+        else:
+            # Couldn't even open the DB to record the failure (bad
+            # PHYSCOM_DB path, locked/corrupt file) -- the pipeline_runs
+            # row will stay "pending" forever with no way to write an
+            # error_message to it, so at least don't let that swallow the
+            # real cause silently. Server logs are the only trace left.
+            print(f"pipeline run {run_id} failed before DB was reachable: {exc!r}")
     finally:
+        if conn is not None:
             try:
                 conn.close()
             except Exception:

src/physcom_web/templates/domains/_formulas_table.html

@@ -0,0 +1,56 @@
+<table id="formulas-table">
+  <thead>
+    <tr>
+      <th>Metric</th>
+      <th>Formula</th>
+      <th></th>
+    </tr>
+  </thead>
+  <tbody>
+    {% for mf in domain.metric_formulas %}
+    <tr>
+      <td>{{ mf.metric_name }}</td>
+      <td><code>{{ mf.formula }}</code></td>
+      <td class="actions">
+        <button class="btn btn-sm"
+                onclick="this.closest('tr').nextElementSibling.style.display='table-row'; this.closest('tr').style.display='none'">
+          Edit
+        </button>
+        <form method="post"
+              hx-post="{{ url_for('domains.formula_delete', domain_id=domain.id, formula_id=mf.id) }}"
+              hx-target="#formulas-section" hx-swap="innerHTML"
+              class="inline-form">
+          <button type="submit" class="btn btn-sm btn-danger">Del</button>
+        </form>
+      </td>
+    </tr>
+    <tr class="edit-row" style="display:none">
+      <form method="post"
+            hx-post="{{ url_for('domains.formula_edit', domain_id=domain.id, formula_id=mf.id) }}"
+            hx-target="#formulas-section" hx-swap="innerHTML">
+        <td><input name="metric_name" value="{{ mf.metric_name }}" required></td>
+        <td><input name="formula" value="{{ mf.formula }}" required></td>
+        <td>
+          <button type="submit" class="btn btn-sm btn-primary">Save</button>
+          <button type="button" class="btn btn-sm"
+                  onclick="this.closest('tr').style.display='none'; this.closest('tr').previousElementSibling.style.display=''">
+            Cancel
+          </button>
+        </td>
+      </form>
+    </tr>
+    {% endfor %}
+  </tbody>
+</table>
+
+<h3>Add Formula</h3>
+<form method="post"
+      hx-post="{{ url_for('domains.formula_add', domain_id=domain.id) }}"
+      hx-target="#formulas-section" hx-swap="innerHTML"
+      class="dep-add-form">
+  <div class="form-row">
+    <input name="metric_name" placeholder="metric name" required>
+    <input name="formula" placeholder='formula, e.g. draw_weight_chosen * 2' required>
+    <button type="submit" class="btn btn-primary">Add</button>
+  </div>
+</form>

src/physcom_web/templates/domains/_free_vars_table.html

@@ -0,0 +1,64 @@
+<table id="free-vars-table">
+  <thead>
+    <tr>
+      <th>Name</th>
+      <th>Order</th>
+      <th>Floor formula</th>
+      <th>Ceiling formula</th>
+      <th></th>
+    </tr>
+  </thead>
+  <tbody>
+    {% for fv in domain.free_variables %}
+    <tr>
+      <td>{{ fv.name }}</td>
+      <td>{{ fv.sort_order }}</td>
+      <td><code>{{ fv.floor_formula }}</code></td>
+      <td><code>{{ fv.ceiling_formula }}</code></td>
+      <td class="actions">
+        <button class="btn btn-sm"
+                onclick="this.closest('tr').nextElementSibling.style.display='table-row'; this.closest('tr').style.display='none'">
+          Edit
+        </button>
+        <form method="post"
+              hx-post="{{ url_for('domains.free_var_delete', domain_id=domain.id, fv_id=fv.id) }}"
+              hx-target="#free-vars-section" hx-swap="innerHTML"
+              class="inline-form">
+          <button type="submit" class="btn btn-sm btn-danger">Del</button>
+        </form>
+      </td>
+    </tr>
+    <tr class="edit-row" style="display:none">
+      <form method="post"
+            hx-post="{{ url_for('domains.free_var_edit', domain_id=domain.id, fv_id=fv.id) }}"
+            hx-target="#free-vars-section" hx-swap="innerHTML">
+        <td><input name="name" value="{{ fv.name }}" required></td>
+        <td><input name="sort_order" type="number" step="1" value="{{ fv.sort_order }}"></td>
+        <td><input name="floor_formula" value="{{ fv.floor_formula }}" required></td>
+        <td><input name="ceiling_formula" value="{{ fv.ceiling_formula }}" required></td>
+        <td>
+          <button type="submit" class="btn btn-sm btn-primary">Save</button>
+          <button type="button" class="btn btn-sm"
+                  onclick="this.closest('tr').style.display='none'; this.closest('tr').previousElementSibling.style.display=''">
+            Cancel
+          </button>
+        </td>
+      </form>
+    </tr>
+    {% endfor %}
+  </tbody>
+</table>
+
+<h3>Add Free Variable</h3>
+<form method="post"
+      hx-post="{{ url_for('domains.free_var_add', domain_id=domain.id) }}"
+      hx-target="#free-vars-section" hx-swap="innerHTML"
+      class="dep-add-form">
+  <div class="form-row">
+    <input name="name" placeholder="name, e.g. draw_weight_chosen" required>
+    <input name="sort_order" type="number" step="1" placeholder="order" value="0">
+    <input name="floor_formula" placeholder='floor formula, e.g. dep("draw_weight", "range_min")' required>
+    <input name="ceiling_formula" placeholder='ceiling formula, e.g. dep("draw_weight", "range_max")' required>
+    <button type="submit" class="btn btn-primary">Add</button>
+  </div>
+</form>

src/physcom_web/templates/domains/detail.html

@@ -33,4 +33,18 @@
 <div id="metrics-section">
   {% include "domains/_metrics_table.html" %}
 </div>
+
+<h2>Free Variables</h2>
+<p class="hint">Quantities pass 2's estimator searches to maximize this domain's composite score. Leave empty for domains where nothing needs sizing.</p>
+
+<div id="free-vars-section">
+  {% include "domains/_free_vars_table.html" %}
+</div>
+
+<h2>Metric Formulas</h2>
+<p class="hint">How each metric's raw value is computed from declared entity properties (via <code>dep(key, constraint_type="provides")</code>) and any free variable above. If this domain declares any formulas, pass 2 uses them instead of the built-in vehicle physics model.</p>
+
+<div id="formulas-section">
+  {% include "domains/_formulas_table.html" %}
+</div>
 {% endblock %}

src/physcom_web/templates/results/_explore_result.html

@@ -1,7 +1,6 @@
 {% if explore_result is none %}
 <p class="empty">No free mass allocation to explore for this combination — its
-actuator's mass isn't a design choice (a physiological or footprint-derived
+platform has no declared mass ceiling to bound the sliders.</p>
-quantity), or the platform has no declared mass ceiling to bound the sliders.</p>
 {% else %}
 {% set r = explore_result %}
 <div class="optimize-summary">
@@ -26,9 +25,11 @@ quantity), or the platform has no declared mass ceiling to bound the sliders.</p
 
 <div class="mass-bar-container" title="platform {{ '%.1f'|format(r.platform_mass) }}kg / actuator {{ '%.1f'|format(r.actuator_mass) }}kg / storage {{ '%.1f'|format(r.storage_mass) }}kg">
   {% set total = r.total_mass %}
+  {% if total > 0 %}
   <div class="mass-bar-seg mass-bar-platform" style="width: {{ (r.platform_mass / total * 100)|round(1) }}%"></div>
   <div class="mass-bar-seg mass-bar-actuator" style="width: {{ (r.actuator_mass / total * 100)|round(1) }}%"></div>
   <div class="mass-bar-seg mass-bar-storage" style="width: {{ (r.storage_mass / total * 100)|round(1) }}%"></div>
+  {% endif %}
 </div>
 <div class="mass-bar-legend">
   <span><span class="mass-swatch mass-bar-platform"></span>platform {{ "%.1f"|format(r.platform_mass) }}kg</span>

src/physcom_web/templates/results/detail.html

@@ -129,7 +129,7 @@
 </div>
 {% endif %}
 
-{% if explore_result is not none %}
+{% if scores %}
 <h2>Explore: Scale the Build</h2>
 <p class="subtitle">
   Purely exploratory — nothing here is saved. Drag a slider to pick a
@@ -141,6 +141,7 @@
   or merely functional one.
 </p>
 <div class="card">
+  {% if explore_result is not none %}
   {% set r = explore_result %}
   <form id="explore-form"
         hx-post="{{ url_for('results.explore', domain_name=domain.name, combo_id=combo.id) }}"
@@ -154,7 +155,7 @@
       <output id="out_platform_mass">{{ "%.1f"|format(r.platform_mass) }}kg</output>
     </div>
     <div class="weight-slider-row">
-      <label for="actuator_mass">actuator (motor size)</label>
+      <label for="actuator_mass">actuator (motor/collector size, or operator count/size for muscle power)</label>
       <input type="range" min="{{ r.actuator_min }}" max="{{ r.actuator_slider_max }}" step="0.1"
              id="actuator_mass" name="actuator_mass" value="{{ r.actuator_mass }}"
              oninput="document.getElementById('out_actuator_mass').textContent = (+this.value).toFixed(1) + 'kg'">
@@ -168,6 +169,7 @@
       <output id="out_storage_mass">{{ "%.1f"|format(r.storage_mass) }}kg</output>
     </div>
   </form>
+  {% endif %}
   <div id="explore-result">
     {% include "results/_explore_result.html" %}
   </div>

tests/test_formula.py

@@ -0,0 +1,108 @@
+"""Tests for the safe formula evaluator."""
+
+import math
+
+import pytest
+
+from physcom.engine.formula import FormulaError, compile_formula, evaluate_formula
+
+
+class TestArithmetic:
+    def test_constant(self):
+        assert evaluate_formula(compile_formula("42"), {}) == 42.0
+
+    def test_basic_ops(self):
+        assert evaluate_formula(compile_formula("2 + 3 * 4"), {}) == 14.0
+        assert evaluate_formula(compile_formula("(2 + 3) * 4"), {}) == 20.0
+        assert evaluate_formula(compile_formula("10 / 4"), {}) == 2.5
+        assert evaluate_formula(compile_formula("2 ** 3"), {}) == 8.0
+
+    def test_unary_minus(self):
+        assert evaluate_formula(compile_formula("-5 + 2"), {}) == -3.0
+
+    def test_variable_lookup(self):
+        result = evaluate_formula(compile_formula("mass * 2"), {"mass": 3.0})
+        assert result == 6.0
+
+    def test_unknown_variable_raises(self):
+        with pytest.raises(FormulaError):
+            evaluate_formula(compile_formula("unknown_var"), {})
+
+    def test_division_by_zero_raises_formula_error(self):
+        with pytest.raises(FormulaError):
+            evaluate_formula(compile_formula("1 / 0"), {})
+
+
+class TestFunctions:
+    def test_default_math_functions(self):
+        assert evaluate_formula(compile_formula("sqrt(16)"), {}) == 4.0
+        assert evaluate_formula(compile_formula("max(1, 2, 3)"), {}) == 3.0
+        assert evaluate_formula(compile_formula("min(1, 2, 3)"), {}) == 1.0
+        assert evaluate_formula(compile_formula("abs(-5)"), {}) == 5.0
+        assert evaluate_formula(compile_formula("exp(0)"), {}) == 1.0
+        assert math.isclose(evaluate_formula(compile_formula("log(exp(1))"), {}), 1.0)
+
+    def test_custom_injected_function(self):
+        formula = compile_formula('dep("power_density", "provides")')
+        result = evaluate_formula(
+            formula, {}, functions={"dep": lambda key, constraint_type: 99.0}
+        )
+        assert result == 99.0
+
+    def test_unknown_function_raises(self):
+        with pytest.raises(FormulaError):
+            evaluate_formula(compile_formula("unknown_fn(1)"), {})
+
+    def test_string_constant_passthrough_to_function(self):
+        formula = compile_formula('dep("mass")')
+        result = evaluate_formula(formula, {}, functions={"dep": lambda key: len(key)})
+        assert result == 4.0
+
+
+class TestSecurity:
+    @pytest.mark.parametrize("source", [
+        "__import__('os').system('echo hi')",
+        "().__class__",
+        "[1, 2, 3]",
+        "{1: 2}",
+        "{1, 2}",
+        "(x for x in [1])",
+        "lambda: 1",
+        "1 if True else 0",
+        "1 == 1",
+        "x.__class__",
+        "x[0]",
+        "(lambda: 1)()",
+        "1; 2",
+    ])
+    def test_disallowed_constructs_rejected(self, source):
+        with pytest.raises(FormulaError):
+            compile_formula(source)
+
+    def test_unregistered_function_name_never_executes(self):
+        """exec/eval/__import__ etc. parse as ordinary Call nodes -- the
+        actual guarantee is that no function name is callable unless it's
+        explicitly in DEFAULT_FUNCTIONS or caller-supplied, checked at
+        evaluate time, not that the bare name is rejected at compile time."""
+        with pytest.raises(FormulaError):
+            evaluate_formula(compile_formula("exec('1')"), {})
+
+    def test_dunder_name_rejected(self):
+        with pytest.raises(FormulaError):
+            compile_formula("__builtins__")
+
+    def test_indirect_call_rejected(self):
+        with pytest.raises(FormulaError):
+            compile_formula("(a + b)(1)")
+
+    def test_invalid_syntax_raises_formula_error(self):
+        with pytest.raises(FormulaError):
+            compile_formula("2 +")
+
+    def test_boolean_constant_rejected(self):
+        with pytest.raises(FormulaError):
+            compile_formula("True")
+
+    def test_large_exponent_overflows_cleanly_not_hangs(self):
+        with pytest.raises(FormulaError):
+            evaluate_formula(compile_formula("9 ** 9 ** 9 ** 9"), {})

tests/test_pipeline_formulas.py

@@ -0,0 +1,113 @@
+"""End-to-end test that pass 2 can estimate a non-transport domain entirely
+from domain-authored formulas, without touching the platform/actuator/
+energy_storage physics model in Pipeline._estimate_physics."""
+
+import pytest
+
+from physcom.engine.constraint_resolver import ConstraintResolver
+from physcom.engine.scorer import Scorer
+from physcom.engine.pipeline import Pipeline
+from physcom.models.domain import Domain, FreeVariable, MetricBound, MetricFormula
+from physcom.models.entity import Dependency, Entity
+
+
+def _build_archery_domain(repo):
+    repo.add_entity(Entity(
+        name="Recurve",
+        dimension="bow",
+        dependencies=[
+            Dependency("physical", "draw_weight", "20", None, "range_min"),
+            Dependency("physical", "draw_weight", "50", None, "range_max"),
+        ],
+    ))
+    repo.add_entity(Entity(
+        name="Carbon",
+        dimension="arrow",
+        dependencies=[
+            Dependency("physical", "arrow_mass", "0.02", None, "provides"),
+        ],
+    ))
+    return repo.add_domain(Domain(
+        name="archery_test",
+        metric_bounds=[
+            MetricBound("drawback_force", weight=0.6, norm_min=0, norm_max=100),
+            MetricBound("range", weight=0.4, norm_min=0, norm_max=300),
+        ],
+        free_variables=[
+            FreeVariable(
+                name="draw_weight_chosen",
+                floor_formula='dep("draw_weight", "range_min")',
+                ceiling_formula='dep("draw_weight", "range_max")',
+                sort_order=0,
+            ),
+        ],
+        metric_formulas=[
+            MetricFormula(metric_name="drawback_force", formula="draw_weight_chosen * 2"),
+            MetricFormula(
+                metric_name="range",
+                formula='draw_weight_chosen * 5 / dep("arrow_mass")',
+            ),
+        ],
+    ))
+
+
+def test_formula_domain_scores_without_platform_actuator_shape(repo):
+    domain = _build_archery_domain(repo)
+    resolver = ConstraintResolver()
+    scorer = Scorer(domain)
+    pipeline = Pipeline(repo, resolver, scorer)
+
+    result = pipeline.run(
+        domain, ["bow", "arrow"], score_threshold=0.01, passes=[1, 2, 3, 5],
+    )
+
+    assert result.total_generated == 1
+    assert result.pass1_failed == 0
+    assert result.pass2_estimated == 1
+    assert result.pass3_above_threshold == 1
+
+    combos = repo.list_combinations()
+    assert len(combos) == 1
+    combo = combos[0]
+    scores = {
+        s["metric_name"]: s["raw_value"]
+        for s in repo.get_combination_scores(combo.id, domain.id)
+    }
+    # Both metrics increase monotonically with draw_weight_chosen and nothing
+    # trades off against it, so the optimizer should push to the declared
+    # ceiling (50) -- confirms _search_free_variables is actually searching,
+    # not just evaluating at the floor.
+    assert scores["drawback_force"] == pytest.approx(100.0, rel=0.02)
+
+
+def test_formula_domain_zero_free_variables_direct_evaluation(repo):
+    """A domain with metric_formulas but no free_variables should evaluate
+    each formula once directly -- no search loop at all."""
+    repo.add_entity(Entity(
+        name="Recurve",
+        dimension="bow",
+        dependencies=[Dependency("physical", "draw_weight", "30", None, "provides")],
+    ))
+    repo.add_entity(Entity(name="Carbon", dimension="arrow"))
+    domain = repo.add_domain(Domain(
+        name="archery_direct_test",
+        metric_bounds=[MetricBound("drawback_force", weight=1.0, norm_min=0, norm_max=100)],
+        metric_formulas=[
+            MetricFormula(metric_name="drawback_force", formula='dep("draw_weight") * 2'),
+        ],
+    ))
+    resolver = ConstraintResolver()
+    scorer = Scorer(domain)
+    pipeline = Pipeline(repo, resolver, scorer)
+
+    result = pipeline.run(
+        domain, ["bow", "arrow"], score_threshold=0.01, passes=[1, 2, 3, 5],
+    )
+
+    assert result.pass2_estimated == 1
+    combos = repo.list_combinations()
+    scores = {
+        s["metric_name"]: s["raw_value"]
+        for s in repo.get_combination_scores(combos[0].id, domain.id)
+    }
+    assert scores["drawback_force"] == 60.0

tests/test_repository.py

@@ -1,7 +1,7 @@
 """Tests for the database repository."""
 
 from physcom.models.entity import Entity, Dependency
-from physcom.models.domain import Domain, MetricBound
+from physcom.models.domain import Domain, FreeVariable, MetricBound, MetricFormula
 
 
 def test_ensure_dimension(repo):
@@ -63,6 +63,62 @@ def test_add_and_get_domain(repo):
     assert loaded.metric_bounds[0].metric_name == "speed"
 
 
+def test_add_domain_with_free_variables_and_formulas(repo):
+    domain = Domain(
+        name="archery_test",
+        metric_bounds=[MetricBound("drawback_force", weight=1.0, norm_min=0, norm_max=500)],
+        free_variables=[
+            FreeVariable(
+                name="draw_weight",
+                floor_formula='dep("draw_weight", "range_min")',
+                ceiling_formula='dep("draw_weight", "range_max")',
+                sort_order=0,
+            ),
+        ],
+        metric_formulas=[
+            MetricFormula(metric_name="drawback_force", formula="draw_weight * 1.5"),
+        ],
+    )
+    saved = repo.add_domain(domain)
+    assert saved.id is not None
+
+    loaded = repo.get_domain("archery_test")
+    assert loaded is not None
+    assert len(loaded.free_variables) == 1
+    assert loaded.free_variables[0].name == "draw_weight"
+    assert loaded.free_variables[0].id is not None
+    assert len(loaded.metric_formulas) == 1
+    assert loaded.metric_formulas[0].formula == "draw_weight * 1.5"
+
+
+def test_free_variable_and_formula_crud(repo):
+    domain = repo.add_domain(Domain(name="crud_test"))
+
+    fv = repo.add_free_variable(
+        domain.id,
+        FreeVariable(name="x", floor_formula="0", ceiling_formula="100", sort_order=0),
+    )
+    mf = repo.add_metric_formula(
+        domain.id, MetricFormula(metric_name="m", formula="x * 2")
+    )
+
+    repo.update_free_variable(
+        fv.id, FreeVariable(name="x", floor_formula="1", ceiling_formula="200", sort_order=0)
+    )
+    repo.update_metric_formula(mf.id, MetricFormula(metric_name="m", formula="x * 3"))
+
+    loaded = repo.get_domain_by_id(domain.id)
+    assert loaded.free_variables[0].floor_formula == "1"
+    assert loaded.free_variables[0].ceiling_formula == "200"
+    assert loaded.metric_formulas[0].formula == "x * 3"
+
+    repo.delete_free_variable(fv.id)
+    repo.delete_metric_formula(mf.id)
+    loaded = repo.get_domain_by_id(domain.id)
+    assert loaded.free_variables == []
+    assert loaded.metric_formulas == []
+
+
 def test_combination_save_and_dedup(repo):
     e1 = repo.add_entity(Entity(name="A", dimension="platform"))
     e2 = repo.add_entity(Entity(name="B", dimension="actuator"))

tests/test_snapshot.py

@@ -7,7 +7,7 @@ import pytest
 from physcom.db.schema import init_db
 from physcom.db.repository import Repository
 from physcom.models.entity import Entity, Dependency
-from physcom.models.domain import Domain, DomainConstraint, MetricBound
+from physcom.models.domain import Domain, DomainConstraint, FreeVariable, MetricBound, MetricFormula
 from physcom.models.combination import Combination
 from physcom.snapshot import export_snapshot, import_snapshot
 
@@ -169,6 +169,46 @@ def test_import_with_combinations(seeded_repo, tmp_path):
     assert len(fresh_combos) == len(data["combinations"])
 
 
+def test_export_import_roundtrip_free_variables_and_formulas(repo, tmp_path):
+    domain = Domain(
+        name="archery_snapshot_test",
+        metric_bounds=[MetricBound("drawback_force", weight=1.0, norm_min=0, norm_max=100)],
+        free_variables=[
+            FreeVariable(
+                name="draw_weight_chosen",
+                floor_formula='dep("draw_weight", "range_min")',
+                ceiling_formula='dep("draw_weight", "range_max")',
+                sort_order=0,
+            ),
+        ],
+        metric_formulas=[
+            MetricFormula(metric_name="drawback_force", formula="draw_weight_chosen * 2"),
+        ],
+    )
+    repo.add_domain(domain)
+
+    data = export_snapshot(repo)
+    exported = next(d for d in data["domains"] if d["name"] == "archery_snapshot_test")
+    assert exported["free_variables"] == [{
+        "name": "draw_weight_chosen", "sort_order": 0,
+        "floor_formula": 'dep("draw_weight", "range_min")',
+        "ceiling_formula": 'dep("draw_weight", "range_max")',
+    }]
+    assert exported["metric_formulas"] == [
+        {"metric_name": "drawback_force", "formula": "draw_weight_chosen * 2"},
+    ]
+
+    conn = init_db(tmp_path / "fresh.db")
+    fresh = Repository(conn)
+    import_snapshot(fresh, data, clear=True)
+
+    loaded = fresh.get_domain("archery_snapshot_test")
+    assert len(loaded.free_variables) == 1
+    assert loaded.free_variables[0].floor_formula == 'dep("draw_weight", "range_min")'
+    assert len(loaded.metric_formulas) == 1
+    assert loaded.metric_formulas[0].formula == "draw_weight_chosen * 2"
+
+
 def test_import_merge_skips_existing_domain(repo):
     """Merge import skips domains that already exist."""
     domain = Domain(