#!/usr/bin/env python3
"""Parametric CAD and presentation renders for the HTMAA benchtop TVAC.

The model intentionally separates engineering geometry (CadQuery/STEP/STL)
from presentation geometry (OpenSCAD/PNG), while generating both from the
same PARAMS dictionary. Measure the purchased chamber, update PARAMS, rerun.

Outputs:
  output/cad/tvac_testbed_assembly.step
  output/cad/stl/*.stl
  output/renders/*.png
  output/render_source/tvac_scene.scad
"""

from __future__ import annotations

import argparse
import json
import math
from pathlib import Path

import cadquery as cq
from cadquery import exporters
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import numpy as np


ROOT = Path(__file__).resolve().parent
OUTPUT = ROOT / "output"

# All dimensions are millimetres. Defaults are faithful to the v2.1 manual,
# with unmeasured COTS dimensions kept conspicuously editable.
PARAMS = {
    "pot_inner_d": 230.0,
    "pot_wall": 1.8,
    "pot_height": 250.0,
    "pot_base": 3.0,
    "pot_rim_outer_d": 250.0,
    "pot_rim_h": 8.0,
    "lid_d": 280.0,
    "lid_t": 20.0,
    "gasket_mean_d": 238.0,
    "gasket_w": 6.0,
    "gasket_depth": 2.5,
    "port_radius": 85.0,
    "pump_bore_d": 18.0,
    "gauge_bore_d": 15.0,
    "db25_cutout_x": 36.0,
    "db25_cutout_y": 14.0,
    "membrane_pocket_d": 72.0,
    "membrane_t": 4.0,
    "deck_d": 188.0,
    "deck_t": 6.0,
    "deck_z": 42.0,
    "post_d": 8.0,
    "post_radius": 77.0,
    "puck_size": 50.0,
    "puck_t": 2.0,
    "puck_center_x": 44.0,
    "cold_puck_z": 82.0,
    "puck_gap": 20.0,
    "occluder_t": 1.5,
    "shaft_d": 6.0,
    "magnet_d": 8.0,
    "magnet_t": 3.0,
    "rotor_d": 48.0,
    "rotor_t": 8.0,
    "coupling_clearance": 0.5,
    "deck_foot_h": 4.0,
    "cold_standoff_d": 5.0,
    "hot_support_post_d": 3.0,
    "hot_support_arm_d": 3.0,
    "occluder_swing_deg": 90.0,
    "occluder_swing_direction": -1.0,
    "occluder_clearance": 2.0,
    "occluder_hub_d": 16.0,
    "occluder_hub_h": 5.0,
}


def polar_xy(radius: float, angle_deg: float) -> tuple[float, float]:
    a = math.radians(angle_deg)
    return radius * math.cos(a), radius * math.sin(a)


def pot(p: dict) -> cq.Workplane:
    ro = p["pot_inner_d"] / 2 + p["pot_wall"]
    ri = p["pot_inner_d"] / 2
    h = p["pot_height"]
    body = cq.Workplane("XY").circle(ro).extrude(h)
    cavity = (
        cq.Workplane("XY")
        .workplane(offset=p["pot_base"])
        .circle(ri)
        .extrude(h + 2)
    )
    shell = body.cut(cavity)
    rim = (
        cq.Workplane("XY")
        .workplane(offset=h - p["pot_rim_h"])
        .circle(p["pot_rim_outer_d"] / 2)
        .circle(ri - 1.0)
        .extrude(p["pot_rim_h"])
    )
    foot = (
        cq.Workplane("XY")
        .circle(ro + 3)
        .circle(ro - 3)
        .extrude(4)
    )
    return shell.union(rim).union(foot)


def lid(p: dict) -> cq.Workplane:
    t = p["lid_t"]
    z0 = p["pot_height"]
    part = cq.Workplane("XY").workplane(offset=z0).circle(p["lid_d"] / 2).extrude(t)

    # DB25, pump, and gauge bores sit on the specified 85 mm circle.
    dbx, dby = polar_xy(p["port_radius"], 150)
    px, py = polar_xy(p["port_radius"], 30)
    gx, gy = polar_xy(p["port_radius"], 270)
    db_cut = (
        cq.Workplane("XY").workplane(offset=z0 - 1)
        .center(dbx, dby).rect(p["db25_cutout_x"], p["db25_cutout_y"])
        .extrude(t + 2)
    )
    pump_cut = (
        cq.Workplane("XY").workplane(offset=z0 - 1)
        .center(px, py).circle(p["pump_bore_d"] / 2).extrude(t + 2)
    )
    gauge_cut = (
        cq.Workplane("XY").workplane(offset=z0 - 1)
        .center(gx, gy).circle(p["gauge_bore_d"] / 2).extrude(t + 2)
    )

    # Top pocket leaves a 4 mm pressure membrane.
    pocket = (
        cq.Workplane("XY").workplane(offset=z0 + p["membrane_t"])
        .circle(p["membrane_pocket_d"] / 2)
        .extrude(t - p["membrane_t"] + 1)
    )

    # Underside gasket groove.
    rmean = p["gasket_mean_d"] / 2
    groove = (
        cq.Workplane("XY").workplane(offset=z0 - 0.1)
        .circle(rmean + p["gasket_w"] / 2)
        .circle(rmean - p["gasket_w"] / 2)
        .extrude(p["gasket_depth"] + 0.1)
    )
    return part.cut(db_cut).cut(pump_cut).cut(gauge_cut).cut(pocket).cut(groove)


def deck(p: dict) -> cq.Workplane:
    z = p["deck_z"]
    part = cq.Workplane("XY").workplane(offset=z).circle(p["deck_d"] / 2).extrude(p["deck_t"])
    part = part.cut(
        cq.Workplane("XY").workplane(offset=z - 1).circle((p["shaft_d"] + 1.2) / 2).extrude(p["deck_t"] + 2)
    )
    for a in (90, 210, 330):
        x, y = polar_xy(p["post_radius"], a)
        hole = (
            cq.Workplane("XY").workplane(offset=z - 1)
            .center(x, y).circle(2.2).extrude(p["deck_t"] + 2)
        )
        part = part.cut(hole)
    return part


def post(p: dict) -> cq.Workplane:
    z0 = p["pot_base"] + p["deck_foot_h"]
    return (
        cq.Workplane("XY")
        .workplane(offset=z0)
        .circle(p["post_d"] / 2)
        .extrude(p["deck_z"] - z0)
    )


def rubber_foot(p: dict) -> cq.Workplane:
    return (
        cq.Workplane("XY")
        .workplane(offset=p["pot_base"])
        .circle(7.0)
        .extrude(p["deck_foot_h"])
    )


def puck(p: dict, z: float) -> cq.Workplane:
    s = p["puck_size"]
    x = p["puck_center_x"]
    board = cq.Workplane("XY").workplane(offset=z).center(x, 0).rect(s, s).extrude(p["puck_t"])
    for dx in (-s / 2 + 5, s / 2 - 5):
        for dy in (-s / 2 + 5, s / 2 - 5):
            hole = (
                cq.Workplane("XY").workplane(offset=z - 1)
                .center(x + dx, dy).circle(1.6).extrude(p["puck_t"] + 2)
            )
            board = board.cut(hole)
    return board


def puck_standoff(p: dict) -> cq.Workplane:
    z0 = p["deck_z"] + p["deck_t"]
    h = p["cold_puck_z"] - z0
    return cq.Workplane("XY").workplane(offset=z0).circle(p["cold_standoff_d"] / 2).extrude(h)


def hot_support_layout(p: dict) -> list[tuple[tuple[float, float], tuple[float, float]]]:
    """Deck-post locations and their corresponding hot-puck mounting holes.

    All three vertical posts are outside the clockwise 0 to -90 degree
    occluder sweep. Short nylon arms connect them to three board holes.
    """
    cx = p["puck_center_x"]
    return [
        ((cx - 20, 32), (cx - 20, 20)),
        ((cx + 36, 20), (cx + 20, 20)),
        ((cx + 36, -20), (cx + 20, -20)),
    ]


def assert_occluder_clearance(p: dict) -> None:
    radius = p["hot_support_post_d"] / 2 + p["occluder_clearance"]
    xmin, xmax = 5 - radius, 72 + radius
    ymin, ymax = -22 - radius, 22 + radius
    end = p["occluder_swing_direction"] * p["occluder_swing_deg"]
    for (px, py), _target in hot_support_layout(p):
        for angle in np.linspace(0, end, int(abs(end)) * 2 + 1):
            a = math.radians(-angle)
            lx = math.cos(a) * px - math.sin(a) * py
            ly = math.sin(a) * px + math.cos(a) * py
            if xmin <= lx <= xmax and ymin <= ly <= ymax:
                raise ValueError(
                    f"Hot support at {(px, py)} intersects the occluder sweep at {angle:.1f} deg"
                )


def occluder(p: dict) -> cq.Workplane:
    cold_top = p["cold_puck_z"] + p["puck_t"]
    hot_bottom = cold_top + p["puck_gap"]
    z = (cold_top + hot_bottom) / 2 - p["occluder_t"] / 2
    # A radial vane sweeps into the overlapping puck area from the center shaft.
    pts = [(5, -22), (68, -22), (72, -14), (72, 14), (68, 22), (5, 22)]
    vane = cq.Workplane("XY").workplane(offset=z).polyline(pts).close().extrude(p["occluder_t"])
    hub_z = (cold_top + hot_bottom) / 2 - p["occluder_hub_h"] / 2
    hub = (
        cq.Workplane("XY").workplane(offset=hub_z)
        .circle(p["occluder_hub_d"] / 2)
        .extrude(p["occluder_hub_h"])
    )
    bore = (
        cq.Workplane("XY").workplane(offset=hub_z - 1)
        .circle((p["shaft_d"] + 0.25) / 2)
        .extrude(p["occluder_hub_h"] + 2)
    )
    return vane.union(hub).cut(bore)


def shaft(p: dict) -> cq.Workplane:
    z0 = p["deck_z"]
    z1 = p["pot_height"] - p["coupling_clearance"]
    return cq.Workplane("XY").workplane(offset=z0).circle(p["shaft_d"] / 2).extrude(z1 - z0)


def rotor(p: dict, z: float, pocket_face: str) -> cq.Workplane:
    hub = cq.Workplane("XY").workplane(offset=z).circle(p["rotor_d"] / 2).extrude(p["rotor_t"])
    hub = hub.cut(
        cq.Workplane("XY").workplane(offset=z - 1).circle((p["shaft_d"] + 0.4) / 2).extrude(p["rotor_t"] + 2)
    )
    for a in (0, 90, 180, 270):
        x, y = polar_xy(15, a)
        pocket_z = z - 0.1 if pocket_face == "bottom" else z + p["rotor_t"] - p["magnet_t"]
        m = (
            cq.Workplane("XY").workplane(offset=pocket_z)
            .center(x, y).circle(p["magnet_d"] / 2).extrude(p["magnet_t"] + 0.2)
        )
        hub = hub.cut(m)
    return hub


def simple_box(x: float, y: float, z: float, origin: tuple[float, float, float]) -> cq.Workplane:
    return cq.Workplane("XY").workplane(offset=origin[2]).center(origin[0], origin[1]).rect(x, y).extrude(z)


def ring(outer_r: float, inner_r: float, h: float, z: float) -> cq.Workplane:
    return cq.Workplane("XY").workplane(offset=z).circle(outer_r).circle(inner_r).extrude(h)


def gasket_mold(p: dict) -> cq.Workplane:
    """Open annular tray for casting the silicone lid gasket."""
    mean_r = p["gasket_mean_d"] / 2
    inner_r = mean_r - p["gasket_w"] / 2
    outer_r = mean_r + p["gasket_w"] / 2
    tray = cq.Workplane("XY").circle(outer_r + 7).circle(inner_r - 7).extrude(6)
    cavity = (
        cq.Workplane("XY").workplane(offset=3)
        .circle(outer_r).circle(inner_r).extrude(4)
    )
    return tray.cut(cavity)


def cylinder_between(a: tuple[float, float, float], b: tuple[float, float, float], r: float) -> cq.Workplane:
    av = np.array(a, dtype=float)
    bv = np.array(b, dtype=float)
    vec = bv - av
    length = float(np.linalg.norm(vec))
    shape = cq.Workplane("XY").circle(r).extrude(length)
    direction = vec / length
    zaxis = np.array([0.0, 0.0, 1.0])
    axis = np.cross(zaxis, direction)
    axis_norm = float(np.linalg.norm(axis))
    dot = float(np.clip(np.dot(zaxis, direction), -1.0, 1.0))
    if axis_norm > 1e-9:
        axis /= axis_norm
        shape = shape.rotate((0, 0, 0), tuple(axis), math.degrees(math.acos(dot)))
    elif dot < 0:
        shape = shape.rotate((0, 0, 0), (1, 0, 0), 180)
    return shape.translate(tuple(av))


def build_parts(p: dict) -> dict[str, cq.Workplane]:
    assert_occluder_clearance(p)
    cold_z = p["cold_puck_z"]
    hot_z = cold_z + p["puck_t"] + p["puck_gap"]
    inner_rotor_z = p["pot_height"] - p["rotor_t"] - p["coupling_clearance"]
    outer_rotor_z = p["pot_height"] + p["membrane_t"] + p["coupling_clearance"]
    parts: dict[str, cq.Workplane] = {
        "pot": pot(p),
        "lid": lid(p),
        "deck": deck(p),
        "cold_puck": puck(p, cold_z),
        "hot_puck": puck(p, hot_z),
        "occluder": occluder(p),
        "shaft": shaft(p),
        "inner_rotor": rotor(p, inner_rotor_z, "top"),
        "outer_rotor": rotor(p, outer_rotor_z, "bottom"),
    }
    for i, a in enumerate((90, 210, 330), start=1):
        x, y = polar_xy(p["post_radius"], a)
        parts[f"deck_post_{i}"] = post(p).translate((x, y, 0))
        parts[f"deck_foot_{i}"] = rubber_foot(p).translate((x, y, 0))
    # Four puck supports, shown as one fused component.
    supports = None
    s = p["puck_size"]
    for dx in (-s / 2 + 5, s / 2 - 5):
        for dy in (-s / 2 + 5, s / 2 - 5):
            st = puck_standoff(p).translate((p["puck_center_x"] + dx, dy, 0))
            supports = st if supports is None else supports.union(st)
    parts["puck_standoffs"] = supports

    # The hot puck is supported independently from the deck. These three
    # slender nylon posts sit outside the commanded occluder sweep; short top
    # arms reach the board holes without crossing the measurement gap.
    deck_top = p["deck_z"] + p["deck_t"]
    support_z = hot_z - 2.0
    for i, (base_xy, target_xy) in enumerate(hot_support_layout(p), start=1):
        bx, by = base_xy
        tx, ty = target_xy
        parts[f"hot_support_post_{i}"] = cylinder_between(
            (bx, by, deck_top), (bx, by, support_z), p["hot_support_post_d"] / 2
        )
        parts[f"hot_support_arm_{i}"] = cylinder_between(
            (bx, by, support_z),
            (tx, ty, support_z),
            p["hot_support_arm_d"] / 2,
        )
        parts[f"hot_support_pin_{i}"] = cylinder_between(
            (tx, ty, support_z),
            (tx, ty, hot_z + p["puck_t"] + 2.0),
            1.4,
        )

    # Explicit N42 magnet solids show pocket direction and the true face gap.
    for label, rotor_z, magnet_z in (
        ("inner", inner_rotor_z, inner_rotor_z + p["rotor_t"] - p["magnet_t"]),
        ("outer", outer_rotor_z, outer_rotor_z),
    ):
        for i, a in enumerate((0, 90, 180, 270), start=1):
            mx, my = polar_xy(15, a)
            parts[f"{label}_magnet_{i}"] = (
                cq.Workplane("XY").workplane(offset=magnet_z)
                .center(mx, my).circle(p["magnet_d"] / 2).extrude(p["magnet_t"])
            )

    # Gasket and recognizable external hardware.
    rmean = p["gasket_mean_d"] / 2
    parts["gasket"] = ring(rmean + p["gasket_w"] / 2, rmean - p["gasket_w"] / 2, 2.2, p["pot_height"] + 0.2)
    dbx, dby = polar_xy(p["port_radius"], 150)
    px, py = polar_xy(p["port_radius"], 30)
    gx, gy = polar_xy(p["port_radius"], 270)
    parts["db25_flange"] = simple_box(48, 25, 4, (dbx, dby, p["pot_height"] + p["lid_t"]))
    parts["db25_body"] = simple_box(34, 13, 9, (dbx, dby, p["pot_height"] + p["lid_t"] + 4))
    parts["pump_hex"] = (
        cq.Workplane("XY").workplane(offset=p["pot_height"] + p["lid_t"])
        .center(px, py).polygon(6, 22).extrude(6)
    )
    parts["pump_barb"] = (
        cq.Workplane("XY").workplane(offset=p["pot_height"] + p["lid_t"] + 6)
        .center(px, py).circle(7).extrude(31)
    )
    parts["gauge_hex"] = (
        cq.Workplane("XY").workplane(offset=p["pot_height"] + p["lid_t"])
        .center(gx, gy).polygon(6, 19).extrude(7)
    )
    parts["gauge_body"] = simple_box(38, 25, 24, (gx, gy, p["pot_height"] + p["lid_t"] + 7))

    # External stepper and bridge bracket, dimensionally recognizable as NEMA17.
    bracket_z = p["pot_height"] + p["lid_t"] + 10
    parts["bracket_left"] = simple_box(6, 64, 48, (-35, 0, bracket_z))
    parts["bracket_right"] = simple_box(6, 64, 48, (35, 0, bracket_z))
    parts["bracket_top"] = simple_box(76, 64, 5, (0, 0, bracket_z + 43))
    parts["stepper_body"] = simple_box(42, 42, 40, (0, 0, bracket_z + 48))
    parts["stepper_face"] = cq.Workplane("XY").workplane(offset=bracket_z + 45).circle(11).extrude(3)
    parts["stepper_shaft"] = cq.Workplane("XY").workplane(offset=p["pot_height"] + p["membrane_t"] + 6).circle(2.5).extrude(52)

    # PCB traces, sensors and compact wiring harness make the instrument legible.
    for label, z in (("cold", cold_z), ("hot", hot_z)):
        for i, y in enumerate(range(-18, 19, 6)):
            parts[f"{label}_trace_{i}"] = simple_box(40, 1.2, 0.35, (p["puck_center_x"], y, z + p["puck_t"]))
        for i, tx in enumerate((-14, 14)):
            parts[f"{label}_sensor_{i}"] = simple_box(4, 4, 1.5, (p["puck_center_x"] + tx, -20, z + p["puck_t"] + 0.35))
    wire_paths = [
        ((dbx, dby, p["pot_height"] - 4), (-70, 55, 150), 1.5, "wire_power_1"),
        ((-70, 55, 150), (20, 28, 115), 1.5, "wire_power_2"),
        ((20, 28, 115), (p["puck_center_x"], 24, cold_z + 3), 1.5, "wire_power_3"),
        ((dbx, dby, p["pot_height"] - 6), (-62, 48, 145), 1.2, "wire_signal_1"),
        ((-62, 48, 145), (p["puck_center_x"], 24, hot_z + 3), 1.2, "wire_signal_2"),
    ]
    for a, b, radius, name in wire_paths:
        parts[name] = cylinder_between(a, b, radius)
    return parts


COLORS = {
    "pot": cq.Color(0.56, 0.62, 0.67, 1.0),
    "lid": cq.Color(0.55, 0.87, 1.0, 0.35),
    "deck": cq.Color(0.18, 0.22, 0.27, 1.0),
    "cold_puck": cq.Color(0.03, 0.33, 0.26, 1.0),
    "hot_puck": cq.Color(0.03, 0.33, 0.26, 1.0),
    "occluder": cq.Color(0.92, 0.58, 0.10, 1.0),
    "shaft": cq.Color(0.78, 0.80, 0.82, 1.0),
    "inner_rotor": cq.Color(0.12, 0.14, 0.18, 1.0),
    "outer_rotor": cq.Color(0.12, 0.14, 0.18, 1.0),
    "stepper_body": cq.Color(0.10, 0.12, 0.15, 1.0),
    "gasket": cq.Color(0.10, 0.10, 0.12, 1.0),
    "db25_flange": cq.Color(0.70, 0.72, 0.75, 1.0),
    "db25_body": cq.Color(0.05, 0.06, 0.08, 1.0),
    "pump_hex": cq.Color(0.68, 0.70, 0.73, 1.0),
    "pump_barb": cq.Color(0.68, 0.70, 0.73, 1.0),
    "gauge_hex": cq.Color(0.68, 0.70, 0.73, 1.0),
    "gauge_body": cq.Color(0.06, 0.07, 0.09, 1.0),
}


def color_rgba(name: str) -> tuple[float, float, float, float]:
    if name == "lid": return (0.33, 0.78, 0.95, 0.28)
    if name == "pot": return (0.58, 0.62, 0.67, 1.0)
    if name == "gasket": return (0.06, 0.07, 0.08, 1.0)
    if name in {"cold_puck", "hot_puck"}: return (0.015, 0.31, 0.22, 1.0)
    if "trace" in name: return (0.92, 0.42, 0.09, 1.0)
    if "cold_sensor" in name: return (0.08, 0.58, 0.92, 1.0)
    if "hot_sensor" in name: return (0.95, 0.16, 0.05, 1.0)
    if name == "occluder": return (1.0, 0.42, 0.04, 1.0)
    if name.startswith("wire_power"): return (0.93, 0.17, 0.07, 1.0)
    if name.startswith("wire_signal"): return (0.05, 0.44, 0.93, 1.0)
    if name.startswith("hot_support"): return (0.88, 0.88, 0.82, 1.0)
    if name.startswith("deck_foot"): return (0.05, 0.05, 0.06, 1.0)
    if "magnet" in name: return (0.74, 0.78, 0.83, 1.0)
    if name in {"deck", "inner_rotor", "outer_rotor"}: return (0.10, 0.13, 0.17, 1.0)
    if name.startswith("bracket") or name == "stepper_body": return (0.08, 0.10, 0.13, 1.0)
    if name == "db25_body" or name == "gauge_body": return (0.04, 0.05, 0.07, 1.0)
    return (0.66, 0.69, 0.72, 1.0)


def export_cad(p: dict) -> None:
    cad_dir = OUTPUT / "cad"
    stl_dir = cad_dir / "stl"
    stl_dir.mkdir(parents=True, exist_ok=True)
    parts = build_parts(p)
    assy = cq.Assembly(name="HTMAA_TVAC_Testbed")
    fabricated_assy = cq.Assembly(name="HTMAA_TVAC_Fabricated_Subsystem")
    fabricated_exact = {
        "lid", "deck", "cold_puck", "hot_puck", "occluder",
        "inner_rotor", "outer_rotor", "puck_standoffs", "gasket",
    }
    for name, shape in parts.items():
        rgba = color_rgba(name)
        color = COLORS.get(name, cq.Color(*rgba))
        assy.add(shape, name=name, color=color)
        exporters.export(shape, str(stl_dir / f"{name}.stl"), tolerance=0.08, angularTolerance=0.2)
        if name in fabricated_exact or name.startswith("hot_support"):
            fabricated_assy.add(shape, name=name, color=color)
    assy.save(str(cad_dir / "tvac_testbed_assembly.step"))
    fabricated_assy.add(
        gasket_mold(p).translate((p["lid_d"] / 2 + 150, 0, 0)),
        name="gasket_mold_print",
        color=cq.Color(0.55, 0.35, 0.72, 1.0),
    )
    fabricated_assy.save(str(cad_dir / "tvac_fabricated_parts.step"))
    exporters.export(
        gasket_mold(p), str(stl_dir / "gasket_mold_print.stl"),
        tolerance=0.08, angularTolerance=0.2,
    )
    inner_top = p["pot_height"] - p["coupling_clearance"]
    outer_bottom = p["pot_height"] + p["membrane_t"] + p["coupling_clearance"]
    hot_z = p["cold_puck_z"] + p["puck_t"] + p["puck_gap"]
    post_length_m = (hot_z - 2.0 - (p["deck_z"] + p["deck_t"])) / 1000
    post_area_m2 = math.pi * (p["hot_support_post_d"] / 2000) ** 2
    nylon_k = 0.25
    design_report = {
        "magnet_face_gap_mm": outer_bottom - inner_top,
        "membrane_thickness_mm": p["membrane_t"],
        "occluder_sweep_deg": [0, p["occluder_swing_direction"] * p["occluder_swing_deg"]],
        "occluder_collision_check": "passed",
        "hot_puck_support": "three independent deck-mounted nylon posts with short top arms",
        "estimated_hot_support_vertical_conductance_W_per_K": 3 * nylon_k * post_area_m2 / post_length_m,
        "estimated_hot_support_vertical_loss_at_30K_W": 3 * nylon_k * post_area_m2 / post_length_m * 30,
        "note": "Estimate excludes short horizontal arms and fastener contacts; measure final parasitic conductance.",
    }
    (cad_dir / "design_checks.json").write_text(json.dumps(design_report, indent=2) + "\n")
    (cad_dir / "parameters.json").write_text(json.dumps(p, indent=2) + "\n")


def scad_source(p: dict) -> str:
    # Presentation model. Fine PCB traces, wiring, fittings and brackets are
    # intentionally visualization-level details; pressure-boundary dimensions
    # come directly from PARAMS.
    return f'''// Generated by testbed_pot_cad.py; edit PARAMS in the Python source.
$fn = 96;
SCENE = is_undef(SCENE) ? "hero" : SCENE;

pot_id={p['pot_inner_d']}; pot_wall={p['pot_wall']}; pot_h={p['pot_height']};
pot_ro=pot_id/2+pot_wall; rim_ro={p['pot_rim_outer_d']/2};
lid_d={p['lid_d']}; lid_t={p['lid_t']}; deck_d={p['deck_d']}; deck_z={p['deck_z']};
puck={p['puck_size']}; puck_t={p['puck_t']}; puck_x={p['puck_center_x']};
cold_z={p['cold_puck_z']}; gap={p['puck_gap']}; shaft_d={p['shaft_d']};

module steel(c=[0.50,0.56,0.62,1]) {{ color(c) children(); }}
module acrylic() {{ color([0.48,0.86,1.0,0.42]) children(); }}
module pcb() {{ color([0.02,0.32,0.23,1]) children(); }}
module black() {{ color([0.055,0.065,0.08,1]) children(); }}
module orange() {{ color([1.0,0.47,0.06,1]) children(); }}
module copper() {{ color([0.92,0.46,0.13,1]) children(); }}
module silicone() {{ color([0.12,0.12,0.14,1]) children(); }}

module rounded_box(size=[10,10,4], r=1) {{
  hull() for (x=[-size[0]/2+r,size[0]/2-r], y=[-size[1]/2+r,size[1]/2-r])
    translate([x,y,0]) cylinder(r=r,h=size[2]);
}}

module pot(full=true) {{
  module raw_pot() {{
    steel() difference() {{
      union() {{
        cylinder(r=pot_ro,h=pot_h);
        translate([0,0,pot_h-8]) cylinder(r=rim_ro,h=8);
        cylinder(r=pot_ro+3,h=4);
      }}
      translate([0,0,3]) cylinder(r=pot_id/2,h=pot_h+10);
    }}
  }}
  if (full) raw_pot();
  else intersection() {{ raw_pot(); translate([-180,0,-5]) cube([360,180,280]); }}
}}

module gasket(z=pot_h+0.3) {{
  silicone() translate([0,0,z]) difference() {{
    cylinder(r={p['gasket_mean_d']/2+p['gasket_w']/2},h=2.2);
    translate([0,0,-1]) cylinder(r={p['gasket_mean_d']/2-p['gasket_w']/2},h=4.2);
  }}
}}

module lid(z=pot_h, show_fittings=true) {{
  acrylic() translate([0,0,z]) difference() {{
    cylinder(r=lid_d/2,h=lid_t);
    translate([0,0,{p['membrane_t']}]) cylinder(r={p['membrane_pocket_d']/2},h=lid_t+2);
    translate([{polar_xy(p['port_radius'],150)[0]:.4f},{polar_xy(p['port_radius'],150)[1]:.4f},-1]) cube([{p['db25_cutout_x']},{p['db25_cutout_y']},lid_t+2],center=true);
    translate([{polar_xy(p['port_radius'],30)[0]:.4f},{polar_xy(p['port_radius'],30)[1]:.4f},-1]) cylinder(r={p['pump_bore_d']/2},h=lid_t+2);
    translate([{polar_xy(p['port_radius'],270)[0]:.4f},{polar_xy(p['port_radius'],270)[1]:.4f},-1]) cylinder(r={p['gauge_bore_d']/2},h=lid_t+2);
  }}
  if (show_fittings) {{
    db25(z+lid_t);
    pump_fitting(z+lid_t);
    gauge_fitting(z+lid_t);
    coupling(z);
    stepper(z+lid_t+10);
  }}
}}

module db25(z) {{
  x={polar_xy(p['port_radius'],150)[0]:.4f}; y={polar_xy(p['port_radius'],150)[1]:.4f};
  translate([x,y,z]) rotate([0,0,150]) {{
    steel([0.72,0.74,0.76,1]) rounded_box([48,25,4],3);
    black() translate([0,0,3.5]) rounded_box([34,12,8],2);
    color([0.85,0.72,0.22,1]) for (r=[-3,-1,1,3], c=[-5:1:5])
      translate([c*2.4+(r%2)*0.6,r*1.3,11]) cylinder(r=.55,h=2);
  }}
}}

module hose_barb(r=8,h=34) {{
  steel([0.72,0.74,0.76,1]) union() {{
    cylinder(r=r+4,h=5,$fn=6);
    translate([0,0,5]) cylinder(r=r,h=h-5);
    for(z=[11:6:h-2]) translate([0,0,z]) cylinder(r1=r+1.8,r2=r,h=3);
  }}
}}

module pump_fitting(z) {{
  translate([{polar_xy(p['port_radius'],30)[0]:.4f},{polar_xy(p['port_radius'],30)[1]:.4f},z]) hose_barb(7,38);
}}

module gauge_fitting(z) {{
  translate([{polar_xy(p['port_radius'],270)[0]:.4f},{polar_xy(p['port_radius'],270)[1]:.4f},z]) {{
    hose_barb(5.5,18);
    translate([0,0,18]) steel([0.68,0.70,0.72,1]) cylinder(r=16,h=8,$fn=6);
    translate([0,0,26]) black() rounded_box([38,23,22],3);
    translate([0,-11.7,30]) color([0.15,0.72,0.90,1]) cube([25,1,11],center=true);
  }}
}}

module magnet_set(z) {{
  color([0.77,0.80,0.84,1]) for(a=[0,90,180,270])
    rotate([0,0,a]) translate([15,0,z]) cylinder(r={p['magnet_d']/2},h={p['magnet_t']});
}}

module coupling(lidz=pot_h) {{
  // Internal and external rotors straddle the 4 mm acrylic membrane.
  black() translate([0,0,lidz-13]) cylinder(r={p['rotor_d']/2},h={p['rotor_t']});
  magnet_set(lidz-12);
  black() translate([0,0,lidz+{p['membrane_t']}+1]) cylinder(r={p['rotor_d']/2},h={p['rotor_t']});
  magnet_set(lidz+{p['membrane_t']}+2);
}}

module stepper(z) {{
  // Bridge bracket leaves the membrane pocket visible.
  steel([0.16,0.19,0.23,1]) {{
    translate([-38,-32,z]) cube([6,64,48]);
    translate([32,-32,z]) cube([6,64,48]);
    translate([-38,-32,z+43]) cube([76,64,5]);
  }}
  black() translate([0,0,z+48]) rounded_box([42,42,40],3);
  steel([0.72,0.74,0.76,1]) translate([0,0,z+88]) cylinder(r=11,h=2);
  steel([0.72,0.74,0.76,1]) translate([0,0,z+34]) cylinder(r=2.5,h=18);
  color([0.28,0.30,0.34,1]) for(a=[45,135,225,315]) rotate([0,0,a]) translate([23,0,z+45]) cylinder(r=2,h=5);
}}

module deck(z=deck_z) {{
  steel([0.18,0.22,0.27,1]) translate([0,0,z]) difference() {{
    cylinder(r=deck_d/2,h={p['deck_t']});
    translate([0,0,-1]) cylinder(r=(shaft_d+1.2)/2,h={p['deck_t']+2});
  }}
  for(a=[90,210,330]) rotate([0,0,a]) translate([{p['post_radius']},0,5]) {{
    steel([0.34,0.36,0.39,1]) cylinder(r={p['post_d']/2},h=deck_z-5);
    black() translate([0,0,-4]) cylinder(r=7,h=5);
  }}
}}

module board(z, hot=false) {{
  translate([puck_x,0,z]) {{
    pcb() rounded_box([puck,puck,puck_t],2.5);
    // Copper serpentine reads clearly at website scale.
    copper() for(y=[-18:-6:18]) translate([-20,y,puck_t+.05]) cube([40,1.25,.28]);
    copper() for(y=[-18,-6,6,18]) translate([19,y+3,puck_t+.05]) cube([1.25,6,.28]);
    copper() for(y=[-12,0,12]) translate([-20,y+3,puck_t+.05]) cube([1.25,6,.28]);
    color(hot ? [0.96,0.18,0.08,1] : [0.12,0.70,0.92,1])
      for(x=[-14,14]) translate([x,-20,puck_t+.2]) cube([4,4,1.5],center=true);
    black() translate([0,18,puck_t+.3]) rounded_box([15,6,3],1);
    color([0.88,0.88,0.90,1]) for(x=[-20,20],y=[-20,20]) translate([x,y,-.5]) cylinder(r=1.6,h=4);
  }}
}}

module pucks() {{
  // Lower supports and the upper four-column tower make the load path legible.
  for(x=[puck_x-20,puck_x+20],y=[-20,20]) {{
    steel([0.72,0.74,0.76,1]) translate([x,y,deck_z+{p['deck_t']}]) cylinder(r=2.5,h=cold_z-deck_z-{p['deck_t']});
    steel([0.72,0.74,0.76,1]) translate([x,y,cold_z+puck_t]) cylinder(r=2.5,h=gap);
  }}
  board(cold_z,false);
  board(cold_z+puck_t+gap,true);
}}

module occluder() {{
  z=cold_z+puck_t+gap/2-{p['occluder_t']/2};
  orange() translate([0,0,z]) linear_extrude(height={p['occluder_t']})
    polygon([[5,-22],[68,-22],[72,-14],[72,14],[68,22],[5,22]]);
  steel([0.76,0.78,0.80,1]) translate([0,0,deck_z]) cylinder(r=shaft_d/2,h=pot_h-deck_z-5);
}}

module wire_segment(a,b,r=1.25,c=[0.9,0.2,0.1,1]) {{
  color(c) hull() {{ translate(a) sphere(r=r,$fn=18); translate(b) sphere(r=r,$fn=18); }}
}}

module harness() {{
  db=[{polar_xy(p['port_radius'],150)[0]:.4f},{polar_xy(p['port_radius'],150)[1]:.4f},pot_h-4];
  wire_segment(db,[-70,55,150],1.5,[0.95,0.28,0.10,1]);
  wire_segment([-70,55,150],[20,28,115],1.5,[0.95,0.28,0.10,1]);
  wire_segment([20,28,115],[puck_x,24,cold_z+3],1.5,[0.95,0.28,0.10,1]);
  wire_segment(db,[-62,48,145],1.25,[0.12,0.55,0.95,1]);
  wire_segment([-62,48,145],[puck_x,24,cold_z+puck_t+gap+3],1.25,[0.12,0.55,0.95,1]);
  black() translate([puck_x,27,cold_z+1]) rounded_box([18,8,6],1);
}}

module internals() {{ deck(); pucks(); occluder(); harness(); }}

module hero() {{
  pot(true); gasket(); lid();
}}

module cutaway() {{
  pot(false); gasket(); lid(); internals();
}}

module exploded() {{
  pot(false);
  translate([0,0,75]) internals();
  translate([0,0,150]) {{ gasket(); lid(); }}
}}

module detail() {{
  deck(); pucks(); occluder(); harness();
}}

module lid_top() {{
  gasket(); lid(pot_h,false);
}}

if (SCENE=="hero") hero();
else if (SCENE=="cutaway") cutaway();
else if (SCENE=="exploded") exploded();
else if (SCENE=="detail") detail();
else if (SCENE=="lid") lid_top();
'''


RENDERS = {
    "01_hero_isometric": ("hero", "390,-510,365,0,0,150"),
    "02_cutaway": ("cutaway", "385,-500,300,0,0,145"),
    "03_exploded": ("exploded", "500,-650,560,0,0,250"),
    "04_experiment_detail": ("detail", "245,-330,215,30,0,105"),
    "05_lid_top": ("lid", "0,-40,720,0,0,260"),
}


def tessellated(shape: cq.Workplane, tolerance: float = 1.0) -> np.ndarray:
    vertices, faces = shape.val().tessellate(tolerance, 0.35)
    xyz = np.array([[v.x, v.y, v.z] for v in vertices], dtype=float)
    return xyz[np.array(faces, dtype=int)]


def mesh_payload(shape: cq.Workplane, rgba: tuple[float, float, float, float]) -> tuple[np.ndarray, np.ndarray]:
    triangles = tessellated(shape)
    normals = np.cross(triangles[:, 1] - triangles[:, 0], triangles[:, 2] - triangles[:, 0])
    lengths = np.linalg.norm(normals, axis=1)
    normals[lengths > 0] /= lengths[lengths > 0, None]
    light = np.array([-0.35, -0.42, 0.84])
    intensity = np.clip(0.54 + 0.44 * (normals @ light), 0.32, 1.0)
    base = np.array(rgba[:3])
    facecolors = np.clip(base[None, :] * intensity[:, None] + 0.045, 0, 1)
    colors = np.column_stack([facecolors, np.full(len(facecolors), rgba[3])])
    return triangles, colors


def add_mesh(ax, shape: cq.Workplane, rgba: tuple[float, float, float, float], name: str) -> None:
    triangles, colors = mesh_payload(shape, rgba)
    collection = Poly3DCollection(
        triangles,
        facecolors=colors,
        edgecolors="none",
        linewidths=0,
        zsort="average",
    )
    ax.add_collection3d(collection)


def add_opaque_mesh_batch(ax, items: list[tuple[str, cq.Workplane]]) -> None:
    all_triangles = []
    all_colors = []
    for name, shape in items:
        triangles, colors = mesh_payload(shape, color_rgba(name))
        all_triangles.append(triangles)
        all_colors.append(colors)
    collection = Poly3DCollection(
        np.concatenate(all_triangles, axis=0),
        facecolors=np.concatenate(all_colors, axis=0),
        edgecolors="none",
        linewidths=0,
        zsort="average",
    )
    ax.add_collection3d(collection)


def cutaway_pot(p: dict) -> cq.Workplane:
    # Retain the rear half (positive Y); camera views from negative Y.
    clipping_halfspace = cq.Workplane("XY").box(500, 250, 450, centered=(True, False, False))
    return pot(p).intersect(clipping_halfspace)


def scene_parts(p: dict, scene: str) -> list[tuple[str, cq.Workplane]]:
    parts = build_parts(p)
    external = {
        "lid", "gasket", "outer_rotor", "db25_flange", "db25_body",
        "pump_hex", "pump_barb", "gauge_hex", "gauge_body",
        "bracket_left", "bracket_right", "bracket_top", "stepper_body",
        "stepper_face", "stepper_shaft",
    }
    external.update(name for name in parts if name.startswith("outer_magnet"))
    internal = set(parts) - external - {"pot"}
    if scene == "hero":
        names = ["pot"] + sorted(external)
        return [(name, parts[name]) for name in names]
    if scene == "cutaway":
        visible = [("pot", cutaway_pot(p))]
        visible.extend((name, parts[name]) for name in sorted(set(parts) - {"pot"}))
        return visible
    if scene == "exploded":
        visible = [("pot", cutaway_pot(p))]
        visible.extend((name, parts[name].translate((0, 0, 65))) for name in sorted(internal))
        visible.extend((name, parts[name].translate((0, 0, 150))) for name in sorted(external))
        return visible
    if scene == "detail":
        names = sorted(
            name for name in internal
            if name not in {"inner_rotor", "shaft"}
            and not name.startswith("inner_magnet")
        )
        visible = [(name, parts[name]) for name in names]
        detail_shaft = (
            cq.Workplane("XY").workplane(offset=p["deck_z"])
            .circle(p["shaft_d"] / 2).extrude(86)
        )
        visible.append(("shaft", detail_shaft))
        return visible
    if scene == "lid":
        return [("lid", parts["lid"]), ("gasket", parts["gasket"])]
    raise ValueError(scene)


VIEW_SETTINGS = {
    "hero": dict(elev=23, azim=-58, bounds=(-180, 180, -180, 180, -10, 405)),
    "cutaway": dict(elev=19, azim=-62, bounds=(-180, 180, -180, 180, -10, 405)),
    "exploded": dict(elev=18, azim=-57, bounds=(-190, 190, -190, 190, -10, 590)),
    "detail": dict(elev=30, azim=-63, bounds=(-110, 115, -110, 110, 30, 155)),
    "lid": dict(elev=88, azim=-90, bounds=(-155, 155, -155, 155, 244, 276)),
}


def render_mesh_images(p: dict, render_dir: Path) -> None:
    bg = "#f3f5f7"
    for filename, (scene, _camera) in RENDERS.items():
        cfg = VIEW_SETTINGS[scene]
        fig = plt.figure(figsize=(12.86, 10), dpi=140, facecolor=bg)
        ax = fig.add_subplot(111, projection="3d", facecolor=bg)
        items = scene_parts(p, scene)
        opaque = [(name, shape) for name, shape in items if color_rgba(name)[3] >= 1.0]
        transparent = [(name, shape) for name, shape in items if color_rgba(name)[3] < 1.0]
        if opaque:
            add_opaque_mesh_batch(ax, opaque)
        for name, shape in transparent:
            add_mesh(ax, shape, color_rgba(name), name)

        xmin, xmax, ymin, ymax, zmin, zmax = cfg["bounds"]
        ax.set_xlim(xmin, xmax)
        ax.set_ylim(ymin, ymax)
        ax.set_zlim(zmin, zmax)
        ax.set_box_aspect((xmax - xmin, ymax - ymin, zmax - zmin))
        ax.view_init(elev=cfg["elev"], azim=cfg["azim"], roll=0)
        ax.set_proj_type("persp", focal_length=1.15)
        ax.set_axis_off()
        fig.subplots_adjust(left=0, right=1, bottom=0, top=1)
        target = render_dir / f"{filename}.png"
        temporary = render_dir / f".{filename}.tmp"
        fig.savefig(temporary, format="png", facecolor=bg, bbox_inches="tight", pad_inches=0.02)
        plt.close(fig)
        temporary.replace(target)


def render_images(p: dict) -> None:
    render_dir = OUTPUT / "renders"
    render_dir.mkdir(parents=True, exist_ok=True)
    render_mesh_images(p, render_dir)


def main() -> None:
    parser = argparse.ArgumentParser()
    parser.add_argument("--cad-only", action="store_true")
    parser.add_argument("--renders-only", action="store_true")
    args = parser.parse_args()
    if not args.renders_only:
        export_cad(PARAMS)
    if not args.cad_only:
        render_images(PARAMS)
    print(f"Wrote TVAC CAD package to {OUTPUT}")


if __name__ == "__main__":
    main()