diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index c9df706..7ae1b71 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -65,14 +65,14 @@ jobs:
         if: ${{ matrix.python == '3.12' }}
         run: uvx ruff format --check .
 
-      - name: Tests (pytest; tolerate "no tests collected" while skeletal)
-        run: |
-          set +e
-          uv run --python ${{ matrix.python }} --with pytest pytest -q
-          code=$?
-          set -e
-          # exit 5 = no tests collected yet (skeleton); treat as success
-          if [ "$code" -ne 0 ] && [ "$code" -ne 5 ]; then exit "$code"; fi
+      # Coverage floor: a regression guard, not a target. Today's suite sits at
+      # ~79% line coverage; 75 leaves headroom for legitimately-untestable code
+      # (live transports, notebook plotting) while catching a big untested
+      # feature landing. Raise it as coverage grows; never lower it to merge.
+      - name: Tests (pytest + coverage floor)
+        run: >
+          uv run --python ${{ matrix.python }} --with pytest --with pytest-cov
+          pytest -q --cov=magnetics --cov-report=term --cov-fail-under=75
 
   typecheck:
     name: Typecheck (ty)
diff --git a/CLAUDE.md b/CLAUDE.md
index 85a8e8e..797985a 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -12,9 +12,9 @@ The full **fetch → process → service → GUI** path runs end-to-end for **bo
   `data/datafile/` (gitignored); read back via `magnetics.data.h5source`. The GUI can trigger a
   pull from the left rail (`PullControl` → `POST /api/fetch`).
 - **Process:** `core/spectral.py` (MODESPEC) is real and pure. The SLCONTOUR quasi-stationary fit
-  runs end-to-end via the reference pipeline in `magnetics._slcontour/` (xarray, self-contained
-  OMFIT shim) adapted by `core/qs_bridge` — real K / χ² / modes for shots pulled with the Bp LFS
-  midplane array. A pure `core/quasistationary` port exists but is not yet wired in production (#40).
+  runs end-to-end via the shim-free `core/qs_*` modules (`qs_io_data` → `qs_prep` → `qs_fit`,
+  adapted to nodes by `qs_bridge`; the former `_slcontour` translation, promoted into core) — real
+  K / χ² / modes for shots pulled with the Bp LFS midplane array.
 - **Service:** `service/app.py` — `GET /api/node/{shot}/{node_id}` serves `kind`-nodes from
   `service/nodes.py`; `/api/machines` lists fetched shots (mock fallback when none).
 - **Nodes / seam:** `nodes.py` forwards GUI query params and serves the core's real `mode_number` /
@@ -31,9 +31,10 @@ The full **fetch → process → service → GUI** path runs end-to-end for **bo
   the rotating/MODESPEC nodes and the Sensors view render those shots too.
 
 Known gaps / open work: per-sensor σ from the data layer (the QS fit currently uses a constant σ;
-helicity is computed from Ip·Bt); finishing the pure `core/quasistationary` port and wiring it in
-place of the `_slcontour` reference pipeline (#40); real equilibrium plotting in the Sensors tab
-(#43); Br saddle-loop geometry corrections (#44).
+helicity is computed from Ip·Bt); real equilibrium plotting in the Sensors tab (#43); the
+xarray-heavy `qs_*` modules are excluded from the `ty` typecheck (see `[tool.ty.src]` in
+pyproject.toml); the cursor analyses (`phase_fit` / mode shape) read from a 1–25 kHz cached
+array spectrum, so modes above 25 kHz aren't phase-fittable at the cursor.
 
 ## The API contract is FLEXIBLE — change it, don't fake around it
 The `kind`-node contract (`core/contracts.py` ⇄ `gui/web/src/lib/contract.ts`, plus the
@@ -59,16 +60,27 @@ rather than fabricating data in the GUI. Keep `contracts.py` and `contract.ts` i
 
 ## Layout
 The Python project **is the repo root** (a uv project, served as a webapp). `src/magnetics/`:
-- `core/` — device-agnostic math (geometry, basis, design metrics, `quasistationary`, `spectral`).
+- `core/` — device-agnostic math: `spectral` (MODESPEC), `mode_shape`, `geometry`, the
+  quasi-stationary fit (`qs_io_data`/`qs_prep`/`qs_fit`/`qs_bridge`/`qs_device`/`qs_run`/
+  `qs_plots`), and the `contracts` node shapes.
 - `data/` — sources + `fetch/` (toksearch/mdsthin pulls, cluster orchestration); device configs
   in `data/device/*.json`.
 - `service/` — FastAPI; the built GUI is bundled at `service/webapp/` and served here.
-- `_slcontour/` — the reference SLCONTOUR translation (self-contained OMFIT shim), pending port
-  into `core/quasistationary` (issue #40); **excluded from lint/typecheck** until then.
 
 Tests in `tests/`, maintainer scripts in `scripts/`. `gui/web/` — React + Vite + TypeScript
 frontend (its `dist/` is staged into `service/webapp/` for the wheel).
 
+## Running the tests
+- **Python:** `uv run pytest` from the repo root — the whole suite is offline and deterministic
+  (~5 s; synthetic fixtures are generated at collection time, no real data needed). Live-network
+  tests (GA gateway / PPPL flux) are env-gated and skip unless `MAGNETICS_GA_USER` /
+  `MAGNETICS_FLUX_USER` is set; they are manual-only, never wired into CI.
+- **Frontend:** `cd gui/web && npm run test` (vitest, one-shot; `npm run test:watch` to iterate).
+- **Everything CI runs:** `uv run ruff format --check .` + `uvx ruff check .` + `uv run pytest`
+  + `uv run ty check src/magnetics` (Python), and `npm run lint` + `npm run typecheck` +
+  `npm run test` + `npm run build` (in `gui/web/`). Run these locally before committing —
+  a plain `pytest` + `tsc` pass does NOT cover everything CI checks.
+
 ## Conventions
 - Physics lives in `src/magnetics/core` (pure, device-agnostic, testable); **no physics in the
   service routes**.
diff --git a/gui/web/src/components/PullControl.tsx b/gui/web/src/components/PullControl.tsx
index 9cb7e80..68e3b07 100644
--- a/gui/web/src/components/PullControl.tsx
+++ b/gui/web/src/components/PullControl.tsx
@@ -8,6 +8,7 @@
 // decimation to make a pull seconds instead of minutes.
 import { useEffect, useRef, useState } from "react";
 import { apiBase, startFetch, usingLiveBackend, type DeviceInfo } from "../lib/api";
+import { deviceDefaults } from "../lib/deviceDefaults";
 import { useStore } from "../store";
 
 export default function PullControl() {
@@ -57,21 +58,33 @@ export default function PullControl() {
   const esRef = useRef<EventSource | null>(null);
 
   // Snap backend / sensor-set / window / shot to sensible per-device defaults on a
-  // device change. A tree device (NSTX/KSTAR) has no cluster and no analysis→signal
-  // map, so it uses mdsthin over a NARROW window (its raw signals are ~15 MHz and
-  // seconds long — a wide window is gigabytes). An NSTX-style tree device requires a
-  // named sensor set; KSTAR's transport defaults to its declared arrays when none is
-  // chosen. Called from the device onChange handler — NOT a synchronous effect
-  // (react-hooks/set-state-in-effect).
+  // device change. The decision lives in lib/deviceDefaults.ts (pure + unit-tested,
+  // issue #64); this just applies it to the form state.
   function snapDeviceDefaults(d: DeviceInfo) {
-    const tree = d.access === "mdsplus_tree";
-    setBackend(d.remote_capable ? "remote" : "mdsthin");
-    setSensorSet(tree && !d.needs_ssh_creds ? (d.sensor_sets[0] ?? "") : "");
-    setTmin(tree ? "250" : "1000");
-    setTmax(tree ? "350" : "5000");
-    if (d.default_shot != null) setShot(String(d.default_shot));
+    const def = deviceDefaults(d);
+    setBackend(def.backend);
+    setSensorSet(def.sensorSet);
+    setTmin(def.tmin);
+    setTmax(def.tmax);
+    if (def.shot != null) setShot(def.shot);
   }
 
+  // The store owns `device`, and the left rail's Device picker (App.tsx) changes it
+  // WITHOUT going through this component's select — snapping only in the select's
+  // onChange left the backend/window/sensor-set at the previous device's values (a
+  // rail switch to a tree device then POSTed backend "remote" to a device with no
+  // cluster while the select displayed "mdsthin"). Snap once per device id, from
+  // wherever the change originated.
+  const snappedDevice = useRef<string | null>(null);
+  useEffect(() => {
+    if (!device || snappedDevice.current === device) return;
+    snappedDevice.current = device;
+    const d = devices.find((x) => x.id === device);
+    // eslint-disable-next-line react-hooks/set-state-in-effect -- per-device defaults must follow an external (rail) device change
+    if (d) snapDeviceDefaults(d);
+    // eslint-disable-next-line react-hooks/exhaustive-deps
+  }, [device, devices]);
+
   // close any open pull stream when the component unmounts
   useEffect(() => () => esRef.current?.close(), []);
   // Device is owned by the store; derive its capabilities. A tree device (NSTX/KSTAR)
@@ -150,11 +163,7 @@ export default function PullControl() {
       <h3>Pull a shot (live)</h3>
       {devices.length > 0 && (
         <select className="pull-input" value={device} aria-label="device"
-          onChange={(e) => {
-            const d = devices.find((x) => x.id === e.target.value);
-            setDevice(e.target.value);
-            if (d) snapDeviceDefaults(d);
-          }}>
+          onChange={(e) => setDevice(e.target.value)}>
           {devices.map((d) => (
             <option key={d.id} value={d.id}>{d.name}</option>
           ))}
diff --git a/gui/web/src/components/tabs/QuasiStationaryTab.tsx b/gui/web/src/components/tabs/QuasiStationaryTab.tsx
index 15b186f..8aa463f 100644
--- a/gui/web/src/components/tabs/QuasiStationaryTab.tsx
+++ b/gui/web/src/components/tabs/QuasiStationaryTab.tsx
@@ -58,26 +58,26 @@ function lineTraces(
   node: LineNode,
   opts?: { visible?: boolean[]; opacity?: number[]; palette?: string[] },
 ): Partial<Plotly.PlotData>[] {
-  const sigma = node.meta?.sigma as number[][] | undefined;
   const pal = opts?.palette ?? LINE_PALETTE;
   const traces: Partial<Plotly.PlotData>[] = [];
 
   node.series.forEach((s, i) => {
     const color = pal[i % pal.length];
-    const sig = sigma?.[i];
     const vis = opts?.visible?.[i] !== false;
     const opacity = opts?.opacity?.[i] ?? 1;
 
-    if (sig && vis) {
+    // ±1σ band from the contract's typed lower/upper fields (also what the
+    // HDF5 export writes, so the band on screen matches the downloaded data).
+    if (s.lower && s.upper && vis) {
       traces.push({
         type: "scatter", mode: "lines", x: s.x,
-        y: s.y.map((v, j) => v + sig[j]),
+        y: s.upper,
         line: { width: 0, color }, showlegend: false, hoverinfo: "skip",
         opacity,
       } as Partial<Plotly.PlotData>);
       traces.push({
         type: "scatter", mode: "lines", x: s.x,
-        y: s.y.map((v, j) => v - sig[j]),
+        y: s.lower,
         fill: "tonexty", fillcolor: hexToRgba(color, 0.45 * opacity),
         line: { width: 0, color }, showlegend: false, hoverinfo: "skip",
         opacity,
@@ -154,7 +154,10 @@ function CollapseHeader({
 export default function QuasiStationaryTab({ machine }: { machine: string }) {
   const dark = useDarkMode();
   const fontScale = useStore((s) => s.fontScale);
-  const { cursorMs, setCursorMs, machines } = useStore();
+  // Selective subscriptions: a whole-store destructure re-rendered this heavy tab
+  // on EVERY store change (each keystroke in the left rail's credential fields).
+  const setCursorMs = useStore((s) => s.setCursorMs);
+  const machines = useStore((s) => s.machines);
 
   // ── Analysis settings ─────────────────────────────────────────────
   const [ns, setNs]               = useState("1,2,3");
@@ -171,7 +174,7 @@ export default function QuasiStationaryTab({ machine }: { machine: string }) {
   const [uncertainty, setUncertainty]   = useState("2e-5");
   const [energyFraction, setEnergyFraction] = useState("0.98");
   const [fitBasis, setFitBasis]         = useState("sinusoidal-integral");
-  const [fitCond, setFitCond]           = useState("10.0");
+  const [fitCond, setFitCond]           = useState("1000"); // OMFIT SLCONTOUR inversion cutoff (1/rcond), not the K>10 trust threshold
   const [cutoffLo, setCutoffLo]         = useState("5.0");
   const [cutoffHi, setCutoffHi]         = useState("250.0");
 
@@ -208,6 +211,9 @@ export default function QuasiStationaryTab({ machine }: { machine: string }) {
 
   // ── Deferred fetch: only compute when user clicks Plot ────────────
   const [committedParams, setCommittedParams] = useState<Record<string, string> | null>(null);
+  // Bumped on every Plot click so an IDENTICAL param set still re-runs the fetch —
+  // without it a transient failure was unrecoverable except by jiggling a setting.
+  const [plotNonce, setPlotNonce] = useState(0);
 
   const qsParams = useMemo(() => {
     const p: Record<string, string> = {
@@ -236,10 +242,6 @@ export default function QuasiStationaryTab({ machine }: { machine: string }) {
     uncertainty, energyFraction, fitBasis, fitCond, cutoffLo, cutoffHi, excludedChannels,
   ]);
 
-  useEffect(() => {
-    if (cursorMs === 0) setCursorMs(3140);
-    // eslint-disable-next-line react-hooks/exhaustive-deps
-  }, []);
 
   // Linked time-axis zoom (declared here so the trim-window effect below can reset it).
   const [timeRange, setTimeRange] = useState<[number, number] | null>(null);
@@ -275,8 +277,11 @@ export default function QuasiStationaryTab({ machine }: { machine: string }) {
 
   // null fetchMachine suppresses all useNode fetches until initial commit fires.
   const fetchMachine = committedParams !== null ? machine : null;
-  // Highlight the Plot button when settings have drifted from the last commit.
-  const paramsDirty = committedParams !== null && committedParams !== qsParams;
+  // Highlight the Plot button when settings have drifted from the last commit —
+  // by VALUE (the memo object is recreated on any edit, so a reference compare
+  // flagged "settings changed" even after reverting to identical values).
+  const paramsDirty =
+    committedParams !== null && JSON.stringify(committedParams) !== JSON.stringify(qsParams);
 
   // Per-plot export helpers: a stable filename + the download descriptor (node id +
   // the params that produced it) so each figure exports its own image + HDF5 data.
@@ -286,32 +291,48 @@ export default function QuasiStationaryTab({ machine }: { machine: string }) {
   );
   const xn = useCallback((nodeId: string) => `shot_${machine}_${nodeId}`, [machine]);
 
-  // ── Node fetches — gated by fetchMachine (null until Plot is clicked) ──
-  const { node: qualityRaw } = useNode(fetchMachine, "fit_quality", committedParams ?? {});
+  // ── Node fetches — gated by fetchMachine (null until Plot is clicked).
+  // plotNonce re-runs identical-key fetches on Plot/retry. Errors are captured
+  // per node and surfaced — a failed node must never sit behind an eternal
+  // "loading…" placeholder.
+  const { node: qualityRaw } = useNode(fetchMachine, "fit_quality", committedParams ?? {}, plotNonce);
   const qualityNode = qualityRaw?.kind === "metrics" ? (qualityRaw as MetricsNode) : null;
 
   // Time-series nodes
-  const { node: phiTimeNode }                   = useNode(fetchMachine, "phi_t",      committedParams ?? {});
-  const { node: ampNode, error: ampError }       = useNode(fetchMachine, "amplitude",   committedParams ?? {});
-  const { node: phaseTimeNode }                 = useNode(fetchMachine, "phase_t",     committedParams ?? {});
+  const { node: phiTimeNode, error: phiTError }   = useNode(fetchMachine, "phi_t",     committedParams ?? {}, plotNonce);
+  const { node: ampNode, error: ampError }        = useNode(fetchMachine, "amplitude", committedParams ?? {}, plotNonce);
+  const { node: phaseTimeNode, error: phaseTError } = useNode(fetchMachine, "phase_t", committedParams ?? {}, plotNonce);
 
   // Sensor maps (R-Z cross-section + unrolled φ-θ), signal conditioning, fit quality
   // time series.
-  const { node: sensorRzRaw }    = useNode(fetchMachine, "sensor_map_rz",           committedParams ?? {});
-  const { node: sensorCylRaw }   = useNode(fetchMachine, "sensor_map_cylindrical", committedParams ?? {});
-  const { node: signalRaw }      = useNode(fetchMachine, "signal_conditioning",    committedParams ?? {});
-  const { node: chiSqRaw }       = useNode(fetchMachine, "chi_sq_t",               committedParams ?? {});
-  const { node: fitResRaw }      = useNode(fetchMachine, "fit_residuals",          committedParams ?? {});
-  const { node: svdEnergyRaw, error: svdEnergyError } = useNode(fetchMachine, "svd_energy",    committedParams ?? {});
-  const { node: svdCondRaw, error: svdCondError }     = useNode(fetchMachine, "svd_condition", committedParams ?? {});
+  const { node: sensorRzRaw, error: sensorRzError }   = useNode(fetchMachine, "sensor_map_rz",          committedParams ?? {}, plotNonce);
+  const { node: sensorCylRaw, error: sensorCylError } = useNode(fetchMachine, "sensor_map_cylindrical", committedParams ?? {}, plotNonce);
+  const { node: signalRaw, error: signalError }       = useNode(fetchMachine, "signal_conditioning",    committedParams ?? {}, plotNonce);
+  const { node: chiSqRaw, error: chiSqError }         = useNode(fetchMachine, "chi_sq_t",               committedParams ?? {}, plotNonce);
+  const { node: fitResRaw, error: fitResError }       = useNode(fetchMachine, "fit_residuals",          committedParams ?? {}, plotNonce);
+  const { node: svdEnergyRaw, error: svdEnergyError } = useNode(fetchMachine, "svd_energy",    committedParams ?? {}, plotNonce);
+  const { node: svdCondRaw, error: svdCondError }     = useNode(fetchMachine, "svd_condition", committedParams ?? {}, plotNonce);
 
   // No-data guard: 404 means the shot's HDF5 file hasn't been pulled yet.
   const noData = committedParams !== null && ampError?.includes("fetch failed (404)") === true;
-  // Fit-unavailable guard: a non-404 error means the quasi-stationary fit couldn't run
-  // (most often the shot was pulled for rotating-mode analysis and lacks the Bp
-  // LFS midplane array). Show the reason instead of a perpetual "loading…".
-  const fitUnavailable = committedParams !== null && !noData && ampError != null;
+  // Fit-unavailable guard: ONLY a 422 means the quasi-stationary fit can't run on
+  // this shot (no Bp LFS midplane array — the server's explicit answer). Any other
+  // error (500, timeout, network) is transient and must show as a retryable error,
+  // not a confidently wrong "this shot lacks the QS array" diagnosis.
+  const missingArray = committedParams !== null && ampError?.includes("fetch failed (422)") === true;
+  const fitFailed = committedParams !== null && !noData && !missingArray && ampError != null;
   const fitError = ampError?.replace(/^Error:\s*fetch failed \(\d+\):\s*/, "") ?? "";
+  // Nodes that errored while showing nothing — surfaced as one strip above the plots.
+  const nodeErrors = (
+    [
+      ["φ–t", phiTimeNode, phiTError],
+      ["phase", phaseTimeNode, phaseTError],
+      ["signals", signalRaw, signalError],
+      ["residuals", fitResRaw, fitResError],
+      ["χ²", chiSqRaw, chiSqError],
+      ["sensor map", sensorRzRaw, sensorRzError],
+    ] as const
+  ).filter(([, node, err]) => node === null && err != null);
 
   const sensorRzNode  = sensorRzRaw?.kind  === "line" ? (sensorRzRaw  as LineNode) : null;
   const sensorCylNode = sensorCylRaw?.kind === "line" ? (sensorCylRaw as LineNode) : null;
@@ -928,7 +949,7 @@ export default function QuasiStationaryTab({ machine }: { machine: string }) {
             <span style={{ fontSize: "calc(10px * var(--font-scale))", color: "var(--text-dim)" }}>settings changed</span>
           )}
           <button
-            onClick={() => setCommittedParams(qsParams)}
+            onClick={() => { setCommittedParams(qsParams); setPlotNonce(n => n + 1); }}
             style={{
               fontSize: "calc(11px * var(--font-scale))", padding: "2px 10px", borderRadius: 3, cursor: "pointer",
               background: paramsDirty ? "var(--accent)" : "var(--panel)",
@@ -1018,7 +1039,7 @@ export default function QuasiStationaryTab({ machine }: { machine: string }) {
           The HDF5 file for this shot has not been fetched yet.<br />
           Use the <strong>pull panel</strong> in the left sidebar to fetch the data, then click Plot.
         </div>
-      ) : fitUnavailable ? (
+      ) : missingArray ? (
         <div style={{ padding: 16, border: "1px solid var(--border)", borderRadius: 4,
                       color: "var(--text-dim)", fontSize: "calc(12px * var(--font-scale))", lineHeight: 1.6 }}>
           <strong>No quasi-stationary fit for shot {machine}.</strong><br />
@@ -1027,6 +1048,20 @@ export default function QuasiStationaryTab({ machine }: { machine: string }) {
           channels, or choose a QS-capable shot.<br />
           <span style={{ opacity: 0.7 }}>reason: {fitError}</span>
         </div>
+      ) : fitFailed ? (
+        <div style={{ padding: 16, border: "1px solid var(--border)", borderRadius: 4,
+                      color: "var(--text-dim)", fontSize: "calc(12px * var(--font-scale))", lineHeight: 1.6 }}>
+          <strong>The quasi-stationary fit failed for shot {machine}.</strong><br />
+          <span style={{ opacity: 0.7 }}>{fitError}</span><br />
+          <button
+            onClick={() => setPlotNonce(n => n + 1)}
+            style={{ marginTop: 6, fontSize: "calc(11px * var(--font-scale))", padding: "2px 10px",
+                     borderRadius: 3, cursor: "pointer", background: "var(--accent)", color: "#fff",
+                     border: "1px solid var(--border)" }}
+          >
+            Retry
+          </button>
+        </div>
       ) : (<>
 
       {/* ── Fit channels + sensor map (φ-θ) — collapsed by default, above the main plots ── */}
@@ -1058,20 +1093,41 @@ export default function QuasiStationaryTab({ machine }: { machine: string }) {
                     );
                   })}
                 </div>
-              ) : <div className="placeholder">loading channels…</div>}
+              ) : <div className="placeholder">{signalError ? `channels unavailable: ${signalError.replace(/^Error:\s*/, "")}` : "loading channels…"}</div>}
             </div>
             {/* right — φ-θ unrolled sensor map (R-Z lives in the Sensors tab) */}
             <div style={{ flex: 1, minWidth: 0 }}>
               <div style={{ fontSize: 10, color: "var(--text-dim)", marginBottom: 2 }}>unrolled φ-θ · {channelFilter}</div>
               {sensorCylNode
                 ? <Plot height={300} data={sensorCylData} layout={sensorCylLayout} exportName={xn("sensor_map_cylindrical")} download={dl("sensor_map_cylindrical")} />
-                : <div className="placeholder" style={{ height: 300 }}>loading…</div>
+                : <div className="placeholder" style={{ height: 300 }}>{sensorCylError ? `unavailable: ${sensorCylError.replace(/^Error:\s*/, "")}` : "loading…"}</div>
               }
             </div>
           </div>
         )}
       </div>
 
+      {/* Per-node failures: surfaced instead of leaving sections silently absent
+          or stuck on "loading…" — with a retry that re-runs the identical fetch. */}
+      {nodeErrors.length > 0 && (
+        <div style={{ padding: "6px 10px", border: "1px solid var(--border)", borderRadius: 4,
+                      fontSize: "calc(11px * var(--font-scale))", color: "var(--text-dim)",
+                      display: "flex", alignItems: "center", gap: 8, flexWrap: "wrap" }}>
+          <span>
+            {nodeErrors.map(([label, , err]) =>
+              `${label}: ${(err ?? "").replace(/^Error:\s*fetch failed \(\d+\):\s*/, "")}`).join(" · ")}
+          </span>
+          <button
+            onClick={() => setPlotNonce(n => n + 1)}
+            style={{ fontSize: "calc(10px * var(--font-scale))", padding: "1px 8px", borderRadius: 3,
+                     cursor: "pointer", background: "var(--panel)", color: "var(--text)",
+                     border: "1px solid var(--border)" }}
+          >
+            retry
+          </button>
+        </div>
+      )}
+
       {/* ── Section D+E: Time-series results — PRIMARY, at top ────────── */}
       <div style={{ display: "flex", flexDirection: "column", gap: 10 }}>
         {/* Section 8: Contour φ–t heatmap */}
@@ -1223,7 +1279,7 @@ export default function QuasiStationaryTab({ machine }: { machine: string }) {
           </span>
         </div>
         {!signalNode ? (
-          <div className="placeholder" style={{ height: 240 }}>loading signals…</div>
+          <div className="placeholder" style={{ height: 240 }}>{signalError ? `signals unavailable: ${signalError.replace(/^Error:\s*/, "")}` : "loading signals…"}</div>
         ) : signalStacked && signalPairs ? (
           <div>
             {signalPairs.map((pair, i) => {
@@ -1260,12 +1316,12 @@ export default function QuasiStationaryTab({ machine }: { machine: string }) {
               {/* Residuals — top */}
               {fitResNode
                 ? <Plot height={150} data={fitResData} layout={fitResLayout} onClick={seekTo} onRelayout={handleTimeRelayout} exportName={xn("fit_residuals")} download={dl("fit_residuals")} />
-                : <div className="placeholder" style={{ height: 150 }}>loading residuals…</div>
+                : <div className="placeholder" style={{ height: 150 }}>{fitResError ? `residuals unavailable: ${fitResError.replace(/^Error:\s*/, "")}` : "loading residuals…"}</div>
               }
               {/* Chi² — bottom */}
               {chiSqNode
                 ? <Plot height={130} data={chiSqData} layout={chiSqLayout} onClick={seekTo} onRelayout={handleTimeRelayout} exportName={xn("chi_sq_t")} download={dl("chi_sq_t")} />
-                : <div className="placeholder" style={{ height: 130 }}>loading χ²…</div>
+                : <div className="placeholder" style={{ height: 130 }}>{chiSqError ? `χ² unavailable: ${chiSqError.replace(/^Error:\s*/, "")}` : "loading χ²…"}</div>
               }
             </div>
             <div style={{ width: 190, flexShrink: 0, display: "flex", flexDirection: "column", gap: 8 }}>
@@ -1328,14 +1384,14 @@ export default function QuasiStationaryTab({ machine }: { machine: string }) {
               <div style={{ fontSize: "calc(10px * var(--font-scale))", color: "var(--text-dim)", marginBottom: 2 }}>cross-section (R-Z)</div>
               {sensorRzNode
                 ? <Plot height={220} data={sensorRzData} layout={sensorRzLayout} exportName={xn("sensor_map_rz")} download={dl("sensor_map_rz")} />
-                : <div className="placeholder">loading…</div>
+                : <div className="placeholder">{sensorRzError ? `unavailable: ${sensorRzError.replace(/^Error:\s*/, "")}` : "loading…"}</div>
               }
             </div>
             <div>
               <div style={{ fontSize: "calc(10px * var(--font-scale))", color: "var(--text-dim)", marginBottom: 2 }}>unrolled φ-θ</div>
               {sensorCylNode
                 ? <Plot height={220} data={sensorCylData} layout={sensorCylLayout} exportName={xn("sensor_map_cylindrical")} download={dl("sensor_map_cylindrical")} />
-                : <div className="placeholder">loading…</div>
+                : <div className="placeholder">{sensorCylError ? `unavailable: ${sensorCylError.replace(/^Error:\s*/, "")}` : "loading…"}</div>
               }
             </div>
           </div>
diff --git a/gui/web/src/components/tabs/RotatingTab.tsx b/gui/web/src/components/tabs/RotatingTab.tsx
index fbaa578..4a5b862 100644
--- a/gui/web/src/components/tabs/RotatingTab.tsx
+++ b/gui/web/src/components/tabs/RotatingTab.tsx
@@ -58,7 +58,11 @@ function PanelPlaceholder({ text, height = 200 }: { text: string; height?: numbe
 }
 
 export default function RotatingTab({ machine }: { machine: string }) {
-  const { cursorMs, setCursorMs } = useStore();
+  // Selective subscriptions: a whole-store destructure re-rendered this heavy tab
+  // (every Plotly panel rebuilt) on EVERY store change, including each keystroke
+  // in the left rail's credential fields.
+  const cursorMs = useStore((s) => s.cursorMs);
+  const setCursorMs = useStore((s) => s.setCursorMs);
   // Foreground ink that flips with the theme so the raw dB/dt trace stays visible on
   // the light plot background (it was hard-coded white → invisible in light mode).
   const dark = useStore((s) => s.theme === "dark");
@@ -204,7 +208,7 @@ export default function RotatingTab({ machine }: { machine: string }) {
   // Real toroidal mode-number map n(t,f) — a full-array fit per cell (resolves n=1,2,3,4…
   // that the 2-point estimate aliases away). Backs the "Mode n" toggle, gated server-side.
   // Honors the same resolution knob + band as the power view so the two stay consistent.
-  const { node: modeNumberNode } = useNode(machine, "mode_number", {
+  const { node: modeNumberNode, error: modeNumberError } = useNode(machine, "mode_number", {
     slice_duration: specSliceMs / 1000, fmin, fmax, n_amp_pct: powerGate, n_gate: nGate,
     ...smoothParams,
   });
@@ -476,7 +480,7 @@ export default function RotatingTab({ machine }: { machine: string }) {
         zrange: [-3, 0] as [number, number],
       };
     }
-  }, [hasStaticFiles, specNode, modeNumberNode, syntheticSpecNode, displayMode, fmin, fmax, powerGate, gateFrac]);
+  }, [hasStaticFiles, specNode, modeNumberNode, syntheticSpecNode, displayMode, fmin, fmax, gateFrac]);
 
   // Time scrubber source: keep the scrubber alive when the "n" mode map is still
   // computing (processedSpecNode null) by falling back to the loaded power spectrogram's
@@ -728,9 +732,12 @@ export default function RotatingTab({ machine }: { machine: string }) {
     if (!processedSpecNode) {
       // In "n" mode the heavy array-STFT mode_number node can still be loading after the
       // power spectrogram has already arrived (specLoading is false by then) — show an
-      // explicit indicator instead of a blank tile.
+      // explicit indicator instead of a blank tile. A server ERROR (e.g. 422: no usable
+      // toroidal array) must read as an error, not as fake eternal progress.
       if (displayMode === "n" && hasStaticFiles) {
-        return <div className="placeholder">Computing mode map…</div>;
+        return modeNumberError
+          ? <div className="placeholder">mode map unavailable: {modeNumberError.replace(/^Error:\s*fetch failed \(\d+\):\s*/, "")}</div>
+          : <div className="placeholder">Computing mode map…</div>;
       }
       return null;
     }
diff --git a/gui/web/src/components/tabs/SensorsTab.tsx b/gui/web/src/components/tabs/SensorsTab.tsx
index 90fe264..287645b 100644
--- a/gui/web/src/components/tabs/SensorsTab.tsx
+++ b/gui/web/src/components/tabs/SensorsTab.tsx
@@ -10,7 +10,7 @@
 // The node is a scatter2d (R-Z points) whose `meta` carries the full per-sensor
 // records + the vessel outline; the backend owns every device specific (which
 // family is Bp vs a saddle loop, the wall shape), so this view is device-agnostic.
-import { useMemo, useState } from "react";
+import { useEffect, useMemo, useState } from "react";
 import type * as Plotly from "plotly.js";
 import { useStore } from "../../store";
 import { useNode } from "../../lib/useNode";
@@ -122,6 +122,13 @@ export default function SensorsTab({ machine }: { machine: string }) {
   // User's explicit checkbox overrides; null until the first toggle, before which
   // the default selection (broadest Bp + Br set) is derived from `sets` at render.
   const [userSets, setUserSets] = useState<Record<string, boolean> | null>(null);
+  // Overrides are keyed by SET NAME, which is device-specific: carrying them to a
+  // different machine leaves every checkbox unmatched and renders an empty sensor
+  // scene with no hint why. Reset to the derived default on machine change.
+  useEffect(() => {
+    // eslint-disable-next-line react-hooks/set-state-in-effect -- drop per-device set overrides on machine change
+    setUserSets(null);
+  }, [machine]);
   const selectedSets = useMemo<Record<string, boolean>>(() => {
     if (userSets) return userSets;
     const init: Record<string, boolean> = {};
diff --git a/gui/web/src/lib/contract.ts b/gui/web/src/lib/contract.ts
index c41d747..35321df 100644
--- a/gui/web/src/lib/contract.ts
+++ b/gui/web/src/lib/contract.ts
@@ -6,7 +6,7 @@
 // adding a new analysis on the Python side needs ZERO frontend changes — the
 // physics teams just emit one of these shapes.
 //
-// This mirrors the Python builders (magfit/contracts.py): contour · heatmap ·
+// This mirrors the Python builders (magnetics/core/contracts.py): contour · heatmap ·
 // scatter2d · metrics, plus `line` for time-series traces. Keep the two in sync.
 // ─────────────────────────────────────────────────────────────────────────────
 
diff --git a/gui/web/src/lib/deviceDefaults.test.ts b/gui/web/src/lib/deviceDefaults.test.ts
new file mode 100644
index 0000000..17a1cbe
--- /dev/null
+++ b/gui/web/src/lib/deviceDefaults.test.ts
@@ -0,0 +1,74 @@
+// Per-device pull defaults (issue #64) — the decision the rail/pull device
+// pickers snap to. Shapes mirror the real /api/devices payloads.
+import { describe, expect, it } from "vitest";
+
+import type { DeviceInfo } from "./api";
+import { deviceDefaults } from "./deviceDefaults";
+
+const diiid: DeviceInfo = {
+  id: "diiid",
+  name: "DIII-D",
+  sensor_sets: ["Bp LFS midplane", "Bp All"],
+  access: "ptdata",
+  remote_capable: true,
+  default_shot: 184927,
+  needs_ssh_creds: false,
+};
+
+const nstx: DeviceInfo = {
+  id: "nstx",
+  name: "NSTX/NSTX-U",
+  sensor_sets: ["fastmag toroidal", "fastmag poloidal"],
+  access: "mdsplus_tree",
+  remote_capable: false,
+  default_shot: 204718,
+  needs_ssh_creds: false,
+};
+
+const kstar: DeviceInfo = {
+  id: "kstar",
+  name: "KSTAR",
+  sensor_sets: ["mirnov_toroidal"],
+  access: "mdsplus_tree",
+  remote_capable: false,
+  default_shot: 42477,
+  needs_ssh_creds: true, // own VPN+SSH transport
+};
+
+describe("deviceDefaults", () => {
+  it("DIII-D: fast cluster backend, wide flat-top window, no set override", () => {
+    expect(deviceDefaults(diiid)).toEqual({
+      backend: "remote",
+      sensorSet: "",
+      tmin: "1000",
+      tmax: "5000",
+      shot: "184927",
+    });
+  });
+
+  it("NSTX (tree device): mdsthin + narrow window + first sensor set required", () => {
+    expect(deviceDefaults(nstx)).toEqual({
+      backend: "mdsthin",
+      sensorSet: "fastmag toroidal",
+      tmin: "250",
+      tmax: "350",
+      shot: "204718",
+    });
+  });
+
+  it("KSTAR (own transport): tree window but NO preselected set (its arrays block is the default)", () => {
+    const d = deviceDefaults(kstar);
+    expect(d.backend).toBe("mdsthin"); // never "remote" — no cluster block
+    expect(d.sensorSet).toBe("");
+    expect(d.tmin).toBe("250");
+  });
+
+  it("keeps the user's shot when the device has no default", () => {
+    expect(deviceDefaults({ ...nstx, default_shot: null }).shot).toBeNull();
+    expect(deviceDefaults({ ...nstx, default_shot: undefined }).shot).toBeNull();
+  });
+
+  it("tree device with no sensor sets degrades to an empty set", () => {
+    expect(deviceDefaults({ ...nstx, sensor_sets: [] }).sensorSet).toBe("");
+  });
+});
diff --git a/gui/web/src/lib/deviceDefaults.ts b/gui/web/src/lib/deviceDefaults.ts
new file mode 100644
index 0000000..a4843f9
--- /dev/null
+++ b/gui/web/src/lib/deviceDefaults.ts
@@ -0,0 +1,29 @@
+// Per-device pull-form defaults — the decision PullControl snaps to whenever the
+// selected device changes (from its own select OR the left rail's picker).
+//
+// Extracted as a pure function so it is unit-testable (issue #64): a tree device
+// (NSTX/KSTAR) has no cluster and no analysis→signal map, so it pulls via
+// mdsthin over a NARROW window (its raw signals are ~15 MHz and seconds long — a
+// wide window is gigabytes). An NSTX-style tree device requires a named sensor
+// set; a device with its own SSH transport (KSTAR) defaults to its declared
+// arrays when none is chosen, so no set is preselected there.
+import type { DeviceInfo } from "./api";
+
+export interface DeviceDefaults {
+  backend: "remote" | "mdsthin";
+  sensorSet: string;
+  tmin: string;
+  tmax: string;
+  shot: string | null; // null = keep whatever shot the user typed
+}
+
+export function deviceDefaults(d: DeviceInfo): DeviceDefaults {
+  const tree = d.access === "mdsplus_tree";
+  return {
+    backend: d.remote_capable ? "remote" : "mdsthin",
+    sensorSet: tree && !d.needs_ssh_creds ? (d.sensor_sets[0] ?? "") : "",
+    tmin: tree ? "250" : "1000",
+    tmax: tree ? "350" : "5000",
+    shot: d.default_shot != null ? String(d.default_shot) : null,
+  };
+}
diff --git a/gui/web/src/lib/useNode.ts b/gui/web/src/lib/useNode.ts
index b144eeb..fd9b16e 100644
--- a/gui/web/src/lib/useNode.ts
+++ b/gui/web/src/lib/useNode.ts
@@ -1,38 +1,48 @@
 // Small hook: fetch a named analysis node for a machine and track load/error.
 // Views use this instead of hand-rolling fetch + useState each time.
 //
-// State is keyed by (machine, node, params). Stale-while-revalidate: when the key
-// changes we keep returning the last node while the new fetch is in flight, so a
-// consumer plot never unmounts/blanks mid-refetch (which crashes Plotly when e.g.
-// a time-slider changes the fetch params).
+// State is keyed by (machine, node, params). Stale-while-revalidate is
+// PARAM-scoped only: when params change we keep returning the last node while
+// the new fetch is in flight, so a consumer plot never unmounts/blanks
+// mid-refetch (which crashes Plotly when e.g. a time-slider changes the fetch
+// params). A MACHINE change drops the stale node immediately — the previous
+// shot's data must never render under the new shot's labels (or leak into its
+// PNG/HDF5 export descriptors).
+//
+// `retryKey`: bump to re-run the fetch for an IDENTICAL key — without it, a
+// transient failure was unrecoverable from the UI (the QS "Plot" button with
+// unchanged params re-ran no effect).
 import { useEffect, useState } from "react";
 import { fetchNode } from "./api";
 import type { Node } from "./contract";
 
-interface Entry { key: string; node: Node | null; error: string | null }
+interface Entry { key: string; machine: string | null; node: Node | null; error: string | null }
 
 export function useNode(
   machine: string | null,
   nodeId: string,
   params?: Record<string, string | number>,
+  retryKey: number = 0,
 ) {
   const key = `${machine}::${nodeId}::${params ? JSON.stringify(params) : ""}`;
-  const [entry, setEntry] = useState<Entry>({ key, node: null, error: null });
+  const [entry, setEntry] = useState<Entry>({ key, machine, node: null, error: null });
 
   useEffect(() => {
     if (!machine) return;
     let alive = true;
     fetchNode(machine, nodeId, params)
-      .then((n) => { if (alive) setEntry({ key, node: n, error: null }); })
-      .catch((e) => { if (alive) setEntry({ key, node: null, error: String(e) }); });
+      .then((n) => { if (alive) setEntry({ key, machine, node: n, error: null }); })
+      .catch((e) => { if (alive) setEntry({ key, machine, node: null, error: String(e) }); });
     return () => { alive = false; };
     // eslint-disable-next-line react-hooks/exhaustive-deps
-  }, [machine, nodeId, key]);
+  }, [machine, nodeId, key, retryKey]);
 
   const fresh = entry.key === key;
+  const node = fresh || entry.machine === machine ? entry.node : null;
   return {
-    node: entry.node,                                  // keep stale node during refetch
+    node,                                              // stale only across param changes
     error: fresh ? entry.error : null,
-    loading: entry.node === null && !(fresh && entry.error),
+    loading: node === null && !(fresh && entry.error),
+    refetching: !fresh && node !== null,               // stale shown, new fetch in flight
   };
 }
diff --git a/pyproject.toml b/pyproject.toml
index 83aea35..3e10d35 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -22,6 +22,9 @@ dependencies = [
     "uvicorn[standard]>=0.49.0",
     # per-user data dir for installed (non-checkout) runs; see data/h5source.py
     "platformdirs>=4",
+    # KSTAR VPN/SSH tunnel automation (data/fetch/kstar_transport.py imports it
+    # lazily; without declaring it, a wheel install dies at live-pull time)
+    "pexpect>=4.9",
 ]
 
 [build-system]
diff --git a/run-dev.sh b/run-dev.sh
index ac1d997..5e1e639 100755
--- a/run-dev.sh
+++ b/run-dev.sh
@@ -3,8 +3,8 @@
 # One command to run the GUI for DEVELOPMENT.
 #
 # DEFAULT is LIVE: the GUI talks to the real FastAPI service, and NO mock data is
-# served or renderable in this mode. The rotating-mode (MODESPEC) path serves real
-# analysis from fetched shots; the quasi-stationary fit stream is still a stub. Use
+# served or renderable in this mode. Both the rotating-mode (MODESPEC) and the
+# quasi-stationary (SLCONTOUR) analyses serve real results from fetched shots. Use
 # `static` only for offline frontend work against the demo fixtures.
 #
 #   ./run-dev.sh           live    — FastAPI service (:8000) + GUI (:5173) against it (default)
diff --git a/src/magnetics/core/contracts.py b/src/magnetics/core/contracts.py
index d5c80b0..3bf1f8c 100644
--- a/src/magnetics/core/contracts.py
+++ b/src/magnetics/core/contracts.py
@@ -8,8 +8,11 @@
 
 from __future__ import annotations
 
+import math
 from typing import Any
 
+import numpy as np
+
 
 def _clean(d: dict[str, Any]) -> dict[str, Any]:
     """Drop None-valued optional keys so the JSON matches the TS optionals."""
@@ -98,11 +101,24 @@ def metrics(title, fields, *, meta=None) -> dict:
 def quality_for_k(k: float) -> str:
     """SLCONTOUR condition-number thresholds (warn > 10, error > 20).
 
-    Mirrors `qualityForK` in contract.ts so the GUI's traffic-light coloring and
-    the backend agree.
+    Mirrors `qualityForK` in contract.ts (which tests `!Number.isFinite(K)`) so
+    the GUI's traffic-light coloring and the backend agree — including ±inf.
     """
-    if not (k == k) or k > 20:  # NaN or too ill-conditioned
+    if not math.isfinite(k) or k > 20:  # NaN/±inf or too ill-conditioned
         return "bad"
     if k > 10:
         return "warn"
     return "good"
+
+
+def json_finite(a):
+    """Array → nested lists with every non-finite value as None (JSON null).
+
+    The HTTP layer serializes with ``allow_nan=False``, so a single NaN sample in
+    a raw fetched channel otherwise turns a whole node response into an opaque
+    500. The GUI's Plotly renderers treat null as a gap, which is the honest
+    picture of missing samples."""
+    a = np.asarray(a, dtype=float)
+    out = a.astype(object)
+    out[~np.isfinite(a)] = None
+    return out.tolist()
diff --git a/src/magnetics/core/mode_shape.py b/src/magnetics/core/mode_shape.py
index 2a8d70f..8ddcadb 100644
--- a/src/magnetics/core/mode_shape.py
+++ b/src/magnetics/core/mode_shape.py
@@ -297,16 +297,20 @@ def mac_n_spectrum(
     *,
     n_range: tuple[int, int] = (-6, 6),
 ) -> tuple[NDArray[np.integer], NDArray[np.floating], int]:
-    """MAC of a measured shape vector against ideal pure-mode templates e^{−inφ}.
+    """MAC of a measured shape vector against ideal pure-mode templates e^{+inφ}.
 
     Gives a *shape-based*, geometry-aware toroidal mode-number identification (how
     closely the measured array pattern resembles each pure rotating mode), a useful
-    cross-check on the cross-phase fit. Returns (n_values, mac_values, best_n).
+    cross-check on the cross-phase fit. The measured shape comes from the
+    conj(sig)·ref cross-phases (``shape_vector`` of ``extract_mode_at_frequency`` /
+    ``mode_from_spectrum`` output), where a ``cos(nφ - ωt)`` mode appears as a
+    ``+nφ`` phase ramp, i.e. ``z ∝ e^{+inφ}`` — so matching that template reports
+    the same signed n as the cross-phase fit. Returns (n_values, mac_values, best_n).
     """
     phi = np.deg2rad(np.asarray(angle_deg, dtype=np.float64))
     z = np.asarray(complex_shape, dtype=np.complex128)
     ns = np.arange(n_range[0], n_range[1] + 1)
-    macs = np.array([mac(z, np.exp(-1j * n * phi)) for n in ns])
+    macs = np.array([mac(z, np.exp(1j * n * phi)) for n in ns])
     return ns, macs, int(ns[int(np.argmax(macs))])
 
 
diff --git a/src/magnetics/core/qs_bridge.py b/src/magnetics/core/qs_bridge.py
index 3c42c88..277d4ff 100644
--- a/src/magnetics/core/qs_bridge.py
+++ b/src/magnetics/core/qs_bridge.py
@@ -117,7 +117,8 @@ def fit_to_qs_fit_node(
         theta_grid = np.linspace(0, 360, 49)
 
     z = _reconstruct_grid(fit_ds, phi_grid, theta_grid, t_idx) * _T_TO_G
-    zmax = float(np.nanmax(np.abs(z))) or 1.0
+    finite = np.abs(z[np.isfinite(z)])
+    zmax = float(finite.max()) if finite.size and finite.max() > 0 else 1.0
 
     overlay = None
     if sensor_phis is not None and sensor_thetas is not None:
@@ -147,7 +148,7 @@ def fit_to_qs_fit_node(
             "n": int(round(ns[dominant_idx])),
             "m": int(round(ms[dominant_idx])),
             "condition_number": round(K, 2),
-            "red_chi_sq": round(chi2_t, 3),
+            "red_chi_sq": round(chi2_t, 3) if np.isfinite(chi2_t) else None,
             "shot": str(fit_ds.attrs.get("shot", "")),
         },
     )
@@ -156,8 +157,11 @@ def fit_to_qs_fit_node(
 def fit_to_amplitude_node(fit_ds) -> dict:
     """Mode amplitude ± 1σ vs time for each fitted mode → LineNode.
 
-    meta.sigma[i] is the per-series error band consumed by lineTraces() in the GUI.
-    meta.legend_title is "n" when all poloidal m=0, else "m/n" (mirrors plot_fit_modes).
+    The ±1σ band rides in each series' typed ``lower``/``upper`` fields (the
+    contract's error-band channel, same as the rotating mode-shape nodes) so the
+    GUI band and the HDF5 export are the same data — not an untyped meta side
+    channel the exporter drops. meta.legend_title is "n" when all poloidal m=0,
+    else "m/n" (mirrors plot_fit_modes).
     """
     from ..core import contracts
 
@@ -168,17 +172,19 @@ def fit_to_amplitude_node(fit_ds) -> dict:
     sigmas = fit_ds["fit_sigmas"].values  # [mode, time] complex (constant along time)
 
     series = []
-    sigma_bands = []
     for i, (n, m) in enumerate(zip(ns, ms_vals)):
         amp, amp_err, _, _ = _amp_phase(coeffs[i], sigmas[i])
+        y = amp * _T_TO_G
+        err = amp_err * _T_TO_G
         series.append(
             {
                 "name": _mode_label(n, m),
                 "x": np.round(t_ms, 2).tolist(),
-                "y": np.round(amp * _T_TO_G, 4).tolist(),
+                "y": np.round(y, 4).tolist(),
+                "lower": np.round(y - err, 4).tolist(),
+                "upper": np.round(y + err, 4).tolist(),
             }
         )
-        sigma_bands.append(np.round(amp_err * _T_TO_G, 4).tolist())
 
     legend_title = "n" if np.all(ms_vals == 0) else "m/n"
     K = float(fit_ds.attrs.get("condition_number", fit_ds.attrs.get("raw_cn", 0.0)))
@@ -186,7 +192,6 @@ def fit_to_amplitude_node(fit_ds) -> dict:
         series,
         {"x": "time (ms)", "y": "amplitude (G)"},
         meta={
-            "sigma": sigma_bands,
             "legend_title": legend_title,
             "condition_number": round(K, 2),
             "shot": str(fit_ds.attrs.get("shot", "")),
@@ -215,7 +220,6 @@ def fit_to_phase_t_node(fit_ds) -> dict:
     phase_visible = [bool(a > 0.1 * max_amp) for a in p90_amps]
 
     series = []
-    sigma_bands = []
     for i, (n, m) in enumerate(zip(ns, ms_vals)):
         _, _, phase, phase_err = _amp_phase(coeffs[i], sigmas[i])
         series.append(
@@ -223,16 +227,16 @@ def fit_to_phase_t_node(fit_ds) -> dict:
                 "name": _mode_label(n, m),
                 "x": np.round(t_ms, 2).tolist(),
                 "y": np.round(phase, 3).tolist(),
+                "lower": np.round(phase - phase_err, 3).tolist(),
+                "upper": np.round(phase + phase_err, 3).tolist(),
             }
         )
-        sigma_bands.append(np.round(phase_err, 3).tolist())
 
     K = float(fit_ds.attrs.get("condition_number", fit_ds.attrs.get("raw_cn", 0.0)))
     return contracts.line(
         series,
         {"x": "time (ms)", "y": "phase (deg)"},
         meta={
-            "sigma": sigma_bands,
             "phase_visible": phase_visible,
             "condition_number": round(K, 2),
             "shot": str(fit_ds.attrs.get("shot", "")),
@@ -248,7 +252,7 @@ def fit_to_fit_quality_node(fit_ds) -> dict:
     eff_cn = float(fit_ds.attrs.get("eff_cn", K))
     mean_chi2 = float(np.nanmean(fit_ds["red_chi_sq"].values))
     n_ch = int(fit_ds.sizes.get("channel", fit_ds["fit_sigmas"].shape[0]))
-    fit_cond = float(fit_ds.attrs.get("fit_condition", 10.0))
+    fit_cond = float(fit_ds.attrs.get("fit_condition", 1e3))
     n_modes = int(fit_ds.sizes["mode"])
 
     return contracts.metrics(
@@ -297,7 +301,9 @@ def fit_to_phi_t_node(fit_ds, theta_fixed_deg: float = 0.0, n_phi: int = 73) ->
         z += (coeffs[i][None, :] * phase_phi * phase_theta).real
     z *= _T_TO_G
 
-    zmax = float(np.nanpercentile(np.abs(z), 99)) or 1.0
+    finite = np.abs(z[np.isfinite(z)])
+    zmax = float(np.percentile(finite, 99)) if finite.size else 1.0
+    zmax = zmax or 1.0
     K = float(fit_ds.attrs.get("condition_number", fit_ds.attrs.get("raw_cn", 0.0)))
 
     # z shape: [n_phi, n_time] but ContourNode z is [n_y][n_x]
@@ -378,7 +384,7 @@ def fit_to_svd_condition_node(fit_ds) -> dict:
     w_a = np.linalg.svd(A, compute_uv=False)
     cond = np.abs(w_a[0] / w_a)
     idx = np.arange(1, len(cond) + 1)
-    fit_cond = float(fit_ds.attrs.get("fit_condition", 10.0))
+    fit_cond = float(fit_ds.attrs.get("fit_condition", 1e3))
 
     series = {
         "name": "condition number",
diff --git a/src/magnetics/core/qs_fit.py b/src/magnetics/core/qs_fit.py
index 14113cf..cecc589 100644
--- a/src/magnetics/core/qs_fit.py
+++ b/src/magnetics/core/qs_fit.py
@@ -171,7 +171,7 @@ def fit(
     fit_exclude=(),
     fit_basis="sinusoidal-integral",
     fit_geometry="cylindrical",
-    fit_cond=10.0,
+    fit_cond=1e3,
     sigma_override=None,
     ncenters=6,
     mcenters=1,
@@ -191,6 +191,11 @@ def fit(
         ``'gaussian-point'``, or ``'gaussian-integral'``.
     :param fit_geometry: ``'cylindrical'`` (phi, theta) or ``'vertical'`` (phi, z).
     :param fit_cond: condition-number cutoff for the lstsq inversion (= 1/rcond).
+        Basis directions with condition > fit_cond are zeroed. The default matches
+        OMFIT SLCONTOUR's 1e3 — this is the *inversion* cutoff, NOT the "warn when
+        K > 10" trust threshold (``contracts.quality_for_k``); setting it that low
+        silently truncates directions the reference fit would keep, while making
+        the reported K(eff) <= 10 by construction.
     :param sigma_override: when given, use this uniform measurement uncertainty
         for every channel instead of the per-channel ``signal_sigma`` baked into
         the dataset at load time.
@@ -449,9 +454,15 @@ def _printv(*a):
 
     # ── assemble output Dataset ───────────────────────────────────────────────
     fit_b = np.dot(A, fit_coeffs).real.reshape(ds["channel"].shape[0], -1)
-    ds["fit_signal"] = xr.DataArray(fit_b, coords=ds["signal"].coords, dims=ds["signal"].dims) * ds[
-        "signal_sigma"
-    ].fillna(float(sigma.mean()))
+    # De-normalize by the SAME effective sigma that normalized the RHS `b` above
+    # (it honors sigma_override and NaN fills). Using the dataset's raw
+    # signal_sigma here would scale fit_signal by signal_sigma/override whenever
+    # they differ, corrupting residual and chi_sq while the coefficients stay
+    # correct.
+    sig_eff = xr.DataArray(sigma, coords={"channel": ds["channel"]}, dims=("channel",))
+    ds["fit_signal"] = (
+        xr.DataArray(fit_b, coords=ds["signal"].coords, dims=ds["signal"].dims) * sig_eff
+    )
     ds["fit_ns"] = xr.DataArray(nms_arr[:, 0], coords={"mode": np.arange(len(nms))}, dims=("mode",))
     ds["fit_ms"] = xr.DataArray(nms_arr[:, 1], coords={"mode": np.arange(len(nms))}, dims=("mode",))
     ds["fit_coeffs"] = xr.DataArray(
@@ -462,8 +473,7 @@ def _printv(*a):
     )
 
     ds["residual"] = ds["signal"] - ds["fit_signal"]
-    sig_for_chi = xr.DataArray(sigma, coords={"channel": ds["channel"]}, dims=("channel",))
-    ds["chi_sq"] = ((ds["residual"] / sig_for_chi) ** 2).sum("channel")
+    ds["chi_sq"] = ((ds["residual"] / sig_eff) ** 2).sum("channel")
     nu = max(b.shape[0] - rank_fit, 1)
     ds["red_chi_sq"] = ds["chi_sq"] / nu
 
diff --git a/src/magnetics/core/spectral.py b/src/magnetics/core/spectral.py
index 04a28f6..a41b5cf 100644
--- a/src/magnetics/core/spectral.py
+++ b/src/magnetics/core/spectral.py
@@ -165,6 +165,17 @@ def _cross_spectral_errors(
     return phase_err, amp_err
 
 
+def wrap_angle_deg(angle: float) -> float:
+    """Reduce an angle (deg) to the signed principal value in (-180, 180].
+
+    A mode's cross-phase n·Δφ is only defined modulo 360°, so the 2-point estimator
+    n = round(phase / Δφ) needs the *short-way, signed* probe separation: probes at
+    20° and 340° are 320° apart measured one way, but the wrapped phase sees -40°.
+    Feeding the long-way separation aliases every mode toward n = 0."""
+    wrapped = (float(angle) + 180.0) % 360.0 - 180.0
+    return 180.0 if wrapped == -180.0 else wrapped
+
+
 def cross_spectrum(
     sig1: NDArray[np.floating],
     sig2: NDArray[np.floating],
@@ -180,7 +191,10 @@ def cross_spectrum(
         sig1 (ndarray): first probe time series.
         sig2 (ndarray): second probe time series.
         sample_rate (float): sampling rate (Hz).
-        delta_phi (float | None): toroidal separation (deg); enables mode-number output.
+        delta_phi (float | None): toroidal separation (deg), phi2 - phi1; enables
+            mode-number output. Wrapped internally to the principal value (-180, 180],
+            so either the short- or long-way separation may be passed. Modes are only
+            resolved unaliased for |n| <= 180 / |wrapped delta_phi|.
         nperseg (int | None): Welch segment length for csd/coherence; None uses scipy's
             default (256). Set below the signal length on short slices to keep multiple
             averaging segments — coherence is meaningless from a single segment.
@@ -209,8 +223,9 @@ def cross_spectrum(
             n_segments=k,
         )
 
+    delta_phi = wrap_angle_deg(delta_phi)
     if delta_phi == 0:
-        raise ValueError("delta_phi must be non-zero to compute mode numbers")
+        raise ValueError("delta_phi must be non-zero (mod 360°) to compute mode numbers")
 
     mode = np.rint(phase / delta_phi).astype(np.intp)
 
@@ -323,15 +338,17 @@ def compute_spectrogram(
 
     Engine: one batched, single-precision short-time FFT per probe (no per-window loop),
     with the column count decimated to ``max_columns`` so cost scales with the display,
-    not the record length. Physics: cross-power = |conj(S1)·S2|, frequency-smoothed
-    coherence, n = round(phase / delta_phi), and per-mode RMS — identical to the 2-point
-    definitions in ``cross_spectrum``.
+    not the record length. Physics: cross-power = |conj(S1)·S2| scaled to a one-sided
+    cross-spectral density (csd's ``scaling="density"``), frequency-smoothed coherence,
+    n = round(phase / delta_phi), and per-mode RMS — identical to the 2-point
+    definitions in ``cross_spectrum``, so amplitudes are comparable across the two.
 
     Inputs:
         time (ndarray): sample times (s); assumed uniformly sampled.
         sig1 (ndarray): first probe time series.
         sig2 (ndarray): second probe time series.
-        delta_phi (float): toroidal separation (deg); must be non-zero.
+        delta_phi (float): toroidal separation (deg), phi2 - phi1; must be non-zero
+            mod 360°. Wrapped internally to (-180, 180] like ``cross_spectrum``.
         slice_duration (float): FFT window width (s) — sets the frequency resolution.
         window (str): scipy window name for the taper.
         max_columns (int): cap on spectrogram time columns (decimation lever).
@@ -340,8 +357,9 @@ def compute_spectrogram(
         result (SpectrogramResult): power/coherence/mode_number as (n_times, n_freqs)
             arrays, with time/frequency/rms_by_mode/mode_indices.
     """
+    delta_phi = wrap_angle_deg(delta_phi)
     if delta_phi == 0:
-        raise ValueError("delta_phi must be non-zero to compute mode numbers")
+        raise ValueError("delta_phi must be non-zero (mod 360°) to compute mode numbers")
 
     time = np.asarray(time)
     s1 = np.ascontiguousarray(sig1, dtype=np.float32)
@@ -366,7 +384,17 @@ def compute_spectrogram(
     # order matches cross_spectrum's scipy.signal.csd(sig2, sig1) convention so the
     # spectrogram and the single-window 2-point analysis report the same *signed* n.
     cross = np.conj(spec2) * spec1
-    power = np.abs(cross)
+    # Scale |cross| to a one-sided cross-spectral DENSITY exactly like scipy's
+    # csd(scaling="density"): 1/(fs·Σwin²), doubled off the DC/Nyquist bins. Without
+    # this, rms_by_mode's sqrt(Σ power·df) is off by a data-dependent factor of
+    # thousands vs the 2-point cross_spectrum RMS it claims to match. (Coherence and
+    # phase are ratios/angles — unaffected.)
+    scale = 1.0 / (sample_rate * float(np.sum(win.astype(np.float64) ** 2)))
+    power = np.abs(cross) * scale
+    if n_fft % 2 == 0:
+        power[:, 1:-1] *= 2.0
+    else:
+        power[:, 1:] *= 2.0
     phase = np.rad2deg(np.angle(cross))
 
     # Coherence needs averaging; smooth the auto/cross spectra over frequency bins.
@@ -605,7 +633,7 @@ class ArrayModeSpectrogram:
     time: NDArray[np.floating]  # (n_times,) window-center times (s)
     freq_band: NDArray[np.floating]  # (n_band,) Hz kept
     mode_number: NDArray[np.integer]  # (n_times, n_band) best-fit toroidal n
-    amplitude: NDArray[np.floating]  # (n_times, n_band) |Σ_p Z_p e^{-inφ}|
+    amplitude: NDArray[np.floating]  # (n_times, n_band) |⟨Z, e^{-inφ}⟩| = |Σ_p Z_p e^{+inφ}|
     quality: NDArray[np.floating]  # (n_times, n_band) harmonic energy fraction ∈ [1/M,1]
 
 
@@ -625,11 +653,13 @@ def array_mode_spectrogram(
 
     The per-cell mode-coherence is the *energy fraction* the single best-fit n captures of
     the summed toroidal-harmonic power, ``|R_{n*}|² / Σ_n |R_n|² ∈ [1/M, 1]`` with
-    ``R_n = Σ_p Z_p e^{-inφ_p}`` and ``M = 2·n_max + 1`` candidates: 1 = a pure single-n
+    ``R_n = Σ_p Z_p e^{+inφ_p}`` and ``M = 2·n_max + 1`` candidates: 1 = a pure single-n
     pattern, 1/M = white across harmonics (incoherent noise). This is a spectral-
     concentration ratio (more noise/signal contrast than the resultant length ``|R_{n*}| /
-    Σ_p|Z_p|``, whose noise floor sits higher at ~1/√P). The ``exp(-i n φ)`` sign matches
-    ``mode_from_spectrum`` and the phase fit, so the reported n agrees with them.
+    Σ_p|Z_p|``, whose noise floor sits higher at ~1/√P). ``R_n`` is the inner product
+    ``⟨Z, e^{-inφ}⟩`` (template conjugated): a ``cos(nφ - ωt)`` mode's positive-frequency
+    STFT pattern is ``Z_p ∝ e^{-inφ_p}``, so the projection peaks at the same signed +n
+    that ``cross_spectrum`` and the phase fit report.
 
     Inputs:
         spectrum (ArrayShapeSpectrum): the per-probe complex STFT band.
@@ -641,8 +671,8 @@ def array_mode_spectrogram(
     z = np.asarray(spectrum.spec)  # (P, T, F) complex
     phi = np.deg2rad(np.asarray(angle_deg, dtype=np.float64))
     ns = np.arange(-int(n_max), int(n_max) + 1)
-    basis = np.exp(-1j * ns[:, None] * phi[None, :])  # (M, P)
-    proj = np.einsum("mp,ptf->mtf", basis, z)  # (M, T, F) = R_n
+    basis = np.exp(1j * ns[:, None] * phi[None, :])  # (M, P) = conj(e^{-inφ_p}) templates
+    proj = np.einsum("mp,ptf->mtf", basis, z)  # (M, T, F) = R_n = ⟨Z, e^{-inφ}⟩
     amp = np.abs(proj)
     k = np.argmax(amp, axis=0)  # (T, F)
     peak = np.take_along_axis(amp, k[None], axis=0)[0]  # |R_{n*}| (T, F)
@@ -867,11 +897,14 @@ def fit_toroidal_mode(
     """Fit a toroidal mode number to per-probe phase-vs-angle data.
 
     A rotating mode of toroidal number n imprints a linear phase ramp across the
-    toroidal array: ``phase(phi) = c - n * phi`` (deg). Rather than unwrap the
-    (mod-360) phases — ill-posed for sparse arrays — this scans integer candidates
-    and picks the n whose residual ``phase + n * phi`` clusters most tightly on the
-    circle (largest amplitude-weighted resultant). The intercept c is the angle of
-    that resultant, giving a wrap-free fit line for the GUI to draw.
+    toroidal array. With the conj(sig)·ref cross-phase convention used by
+    ``extract_mode_at_frequency`` / ``mode_from_spectrum``, a ``cos(nφ - ωt)`` mode
+    measures as ``phase(phi) = c + n * phi`` (deg) — the same signed n the 2-point
+    ``cross_spectrum`` reports. Rather than unwrap the (mod-360) phases — ill-posed
+    for sparse arrays — this scans integer candidates and picks the n whose residual
+    ``phase - n * phi`` clusters most tightly on the circle (largest
+    amplitude-weighted resultant). The intercept c is the angle of that resultant,
+    giving a wrap-free fit line for the GUI to draw.
 
     Inputs:
         mode_result (ModeAtFrequencyResult): per-probe phase/amplitude/angle, e.g.
@@ -902,7 +935,7 @@ def fit_toroidal_mode(
     inter = np.empty(len(candidates))  # per-candidate intercept c_n (deg)
     rmag = np.empty(len(candidates))  # per-candidate resultant length
     for j, n in enumerate(candidates):
-        resultant = np.sum(w * np.exp(1j * np.deg2rad(phase + n * phi))) / w.sum()
+        resultant = np.sum(w * np.exp(1j * np.deg2rad(phase - n * phi))) / w.sum()
         rmag[j] = np.abs(resultant)
         inter[j] = np.rad2deg(np.angle(resultant))
     best_j = int(np.argmax(rmag))
@@ -927,7 +960,7 @@ def _wrap180(x):
 
     chi2 = np.array(
         [
-            np.sum((_wrap180(phase + n * phi - inter[j]) / sigma) ** 2)
+            np.sum((_wrap180(phase - n * phi - inter[j]) / sigma) ** 2)
             for j, n in enumerate(candidates)
         ]
     )
diff --git a/src/magnetics/data/fetch/toksearch.py b/src/magnetics/data/fetch/toksearch.py
index 54d6ab4..584291f 100644
--- a/src/magnetics/data/fetch/toksearch.py
+++ b/src/magnetics/data/fetch/toksearch.py
@@ -99,6 +99,26 @@ def _resolve_pointnames(dev, pointnames, shot):
         if pt is None:  # out of range / NotAvailable
             skipped.append(cid)
             continue
+        prev = canonical_of.get(pt)
+        if prev is not None and prev != cid:
+            # A per-era alias collides with another sensor's queried pointname
+            # (e.g. KSTAR \MC1T10 -> \MC1P03 since shot 17377, while \MC1P03 is
+            # itself in the pull). Fetching one node under two canonical ids
+            # previously corrupted the relabel step — both channels got the same
+            # name and the writer silently dropped one. Query the node once,
+            # prefer the sensor literally named `pt` (identity mapping), and
+            # report the losing id as skipped instead of losing it silently.
+            if pt == cid:
+                skipped.append(prev)
+                canonical_of[pt] = cid
+            else:
+                skipped.append(cid)
+            sys.stderr.write(
+                f"pointname collision at shot {shot_i}: {skipped[-1]!r} aliases "
+                f"{pt!r}, which is fetched as {canonical_of[pt]!r}; skipping the "
+                "alias (check the device file's per-era pointname overrides).\n"
+            )
+            continue
         query.append(pt)
         canonical_of[pt] = cid
     return query, canonical_of, skipped
@@ -213,16 +233,44 @@ def _default_progress(frac: float, msg: str) -> None:
 
 # --- channel record -----------------------------------------------------------
 class Channel:
-    """One fetched pointname: time (ms, float64) + data (float32)."""
+    """One fetched pointname: time (ms, float64) + data (float32).
+
+    ``error_kind`` classifies a failure (``ok=False``): ``"no_data"`` — the server
+    answered and the channel genuinely has nothing (cached as missing so the next
+    incremental pull skips it); ``"transport"`` — the connection/tunnel died (NOT
+    cached, so the next pull retries without --force; see issue #63)."""
 
-    __slots__ = ("name", "time", "data", "ok", "error")
+    __slots__ = ("name", "time", "data", "ok", "error", "error_kind")
 
-    def __init__(self, name, time=None, data=None, ok=False, error=""):
+    def __init__(self, name, time=None, data=None, ok=False, error="", error_kind="no_data"):
         self.name = name
         self.time = time
         self.data = data
         self.ok = ok
         self.error = error
+        self.error_kind = error_kind
+
+
+def _error_kind(exc) -> str:
+    """Classify a fetch exception: ``"transport"`` for connection-level failures
+    (socket/tunnel/timeout — worth retrying next run), ``"no_data"`` for a genuine
+    MDSplus no-data / no-node answer (cacheable as missing). Message matching is
+    the fallback for errors wrapped in strings by intermediate layers."""
+    if isinstance(exc, (OSError, EOFError, TimeoutError, ConnectionError)):
+        return "transport"
+    msg = str(exc).lower()
+    transport_words = (
+        "connection",
+        "timed out",
+        "timeout",
+        "broken pipe",
+        "reset by peer",
+        "eof",
+        "socket",
+        "tunnel",
+        "ssh",
+    )
+    return "transport" if any(w in msg for w in transport_words) else "no_data"
 
 
 def _reduce(t: np.ndarray, y: np.ndarray, tmin, tmax, stride: int):
@@ -382,6 +430,7 @@ def _mdsthin_tree_channels(connect, shot, tree_signals, tmin, tmax, progress):
                 conn.openTree(tree, shot)
             except Exception as exc:
                 ch.error = f"open {tree}: {exc}"
+                ch.error_kind = _error_kind(exc)
                 continue
             try:
                 y = np.atleast_1d(conn.get(node).data())
@@ -398,6 +447,7 @@ def _mdsthin_tree_channels(connect, shot, tree_signals, tmin, tmax, progress):
                 ch.error = f"{tree}{node}: degenerate (n={y.size}, t={t.size})"
             except Exception as exc:
                 ch.error = f"{tree}{node}: {exc}"
+                ch.error_kind = _error_kind(exc)
         out.append(ch)
         progress(1.0, f"tree:{name} ({'ok' if ch.ok else 'missing'})")
     return out
@@ -528,7 +578,9 @@ def fetch_batch(conn, batch, use_many):
                         t = _arr(gm.get(f"t{i}"))
                         out.append(Channel(pt, t, y.astype(np.float32, copy=False), ok=True))
                     except Exception as exc:  # per-channel error in the batch
-                        out.append(Channel(pt, ok=False, error=str(exc)))
+                        out.append(
+                            Channel(pt, ok=False, error=str(exc), error_kind=_error_kind(exc))
+                        )
                 tick(batch[-1], k=len(batch))
                 return out, True
             except Exception as exc:  # whole-batch failure -> fall back
@@ -543,7 +595,7 @@ def fetch_batch(conn, batch, use_many):
                 else:
                     out.append(Channel(pt, ok=False, error="no data"))
             except Exception as exc:
-                out.append(Channel(pt, ok=False, error=str(exc)))
+                out.append(Channel(pt, ok=False, error=str(exc), error_kind=_error_kind(exc)))
             tick(pt)
         return out, False
 
@@ -691,6 +743,7 @@ def fetch_all(conn):
                     ch.error = f"degenerate (n={y.size}, t={t.size})"
             except Exception as exc:
                 ch.error = str(exc)
+                ch.error_kind = _error_kind(exc)
             out.append(ch)
             tick(node)
         # plasma / equilibrium tree signals (efit01, ...): first candidate that opens
@@ -714,6 +767,7 @@ def fetch_all(conn):
                     chp.error = f"{tree}{node}: degenerate"
                 except Exception as exc:
                     chp.error = f"{tree}{node}: {exc}"
+                    chp.error_kind = _error_kind(exc)
             out.append(chp)
             progress(1.0, f"tree:{name} ({'ok' if chp.ok else 'missing'})")
         return out
@@ -744,7 +798,9 @@ def _node_tree_map(dev):
     return node_tree, default_tree
 
 
-def _fetch_mds_tree(shot, items, *, connect, tmin, tmax, progress, resample_dt=1e-6):
+def _fetch_mds_tree(
+    shot, items, *, connect, tmin, tmax, progress, resample_dt=1e-6, tree_signals=None
+):
     """Fetch tree nodes one-by-one over a single mdsthin connection.
 
     ``items`` is a list of ``(canonical_id, node, tree, gain)``: open each node's
@@ -759,20 +815,28 @@ def _fetch_mds_tree(shot, items, *, connect, tmin, tmax, progress, resample_dt=1
     ``resample(node, tmin, tmax, dt)`` when both bounds are given: it bounds the wire
     transfer (a raw 2 MHz channel is ~10^7 pts) at a cadence far above rotating-mode
     frequencies (kHz). ``resample_dt`` seconds sets that cadence (default 1 µs =
-    1 MHz). With an open/absent window it falls back to a raw dimension slice."""
+    1 MHz). With an open/absent window it falls back to a raw dimension slice.
 
-    def _slice(expr):
+    ``tree_signals`` — plasma/equilibrium quantities ``{name: [(tree, node), …]}``
+    (Ip, B_T, kappa from the device's "plasma pointnames"): fetched on the same
+    connection after the sensors, first usable candidate wins, mirroring
+    ``_fetch_mdsthin_tree``. Windowed with a raw dimension slice, NOT resample() —
+    these are slow (ms-cadence) traces and resampling them onto the sensors' µs
+    grid would interpolate them ~1000x."""
+
+    def _slice(expr, resample=True):
         if tmin is None and tmax is None:
             return expr
         lo = "*" if tmin is None else repr(float(tmin))
         hi = "*" if tmax is None else repr(float(tmax))
-        if tmin is not None and tmax is not None and resample_dt:
+        if resample and tmin is not None and tmax is not None and resample_dt:
             return f"resample({expr}, {lo}, {hi}, {float(resample_dt)!r})"
         return f"({expr})[{lo} : {hi}]"
 
     conn = connect()
     out: list[Channel] = []
-    n = len(items) or 1
+    tree_signals = tree_signals or {}
+    n = (len(items) + len(tree_signals)) or 1
     current_tree = None
     for i, (canon, node, tree, gain) in enumerate(items, 1):
         try:
@@ -803,8 +867,37 @@ def _slice(expr):
             )
         except Exception as exc:  # noqa: BLE001 — one bad node shouldn't sink the pull
             current_tree = None  # a failed openTree/get can leave the conn unsure
-            out.append(Channel(canon, ok=False, error=str(exc)))
+            out.append(Channel(canon, ok=False, error=str(exc), error_kind=_error_kind(exc)))
         progress(i / n, canon)
+
+    for j, (name, candidates) in enumerate(tree_signals.items(), len(items) + 1):
+        chp = Channel(name, ok=False, error="no data")
+        for tree, node in candidates:
+            try:
+                if tree != current_tree:
+                    conn.openTree(tree, int(shot))
+                    current_tree = tree
+                sliced = _slice(node, resample=False)
+                y = np.atleast_1d(np.asarray(conn.get(sliced).data(), dtype=np.float32))
+                t = np.atleast_1d(np.asarray(conn.get(f"DIM_OF({sliced})").data()))
+                # A slow trace can legitimately hold few samples in a narrow window
+                # (helicity only needs the median sign), so allow n >= 1 here.
+                if y.size < 1 or t.size != y.size or not np.all(np.isfinite(t)):
+                    chp.error = f"{tree}:{node}: degenerate result (n={y.size})"
+                    continue
+                chp = Channel(
+                    name,
+                    (t * 1000.0).astype(np.float64),  # tree time is s; h5 stores ms
+                    y.astype(np.float32, copy=False),
+                    ok=True,
+                )
+                break
+            except Exception as exc:  # noqa: BLE001 — plasma extras must not sink the pull
+                current_tree = None
+                chp.error = f"{tree}:{node}: {exc}"
+                chp.error_kind = _error_kind(exc)
+        out.append(chp)
+        progress(j / n, f"tree:{name} ({'ok' if chp.ok else 'missing'})")
     return out
 
 
@@ -1078,14 +1171,72 @@ def _write_h5(
                     g.attrs[k] = float(v)
 
         # Union new channels with whatever the file already recorded; a name that
-        # is now fetched is removed from the missing list.
+        # is now fetched is removed from the missing list. Transport failures
+        # (dropped tunnel mid-pull, issue #63) are deliberately NOT recorded as
+        # missing — caching them would make every later incremental pull skip
+        # them, leaving the file permanently incomplete without --force.
+        no_data = {c.name for c in missing if getattr(c, "error_kind", "no_data") != "transport"}
         fetched = _attr_names(h5, "channels_fetched") | {c.name for c in got}
-        missing_names = (_attr_names(h5, "channels_missing") | {c.name for c in missing}) - fetched
+        missing_names = (_attr_names(h5, "channels_missing") | no_data) - fetched
         h5.attrs["channels_fetched"] = np.array(sorted(fetched), dtype="S")
         h5.attrs["channels_missing"] = np.array(sorted(missing_names), dtype="S")
     return got, missing
 
 
+def _merge_h5_files(src, dst, *, compression="lzf"):
+    """Fold every channel of shot file ``src`` into the existing shot file ``dst``
+    (same shot/window/decimation — the caller verifies via ``_existing_channels``).
+
+    Channels are re-written through ``_write_h5``'s merge path so time-base dedup,
+    replace-on-conflict, and the fetched/missing bookkeeping behave exactly like a
+    local incremental fetch. Used by the remote backend: the cluster always ships a
+    fresh file holding only the requested selection, which must never replace a
+    local file holding earlier pulls.
+    """
+    import h5py
+
+    channels: list[Channel] = []
+    query_names: dict[str, str] = {}
+    geometry: dict[str, dict] = {}
+    with h5py.File(src, "r") as h5:
+        analysis = h5.attrs.get("analysis", "")
+        analysis = analysis.decode() if isinstance(analysis, bytes) else str(analysis)
+        tmin = h5.attrs.get("tmin")
+        tmax = h5.attrs.get("tmax")
+        stride = int(h5.attrs.get("decimate", 1))
+        tmin = None if isinstance(tmin, (bytes, str)) else tmin
+        tmax = None if isinstance(tmax, (bytes, str)) else tmax
+        for name, g in h5.items():
+            if name == "_timebases" or not isinstance(g, h5py.Group) or "data" not in g:
+                continue
+            channels.append(Channel(name, time=g["time"][()], data=g["data"][()], ok=True))
+            attrs = dict(g.attrs)
+            pn = attrs.pop("pointname", None)
+            if pn is not None:
+                query_names[name] = pn.decode() if isinstance(pn, bytes) else str(pn)
+            geometry[name] = {
+                k: v
+                for k, v in attrs.items()
+                if k != "time_units" and isinstance(v, (int, float, np.integer, np.floating))
+            }
+        for m in sorted(_attr_names(h5, "channels_missing")):
+            channels.append(Channel(m, ok=False, error="missing in remote pull"))
+    return _write_h5(
+        dst,
+        None,  # shot/device attrs already present in dst (merge never rewrites them)
+        analysis,
+        "remote",
+        channels,
+        compression=compression,
+        tmin=tmin,
+        tmax=tmax,
+        stride=stride,
+        query_names=query_names,
+        channel_geometry=geometry,
+        merge=True,
+    )
+
+
 def _existing_channels(path, tmin, tmax, stride):
     """Names already attempted in an existing shot file IF its reduction matches.
 
@@ -1203,31 +1354,10 @@ def fetch_shot(
             progress(0.0, f"{backend}→mdsthin (tree device)")
         backend = "mdsthin"
 
-    if backend == "remote" and not _tree_transport:
-        # Orchestrate a pull on the cluster from here; remote side runs this same
-        # script with --backend toksearch and writes the file we copy back.
-        from . import remote as remote_run
-
-        kw = {} if remote_python is None else {"python": remote_python}
-        return remote_run.run_remote(
-            shot,
-            analysis,
-            host=remote_host,
-            jump=ssh_jump,
-            username=username,
-            password=password,
-            duo=duo,
-            remote_dir=remote_dir,
-            tmin=tmin,
-            tmax=tmax,
-            decimate=decimate,
-            device=device,
-            sensor_set=sensor_set,
-            raw_pointnames=raw_pointnames,
-            local_out_dir=(str(Path(out).parent) if out else None),
-            progress=progress,
-            **kw,
-        )
+    # NOTE: the remote backend is dispatched BELOW, after the output path and the
+    # incremental-merge bookkeeping are resolved — the cluster ships a fresh file
+    # holding only the requested selection, so it must be merged into (never copied
+    # over) an existing local shot file.
 
     # Device config is the source of truth for mdsip addresses; an explicit
     # gateway/server (CLI or caller) overrides the device file. (`dev`, `device_name`,
@@ -1367,6 +1497,44 @@ def fetch_shot(
         sys.stderr.write(f"All {n_skipped} requested signals already in {out}; nothing to fetch.\n")
         return out
 
+    if backend == "remote" and not _tree_transport:
+        # Orchestrate a pull on the cluster from here; remote side runs this same
+        # script with --backend toksearch and writes the file we copy back. When a
+        # compatible local file already exists (merge), the cluster result is staged
+        # beside it and folded in channel-by-channel — copying it over the local file
+        # would silently destroy every previously pulled channel set.
+        from . import remote as remote_run
+
+        kw = {} if remote_python is None else {"python": remote_python}
+        stage_dir = out_path.parent / "_remote_stage" if merge else out_path.parent
+        fetched = remote_run.run_remote(
+            shot,
+            analysis,
+            host=remote_host,
+            jump=ssh_jump,
+            username=username,
+            password=password,
+            duo=duo,
+            remote_dir=remote_dir,
+            tmin=tmin,
+            tmax=tmax,
+            decimate=decimate,
+            device=device,
+            sensor_set=sensor_set,
+            raw_pointnames=raw_pointnames,
+            local_out_dir=str(stage_dir),
+            progress=progress,
+            **kw,
+        )
+        if merge:
+            progress(0.95, f"merging into {out_path.name}")
+            _merge_h5_files(fetched, out)
+            Path(fetched).unlink(missing_ok=True)
+            with contextlib.suppress(OSError):
+                stage_dir.rmdir()
+            return out
+        return fetched
+
     t0 = time.perf_counter()
     if _tree_transport:
         # Tree device reached via its own transport (KSTAR: KFE VPN + nkstar tunnel).
@@ -1407,7 +1575,13 @@ def _connect():
             tmin_s = None if tmin is None else float(tmin) / 1000.0
             tmax_s = None if tmax is None else float(tmax) / 1000.0
             channels = _fetch_mds_tree(
-                shot_i, items, connect=_connect, tmin=tmin_s, tmax=tmax_s, progress=progress
+                shot_i,
+                items,
+                connect=_connect,
+                tmin=tmin_s,
+                tmax=tmax_s,
+                progress=progress,
+                tree_signals=tree_signals,
             )
     elif tree_access:
         # NSTX/NSTX-U: sensors are fastmag tree nodes fetched with a value-window
diff --git a/src/magnetics/service/nodes.py b/src/magnetics/service/nodes.py
index e27028e..0e18e3c 100644
--- a/src/magnetics/service/nodes.py
+++ b/src/magnetics/service/nodes.py
@@ -111,6 +111,39 @@ def tag(names, role):
         tag([n for n, _ in _poloidal_arr(str(shot))], "poloidal array")
     except ValueError:
         pass
+    # The quasi-stationary fit's default sensor set (same resolution as _prep_qs_ds).
+    try:
+        from ..core import qs_device
+
+        dev, _ = _device_key(str(shot))
+        default_cf = "Bp_LFS_midplane"
+        if dev:
+            default_cf = dev.get("qs_default_set") or (
+                "quasi_stationary"
+                if "quasi_stationary" in dev.get("sensor_sets", {})
+                else default_cf
+            )
+        pats = np.atleast_1d(
+            qs_device.resolve_channel_filter(
+                default_cf, h5source.meta(shot).get("device", "DIII-D")
+            )
+        )
+        tag(
+            [nm for nm in all_names if any(re.match(str(p), nm) for p in pats)],
+            "QS fit array",
+        )
+    except Exception:  # noqa: BLE001 — QS not applicable to this shot/device
+        pass
+    # Plasma context channels: not part of any probe array, but consumed by the
+    # analyses — dropping them silently degrades the physics (θ* falls back to
+    # geometric θ; QS helicity falls back to its default) with no error anywhere.
+    for nm, role in (
+        ("kappa", "θ* elongation correction (poloidal nodes)"),
+        ("ip", "QS helicity (sign of Ip·Bt)"),
+        ("bt", "QS helicity (sign of Ip·Bt)"),
+    ):
+        if nm in all_names:
+            tag([nm], role)
 
     used = [{"name": nm, "roles": roles[nm]} for nm in all_names if nm in roles]
     unused = [nm for nm in all_names if nm not in roles]
@@ -126,7 +159,15 @@ def tag(names, role):
 def refresh() -> None:
     """Forget cached state (call after a new fetch writes a file)."""
     h5source.refresh()
-    for fn in (_spec_result, _stack_cached, _array_spectrum, _array_mode_spec, _qs_run, _dev_geom):
+    for fn in (
+        _spec_result,
+        _stack_cached,
+        _array_spectrum,
+        _array_mode_spec,
+        _qs_run,
+        _dev_geom,
+        _real_theta,  # channel→θ map derives from channel_names: stale after a re-pull
+    ):
         fn.cache_clear()
 
 
@@ -339,21 +380,45 @@ def _named_sets(dg, substr: str) -> list[str]:
 
 
 # ── spectrogram: real 2-point MODESPEC cross-spectrogram ─────────────────────
+_PAIR_N_MAX = 5  # the pair must resolve |n| <= this without aliasing (matches the n-map)
+
+
+def _best_pair(arr) -> tuple[tuple[str, float], tuple[str, float]] | None:
+    """Best 2-point probe pair from a (name, phi)-sorted array: the widest *wrapped*
+    separation that still resolves |n| <= _PAIR_N_MAX unaliased. round(phase/Δφ) is
+    exact for all |n| <= N iff (N + 0.5)·|Δφ| <= 180, and within that limit a wider
+    Δφ dilutes phase noise (σ_n = σ_φ/Δφ) — so: widest under the limit; if the array
+    has no such pair, the smallest non-zero separation (best available n range)."""
+    limit = 180.0 / (_PAIR_N_MAX + 0.5)
+    best = smallest = None
+    for i in range(len(arr)):
+        for j in range(i + 1, len(arr)):
+            sep = abs(spectral.wrap_angle_deg(arr[j][1] - arr[i][1]))
+            if sep == 0:
+                continue
+            if sep <= limit and (best is None or sep > best[0]):
+                best = (sep, arr[i], arr[j])
+            if smallest is None or sep < smallest[0]:
+                smallest = (sep, arr[i], arr[j])
+    pick = best or smallest
+    return None if pick is None else (pick[1], pick[2])
+
+
 def _pick_pair(shot) -> tuple[tuple[str, float], tuple[str, float]]:
     """Two toroidally-separated probes for the 2-point cross-spectrogram.
     Prefer the fast Mirnov dB/dt array (DIII-D MPI_BDOT / an NSTX toroidal set / a
     KSTAR declared toroidal array), then integrated Bp. Returns ((name1, phi1),
-    (name2, phi2)) with the widest non-zero separation."""
+    (name2, phi2)) chosen by ``_best_pair`` (aliasing-safe, noise-optimal)."""
     dev, arrays = _arrays(shot)
     if arrays:  # device declares an explicit toroidal set (KSTAR) → use it precisely
-        arr = _set_channels(dev, arrays["toroidal"], shot)
-        if len(arr) >= 2 and arr[0][1] != arr[-1][1]:
-            return arr[0], arr[-1]
+        pair = _best_pair(_set_channels(dev, arrays["toroidal"], shot))
+        if pair is not None:
+            return pair
         raise ValueError("need two toroidally-separated probes for a spectrogram")
     for families in _pick_pair_prefs(shot):
-        arr = _array_channels(shot, families)  # (name, phi), sorted by phi
-        if len(arr) >= 2 and arr[0][1] != arr[-1][1]:
-            return arr[0], arr[-1]
+        pair = _best_pair(_array_channels(shot, families))  # (name, phi), sorted by phi
+        if pair is not None:
+            return pair
     raise ValueError("need two toroidally-separated probes for a spectrogram")
 
 
@@ -394,7 +459,7 @@ def _spec_result(shot: str, slice_duration: float, coherence_smooth: int, max_co
         coherence_smooth=coherence_smooth,
         max_columns=max_columns,
     )
-    return res, (n1, n2), round(float(phi2 - phi1), 1)
+    return res, (n1, n2), round(spectral.wrap_angle_deg(phi2 - phi1), 1)
 
 
 def _prep_spec(shot, params):
@@ -503,7 +568,12 @@ def _mode_number(shot, params=None) -> dict:
     names = tuple(n for n, _ in arr)
     phis = tuple(float(p) for _, p in arr)
     sd = _f(params, "slice_duration", 0.002)  # honor the resolution knob (500 Hz default)
-    ms = _array_mode_spec(str(shot), names, phis, sd)
+    # Compute band follows the requested display band, quantized to 50-kHz steps so
+    # crops within a step reuse the cached STFT (issue #85: the ceiling was pinned at
+    # 50 kHz, so raising the GUI's f_max never showed data above it).
+    fmax_req = _f(params, "fmax")
+    band_hi = 50_000.0 * max(1.0, float(np.ceil((fmax_req or 50.0) / 50.0)))
+    ms = _array_mode_spec(str(shot), names, phis, sd, band_hi)
 
     # Optional 2-D Gaussian pre-smoothing on the NATIVE n-map grid (before the display-time
     # decimation below). σ is in grid cells, so a coherent mode survives the gate across its
@@ -619,12 +689,13 @@ def _toroidal_grid(shot):
 # ── contour: raw δBp(φ, t) — x=φ, y=time (QS contour hero plot) ──────────────
 def _contour(shot, params=None) -> dict:
     t_sub, phi_grid, z, phis, n_ch = _toroidal_grid(shot)
-    zmax = float(np.nanmax(np.abs(z))) or 1.0
+    finite = np.abs(z[np.isfinite(z)])
+    zmax = float(finite.max()) if finite.size and finite.max() > 0 else 1.0
     overlay = {"points": [{"x": float(p), "y": float(t_sub[0])} for p in phis], "symbol": "square"}
     return contracts.contour(
         phi_grid.tolist(),
         t_sub.tolist(),
-        z.tolist(),
+        contracts.json_finite(z),
         {"x": "φ (deg)", "y": "time (ms)", "z": "δBp (G)"},
         zrange=[-zmax, zmax],
         overlay=overlay,
@@ -674,7 +745,7 @@ def _toroidal_stripes(shot, params=None) -> dict:
     return contracts.heatmap(
         t_sub.tolist(),
         ang.tolist(),
-        z.tolist(),
+        contracts.json_finite(z),
         {"x": "time (ms)", "y": "φ (deg)", "z": "δBp (a.u.)"},
         discrete=False,
         meta={
@@ -698,7 +769,7 @@ def _poloidal_stripes(shot, params=None) -> dict:
     return contracts.heatmap(
         t_sub.tolist(),
         ang.tolist(),
-        z.tolist(),
+        contracts.json_finite(z),
         {"x": "time (ms)", "y": "θ (deg)", "z": "δBp (a.u.)"},
         discrete=False,
         meta={
@@ -730,7 +801,7 @@ def _raw_trace(shot, params=None) -> dict:
         sel = np.arange(max(0, c - 2000), min(t_ms.size, c + 2000))
     if sel.size > 2000:  # keep the line light
         sel = sel[np.linspace(0, sel.size - 1, 2000).astype(int)]
-    series = [{"name": name, "x": t_ms[sel].tolist(), "y": d[sel].tolist()}]
+    series = [{"name": name, "x": t_ms[sel].tolist(), "y": contracts.json_finite(d[sel])}]
     return contracts.line(
         series,
         {"x": "time (ms)", "y": "dB/dt (a.u.)"},
@@ -833,15 +904,18 @@ def _array_spectrum(shot, names):
 
 
 @lru_cache(maxsize=2)
-def _array_mode_spec(shot, names, phis, slice_duration):
+def _array_mode_spec(shot, names, phis, slice_duration, band_hi_hz=50_000.0):
     """Toroidal-|n|-resolved spectrogram from a *dedicated* full-array STFT computed at
-    the requested frequency resolution over a fixed 0–50 kHz band — so the n-map refines
-    with the resolution knob and spans the same band as the power view, rather than being
-    locked to the cursor-analysis spectrum's 1–25 kHz / 1 kHz grid.
+    the requested frequency resolution over 0–``band_hi_hz`` — so the n-map refines
+    with the resolution knob and follows the power view's band, rather than being
+    locked to the cursor-analysis spectrum's 1–25 kHz / 1 kHz grid. The caller
+    quantizes ``band_hi_hz`` to 50-kHz steps, so band crops *within* a step stay a
+    cheap post-op on the cached STFT (issue #85: a fixed 50 kHz ceiling silently
+    truncated the map no matter how far the GUI's f_max was raised).
 
     Heavier than reusing ``_array_spectrum`` (its own 14-probe STFT), so ``maxsize`` is
-    small and only the resolution changes it — band cropping is a cheap post-op in the
-    node. ``names``/``phis`` are tuples so the call is hashable; cleared by ``refresh``."""
+    small and only the resolution/band change it. ``names``/``phis`` are tuples so the
+    call is hashable; cleared by ``refresh``."""
     t_ms, mat = _stack(shot, names)
     # Cap STFT columns at ~2000: the node decimates the display to 1500 anyway, so the
     # native fine hop would just inflate the per-cell projection (an (n, t, f) einsum)
@@ -850,7 +924,7 @@ def _array_mode_spec(shot, names, phis, slice_duration):
         mat,
         np.asarray(t_ms, dtype=float) * 1e-3,
         fmin=0.0,
-        fmax=50_000.0,
+        fmax=float(band_hi_hz),
         slice_duration=slice_duration,
         max_columns=2000,
     )
@@ -860,7 +934,7 @@ def _array_mode_spec(shot, names, phis, slice_duration):
 # ── toroidal mode at one frequency/cursor (shared by phase_fit & mode_shape) ──
 def _toroidal_arr(shot):
     """The toroidal (midplane) array for the n-fit: ONE consistent probe type, all at
-    θ≈0 so the phase is a clean −nφ ramp. Prefer the fast-Mirnov dB/dt array; a "both"
+    θ≈0 so the phase is a clean +nφ ramp. Prefer the fast-Mirnov dB/dt array; a "both"
     pull also brings the integrated-Bp (MPID) family and the off-midplane *poloidal*
     probes — mixing those in (different units, 90° dB/dt-vs-B offset, m·θ dependence)
     scrambles the fit, so they're excluded."""
@@ -909,12 +983,13 @@ def _toroidal_mode(shot, params):
 
 # ── phase_fit: phase-vs-φ at one frequency, at the GUI time cursor ────────────
 def _wrapped_ramp(intercept_deg, slope_n) -> dict:
-    """Fitted line phase(a) = (c − n·a) mod 360 over a∈[0,360], WRAPPED so it traces
-    the same |n| sawteeth as the (wrapped) data instead of one line shooting off-axis.
-    A null break is inserted at each 0/360 wrap so the polyline doesn't draw a vertical
+    """Fitted line phase(a) = (c + n·a) mod 360 over a∈[0,360] (the conj(sig)·ref
+    cross-phase ramp ``fit_toroidal_mode`` fits), WRAPPED so it traces the same |n|
+    sawteeth as the (wrapped) data instead of one line shooting off-axis. A null
+    break is inserted at each 0/360 wrap so the polyline doesn't draw a vertical
     jump across the panel."""
     a = np.linspace(0.0, 360.0, 361)
-    y = (intercept_deg - slope_n * a) % 360.0
+    y = (intercept_deg + slope_n * a) % 360.0
     fx: list = []
     fy: list = []
     prev = None
@@ -1081,7 +1156,7 @@ def _toroidal_n(shot, t0_s, f_khz):
 
 def _poloidal_mode(shot, params):
     """Per-probe phase/amplitude across the poloidal array vs θ. The probes span φ, so
-    the toroidal nφ ramp is removed (phase += n·φ → −m·θ + const) using the toroidal n
+    the toroidal +nφ ramp is removed (phase −= n·φ → m·θ + const) using the toroidal n
     at the same (t0, f). The probe angles are mapped to the elongation-corrected θ*
     (using the EFIT κ at the cursor) when available, so an `m` mode is a clean sinusoid
     rather than a κ-distorted one. Returns (arr, mode, f_khz, t0_ms, kappa)."""
@@ -1097,7 +1172,7 @@ def _poloidal_mode(shot, params):
     spec = _array_spectrum(str(shot), tuple(n for n, _ in arr))
     t0_s = (t0_ms * 1e-3) if t0_ms is not None else float(spec.time[spec.time.size // 2])
     mode = spectral.mode_from_spectrum(spec, thetas, t0_s, f_khz * 1e3)
-    detrended = (mode.phase + _toroidal_n(str(shot), t0_s, f_khz) * pphis) % 360.0
+    detrended = (mode.phase - _toroidal_n(str(shot), t0_s, f_khz) * pphis) % 360.0
     mode = dataclasses.replace(mode, phase=detrended)
     return arr, mode, f_khz, t0_ms, kappa
 
@@ -1421,7 +1496,9 @@ def _prep_qs_ds(shot, params):
     cutoff_hi = float(params.get("cutoff_hi", 250.0)) if params else 250.0
     energy = float(params.get("energy", 0.98)) if params else 0.98
     fit_basis = params.get("fit_basis", "sinusoidal-integral") if params else "sinusoidal-integral"
-    fit_cond = float(params.get("fit_cond", 10.0)) if params else 10.0
+    # 1e3 = OMFIT SLCONTOUR's inversion cutoff (1/rcond); the "warn when K > 10"
+    # trust threshold lives in contracts.quality_for_k, not here.
+    fit_cond = float(params.get("fit_cond", 1e3)) if params else 1e3
     sigma_str = params.get("sigma") if params else None
     sigma = float(sigma_str) if sigma_str is not None else None
 
diff --git a/tests/conftest.py b/tests/conftest.py
index a945cd6..55607a5 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -127,14 +127,16 @@ def synthetic_n2():
     t = np.linspace(0, duration, int(fs * duration), endpoint=False)
     f_mode = 3_000.0
     n = 2
-    phi1, phi2 = 30.0, 63.0  # Δφ = -33°, mimics the real 307/340 pair separation
+    phi1, phi2 = 30.0, 63.0  # Δφ = +33°, mimics the real 307/340 pair separation
     sig1 = np.sin(2 * np.pi * f_mode * t - np.deg2rad(n * phi1))
     sig2 = np.sin(2 * np.pi * f_mode * t - np.deg2rad(n * phi2))
     return {
         "time": t,
         "sig1": sig1,
         "sig2": sig2,
-        "delta_phi": phi1 - phi2,
+        # φ2 − φ1: the service convention (nodes._spec_result), so estimators
+        # recover the SIGNED n of a cos(nφ − ωt) mode as +n.
+        "delta_phi": phi2 - phi1,
         "fs": float(fs),
         "f_mode": f_mode,
         "n_true": n,
diff --git a/tests/test_cli.py b/tests/test_cli.py
index 566cfef..e789fa2 100644
--- a/tests/test_cli.py
+++ b/tests/test_cli.py
@@ -35,7 +35,7 @@ def test_help_exits_clean(capsys):
     with pytest.raises(SystemExit) as exc:
         cli.main(["--help"])
     assert exc.value.code == 0
-    assert "magnetics GUI" in capsys.readouterr().out or True  # help text printed
+    assert "magnetics GUI" in capsys.readouterr().out  # help text printed
 
 
 def test_data_dir_env_override(tmp_path, monkeypatch):
diff --git a/tests/test_contract_meta.py b/tests/test_contract_meta.py
index 35a03e6..a1bb8f7 100644
--- a/tests/test_contract_meta.py
+++ b/tests/test_contract_meta.py
@@ -33,12 +33,17 @@ def test_signal_conditioning_meta_pairs(synthetic_shot):
 
 def test_amplitude_meta_fields(synthetic_shot):
     meta = _meta(synthetic_shot, "amplitude")
-    assert "sigma" in meta and "legend_title" in meta
+    assert "legend_title" in meta
+    # the ±1σ band rides in the typed per-series lower/upper, not meta.sigma
+    node = nodes.build_node(synthetic_shot, "amplitude")
+    assert all("lower" in s and "upper" in s for s in node["series"])
 
 
 def test_phase_t_meta_fields(synthetic_shot):
     meta = _meta(synthetic_shot, "phase_t")
-    assert "sigma" in meta and "phase_visible" in meta
+    assert "phase_visible" in meta
+    node = nodes.build_node(synthetic_shot, "phase_t")
+    assert all("lower" in s and "upper" in s for s in node["series"])
 
 
 def test_phase_fit_reports_n_estimate(synthetic_shot):
diff --git a/tests/test_contract_shapes.py b/tests/test_contract_shapes.py
index 6eb8c54..f98994b 100644
--- a/tests/test_contract_shapes.py
+++ b/tests/test_contract_shapes.py
@@ -1,7 +1,6 @@
 """Smoke tests for the mock seam: every generator emits the CONTRACT.md shapes.
 
-Imports only magnetics.service.mock (numpy), not the FastAPI app — so these run
-with the core install, no `service` extra needed.
+Imports only magnetics.service.mock (numpy), not the FastAPI app.
 """
 
 from magnetics.service import mock
diff --git a/tests/test_fetch_api.py b/tests/test_fetch_api.py
new file mode 100644
index 0000000..7aca6e1
--- /dev/null
+++ b/tests/test_fetch_api.py
@@ -0,0 +1,224 @@
+"""The live-pull job machinery: POST /api/fetch → background worker → status/SSE.
+
+This is the most stateful, most concurrent code in the service (job registry +
+worker thread + progress stream) and the only HTTP surface the GUI's PullControl
+talks to. The real fetcher is monkeypatched — no network, no data files written —
+so these tests pin the ORCHESTRATION: request forwarding, job lifecycle, error
+mapping, cache refresh, and the SSE frame protocol.
+"""
+
+from __future__ import annotations
+
+import json
+import threading
+import time
+
+import pytest
+from fastapi.testclient import TestClient
+
+from magnetics.data.fetch import toksearch
+from magnetics.service import nodes
+from magnetics.service.app import app
+
+client = TestClient(app)
+
+
+def _wait_for(jid, *, until=("done", "error"), timeout=5.0):
+    """Poll GET /api/fetch/{jid} until the job leaves 'running' (or timeout)."""
+    deadline = time.monotonic() + timeout
+    while time.monotonic() < deadline:
+        r = client.get(f"/api/fetch/{jid}")
+        assert r.status_code == 200
+        job = r.json()
+        if job["status"] in until:
+            return job
+        time.sleep(0.02)
+    raise AssertionError(f"job {jid} still running after {timeout}s: {job}")
+
+
+def test_post_fetch_happy_path_forwards_request_and_refreshes(monkeypatch, synthetic_shot):
+    calls = {}
+    refreshed = {"nodes": 0}
+
+    def fake_fetch_shot(shot, analysis, *, progress, **kw):
+        calls["shot"], calls["analysis"], calls["kw"] = shot, analysis, kw
+        progress(0.5, "pulling")
+        return "/tmp/shot_123456.h5"
+
+    real_refresh = nodes.refresh
+
+    def counting_refresh():
+        refreshed["nodes"] += 1
+        real_refresh()
+
+    monkeypatch.setattr(toksearch, "fetch_shot", fake_fetch_shot)
+    monkeypatch.setattr(nodes, "refresh", counting_refresh)
+
+    r = client.post(
+        "/api/fetch",
+        json={
+            "shot": 123456,
+            "analysis": "rotating",
+            "backend": "mdsthin",
+            "tmin": 1000,
+            "tmax": 2000,
+            "device": "diiid",
+            "sensor_set": "Bp LFS midplane",
+        },
+    )
+    assert r.status_code == 200
+    jid = r.json()["job_id"]
+
+    job = _wait_for(jid)
+    assert job["status"] == "done"
+    assert job["progress"] == 1.0
+    # result contract PullControl relies on: shot is a STRING id + fresh machine list
+    assert job["result"]["ok"] is True
+    assert job["result"]["shot"] == "123456"
+    assert job["result"]["file"] == "/tmp/shot_123456.h5"
+    assert any(m["id"] == synthetic_shot for m in job["result"]["machines"])
+    # the request reached the fetcher with the caller's selection...
+    assert calls["shot"] == 123456 and calls["analysis"] == "rotating"
+    assert calls["kw"]["backend"] == "mdsthin"
+    assert calls["kw"]["tmin"] == 1000 and calls["kw"]["tmax"] == 2000
+    assert calls["kw"]["device"] == "diiid"
+    assert calls["kw"]["sensor_set"] == "Bp LFS midplane"
+    # ...unset optionals were NOT passed (the fetcher's own defaults must apply)
+    assert "raw_pointnames" not in calls["kw"]
+    # a successful pull must invalidate the node caches (stale-θ regression family)
+    assert refreshed["nodes"] == 1
+
+
+def test_post_fetch_custom_signals_become_raw_pointnames(monkeypatch):
+    seen = {}
+
+    def fake_fetch_shot(shot, analysis, *, progress, **kw):
+        seen.update(kw)
+        return "/tmp/x.h5"
+
+    monkeypatch.setattr(toksearch, "fetch_shot", fake_fetch_shot)
+    r = client.post("/api/fetch", json={"shot": 1, "signals": ["betan", "li"]})
+    job = _wait_for(r.json()["job_id"])
+    assert job["status"] == "done"
+    assert seen["raw_pointnames"] == ["betan", "li"]
+
+
+@pytest.mark.parametrize(
+    "raiser, expect",
+    [
+        (lambda: (_ for _ in ()).throw(RuntimeError("tunnel died")), "fetch failed: tunnel died"),
+        (
+            lambda: (_ for _ in ()).throw(SystemExit("Missing dependency: mdsthin")),
+            "Missing dependency: mdsthin",
+        ),
+    ],
+)
+def test_post_fetch_maps_failures_to_job_error(monkeypatch, raiser, expect):
+    def fake_fetch_shot(shot, analysis, *, progress, **kw):
+        next(raiser())
+
+    monkeypatch.setattr(toksearch, "fetch_shot", fake_fetch_shot)
+    r = client.post("/api/fetch", json={"shot": 1})
+    job = _wait_for(r.json()["job_id"])
+    assert job["status"] == "error"
+    assert expect in job["error"]
+    assert job["result"] is None
+
+
+def test_fetch_status_unknown_job_404s():
+    assert client.get("/api/fetch/nonexistent0").status_code == 404
+    assert client.get("/api/fetch/nonexistent0/stream").status_code == 404
+
+
+def test_progress_callback_updates_job_live(monkeypatch):
+    """The worker's progress callback must be visible through GET /api/fetch/{jid}
+    while the pull is still running — this is the GUI's per-channel progress bar."""
+    release = threading.Event()
+
+    def fake_fetch_shot(shot, analysis, *, progress, **kw):
+        progress(0.42, "MPI66M307D")
+        assert release.wait(5.0), "test never released the worker"
+        return "/tmp/x.h5"
+
+    monkeypatch.setattr(toksearch, "fetch_shot", fake_fetch_shot)
+    r = client.post("/api/fetch", json={"shot": 1})
+    jid = r.json()["job_id"]
+    try:
+        deadline = time.monotonic() + 5.0
+        while time.monotonic() < deadline:
+            job = client.get(f"/api/fetch/{jid}").json()
+            if job["progress"] == 0.42:
+                break
+            time.sleep(0.02)
+        assert job["status"] == "running"
+        assert job["progress"] == 0.42 and job["msg"] == "MPI66M307D"
+    finally:
+        release.set()
+    assert _wait_for(jid)["status"] == "done"
+
+
+def test_stream_emits_terminal_sse_frame(monkeypatch):
+    """The SSE stream must emit `data: <json>` frames and terminate with the
+    done-frame carrying the result — the exact protocol PullControl's EventSource
+    parses."""
+
+    def fake_fetch_shot(shot, analysis, *, progress, **kw):
+        return "/tmp/x.h5"
+
+    monkeypatch.setattr(toksearch, "fetch_shot", fake_fetch_shot)
+    r = client.post("/api/fetch", json={"shot": 77})
+    jid = r.json()["job_id"]
+    _wait_for(jid)  # let it finish so the stream is a single terminal frame
+
+    frames = []
+    with client.stream("GET", f"/api/fetch/{jid}/stream") as resp:
+        assert resp.status_code == 200
+        assert resp.headers["content-type"].startswith("text/event-stream")
+        for line in resp.iter_lines():
+            if line.startswith("data: "):
+                frames.append(json.loads(line[len("data: ") :]))
+    assert frames, "stream produced no data frames"
+    last = frames[-1]
+    assert last["status"] == "done"
+    assert last["progress"] == 1.0
+    assert last["result"]["shot"] == "77"
+
+
+def test_stream_error_frame_carries_the_reason(monkeypatch):
+    def fake_fetch_shot(shot, analysis, *, progress, **kw):
+        raise RuntimeError("boom")
+
+    monkeypatch.setattr(toksearch, "fetch_shot", fake_fetch_shot)
+    r = client.post("/api/fetch", json={"shot": 1})
+    jid = r.json()["job_id"]
+    _wait_for(jid, until=("error",))
+    with client.stream("GET", f"/api/fetch/{jid}/stream") as resp:
+        line = next(ln for ln in resp.iter_lines() if ln.startswith("data: "))
+    frame = json.loads(line[len("data: ") :])
+    assert frame["status"] == "error"
+    assert "boom" in frame["error"]
+
+
+def test_post_fetch_validation_rejects_missing_shot():
+    assert client.post("/api/fetch", json={"analysis": "rotating"}).status_code == 422
+
+
+def test_job_registry_is_isolated_per_job(monkeypatch):
+    """Two jobs must not share state (the registry is a plain dict guarded by a
+    lock; a key collision or shared entry would cross progress between pulls)."""
+
+    def ok(shot, analysis, *, progress, **kw):
+        return f"/tmp/shot_{shot}.h5"
+
+    def bad(shot, analysis, *, progress, **kw):
+        raise RuntimeError("nope")
+
+    monkeypatch.setattr(toksearch, "fetch_shot", ok)
+    j1 = client.post("/api/fetch", json={"shot": 1}).json()["job_id"]
+    _wait_for(j1)
+    monkeypatch.setattr(toksearch, "fetch_shot", bad)
+    j2 = client.post("/api/fetch", json={"shot": 2}).json()["job_id"]
+    _wait_for(j2, until=("error",))
+    assert j1 != j2
+    assert client.get(f"/api/fetch/{j1}").json()["status"] == "done"
+    assert client.get(f"/api/fetch/{j2}").json()["status"] == "error"
diff --git a/tests/test_kstar_routing.py b/tests/test_kstar_routing.py
index 7e864bf..8decc12 100644
--- a/tests/test_kstar_routing.py
+++ b/tests/test_kstar_routing.py
@@ -4,9 +4,11 @@
 pull was sent to the GA `omega` ssh alias and failed with "SSH to the cluster failed".
 """
 
+import numpy as np
 import pytest
 
 from magnetics.data.fetch import kstar_transport, remote, toksearch
+from magnetics.data.fetch.toksearch import _fetch_mds_tree, _resolve_pointnames
 
 
 class _TransportHit(Exception):
@@ -69,3 +71,107 @@ def fake_remote(*a, **k):
         force=True,
     )
     assert hit["remote"] is True, "DIII-D remote pull no longer reaches run_remote"
+
+
+# ---------------------------------------------------------------------------
+# tree_signals (Ip / B_T / kappa) must be fetched on the KSTAR tree path, and a
+# per-era pointname alias that collides with another sensor must not silently
+# drop a channel. Both run fully offline against stub connections.
+# ---------------------------------------------------------------------------
+
+
+class _DefaultConn:
+    """Stub mdsthin connection: every get() returns a valid (data, DIM_OF) pair
+    unless the expression is in `fail` (then it raises). Records expressions."""
+
+    def __init__(self, fail=()):
+        self.fail = set(fail)
+        self.exprs = []
+        self.opened = []
+
+    def openTree(self, tree, shot):  # noqa: N802 (MDSplus API name)
+        self.opened.append(tree)
+
+    def get(self, expr):
+        self.exprs.append(expr)
+        if any(f in expr for f in self.fail):
+            raise RuntimeError(f"TreeNNF: {expr}")
+
+        class _R:
+            @staticmethod
+            def data():
+                if expr.startswith("DIM_OF("):
+                    return np.linspace(0.25, 0.35, 64)
+                return np.ones(64, np.float32)
+
+        return _R()
+
+
+def test_tree_signals_are_fetched_on_the_kstar_path():
+    """Regression: the _tree_transport branch never passed tree_signals to
+    _fetch_mds_tree, so Ip/B_T/kappa were silently dropped for every KSTAR shot —
+    absent from channels_fetched AND channels_missing — and helicity fell back to
+    its default."""
+    conn = _DefaultConn()
+    channels = _fetch_mds_tree(
+        12345,
+        [("\\MC1T01", "\\MC1T01", "kstar", 1.0)],
+        connect=lambda: conn,
+        tmin=0.25,
+        tmax=0.35,
+        progress=lambda *_: None,
+        tree_signals={
+            "Ip": [("PCS_KSTAR", "\\PC:PCRC03*(-1)*(1e-6)")],
+            "kappa": [("efitrt1", "\\kappa"), ("efitrt1", "\\top.results.aeqdsk:kappa")],
+        },
+    )
+    by_name = {c.name: c for c in channels}
+    assert by_name["Ip"].ok and by_name["kappa"].ok
+    # tree time is seconds; stored time must be ms
+    np.testing.assert_allclose(by_name["Ip"].time[0], 250.0)
+    # plasma traces are window-SLICED, never resample()d onto the sensor µs grid
+    ip_exprs = [e for e in conn.exprs if "PCRC03" in e and not e.startswith("DIM_OF")]
+    assert ip_exprs == ["(\\PC:PCRC03*(-1)*(1e-6))[0.25 : 0.35]"]
+    # the sensor still goes through resample()
+    assert any(e.startswith("resample(\\MC1T01") for e in conn.exprs)
+
+
+def test_tree_signal_candidate_fallback_and_missing_are_reported():
+    conn = _DefaultConn(fail=("\\kappa",))  # first kappa candidate fails
+    channels = _fetch_mds_tree(
+        12345,
+        [],
+        connect=lambda: conn,
+        tmin=None,
+        tmax=None,
+        progress=lambda *_: None,
+        tree_signals={
+            "kappa": [("efitrt1", "\\kappa"), ("efitrt1", "\\top.results.aeqdsk:kappa")],
+            "B_T": [("PCS_KSTAR", "\\kappa_not_there")],  # also matches fail
+        },
+    )
+    by_name = {c.name: c for c in channels}
+    assert by_name["kappa"].ok  # second candidate won
+    assert not by_name["B_T"].ok  # recorded as missing, not lost
+    assert "B_T" in {c.name for c in channels if not c.ok}
+
+
+def test_alias_collision_is_skipped_and_identity_wins():
+    """Regression: kstar.json maps \\MC1T10 -> pointname \\MC1P03 (since shot 17377)
+    while \\MC1P03 is itself in the pull; both fetched channels got relabeled to the
+    same name and the HDF5 writer silently dropped one. The node must be queried
+    once, the sensor literally named \\MC1P03 keeps it, and the alias is reported
+    as skipped — in either request order."""
+    dev = {
+        "sensors": {
+            "\\MC1P03": {"segments": [{"since_shot": 0, "phi": 10.0}]},
+            "\\MC1T10": {
+                "segments": [{"since_shot": 17377, "phi": 160.5, "pointname": "\\MC1P03"}]
+            },
+        }
+    }
+    for order in (["\\MC1T10", "\\MC1P03"], ["\\MC1P03", "\\MC1T10"]):
+        query, canon, skipped = _resolve_pointnames(dev, order, 20000)
+        assert query == ["\\MC1P03"], order
+        assert canon == {"\\MC1P03": "\\MC1P03"}, order
+        assert skipped == ["\\MC1T10"], order
diff --git a/tests/test_mdsthin_fetch.py b/tests/test_mdsthin_fetch.py
new file mode 100644
index 0000000..f7e37fe
--- /dev/null
+++ b/tests/test_mdsthin_fetch.py
@@ -0,0 +1,210 @@
+"""Offline tests for the DIII-D mdsthin fetch body (`_fetch_mdsthin`).
+
+The getMany batching, its per-channel fallback, and the transport-vs-no-data
+error classification are the guts of every laptop pull, but they only ran
+against a live gateway before. A stub `mdsthin.Connection` (per-pointname
+scripted responses) exercises them without a network; `tcp=True` skips the SSH
+tunnel so no subprocess is spawned.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+import mdsthin
+
+from magnetics.data.fetch import toksearch
+
+T = np.linspace(0.0, 10.0, 50)
+Y = np.ones(50, np.float32)
+
+
+class _Desc:
+    def __init__(self, arr):
+        self._a = np.asarray(arr)
+
+    def data(self):
+        return self._a
+
+
+class _StubGetMany:
+    def __init__(self, conn):
+        self.conn = conn
+        self.named: dict[str, tuple[str, str]] = {}  # d0/t0 -> (pointname, expr)
+
+    def append(self, name, expr, pt, shot):
+        self.named[name] = (pt, expr)
+        self.conn.exprs.append(expr)
+
+    def execute(self):
+        self.conn.getmany_execs += 1
+        if self.conn.fail_batch:
+            raise RuntimeError("getMany not supported by this server")
+
+    def get(self, name):
+        pt, _ = self.named[name]
+        kind, *payload = self.conn.responses.get(pt, ("nodata",))
+        if kind == "raise":
+            raise payload[0]
+        if kind == "nodata":
+            return _Desc([0])  # ptdata2's no-data sentinel: length-1 [0]
+        data, time = payload
+        return _Desc(data if name.startswith("d") else time)
+
+
+class _StubConn:
+    """Scripted mdsthin connection. responses: pointname -> ("ok", data, time) |
+    ("nodata",) | ("raise", exc). `fail_batch=True` makes getMany.execute raise,
+    forcing the per-channel fallback."""
+
+    def __init__(self, responses, fail_batch=False):
+        self.responses = responses
+        self.fail_batch = fail_batch
+        self.exprs: list[str] = []
+        self.getmany_execs = 0
+        self.per_channel_calls = 0
+        self._last_time = None
+
+    def getMany(self):  # noqa: N802 (MDSplus API name)
+        return _StubGetMany(self)
+
+    def get(self, expr, *args):
+        self.exprs.append(expr)
+        if expr.startswith("_s = "):
+            self.per_channel_calls += 1
+            pt = args[0]
+            kind, *payload = self.responses.get(pt, ("nodata",))
+            if kind == "raise":
+                raise payload[0]
+            if kind == "nodata":
+                self._last_time = None
+                return _Desc([0])
+            data, time = payload
+            self._last_time = time
+            return _Desc(data)
+        if expr == "dim_of(_s)":
+            return _Desc(self._last_time)
+        raise AssertionError(f"unexpected TDI expr: {expr}")
+
+
+def _fetch(monkeypatch, conn, pointnames, **kw):
+    monkeypatch.setattr(mdsthin, "Connection", lambda addr: conn)
+    defaults = dict(
+        username="user",
+        gateway=None,
+        server="atlas.gat.com:8000",
+        tcp=True,  # direct dial → no SSH tunnel subprocess
+        tmin=1000.0,
+        tmax=2000.0,
+        stride=1,
+        workers=1,
+        batch_size=2,
+        progress=lambda *a: None,
+    )
+    defaults.update(kw)
+    return toksearch._fetch_mdsthin(184927, pointnames, **defaults)
+
+
+def test_getmany_happy_path_batches_and_orders(monkeypatch):
+    conn = _StubConn({p: ("ok", Y * i, T) for i, p in enumerate(["A", "B", "C"], 1)})
+    out = _fetch(monkeypatch, conn, ["A", "B", "C"])
+    # deterministic caller order regardless of batch/worker completion order
+    assert [c.name for c in out] == ["A", "B", "C"]
+    assert all(c.ok for c in out)
+    np.testing.assert_array_equal(out[1].data, Y * 2)
+    np.testing.assert_array_equal(out[0].time, T)
+    # 3 channels at batch_size=2 → 2 getMany round trips, no per-channel calls
+    assert conn.getmany_execs == 2
+    assert conn.per_channel_calls == 0
+    # the server-side reduction suffix is on every expression (window pushed remote)
+    assert all("[1000.0 : 2000.0]" in e for e in conn.exprs if "ptdata2" in e)
+
+
+def test_stride_appends_index_subscript(monkeypatch):
+    conn = _StubConn({"A": ("ok", Y, T)})
+    _fetch(monkeypatch, conn, ["A"], stride=4)
+    assert any("[0 : * : 4]" in e for e in conn.exprs), conn.exprs
+
+
+def test_no_data_channel_is_missing_with_no_data_kind(monkeypatch):
+    conn = _StubConn({"A": ("ok", Y, T), "B": ("nodata",)})
+    out = {c.name: c for c in _fetch(monkeypatch, conn, ["A", "B"])}
+    assert out["A"].ok
+    assert not out["B"].ok
+    assert out["B"].error == "no data"
+    assert out["B"].error_kind == "no_data"  # cacheable as missing
+
+
+def test_per_channel_error_in_batch_is_classified(monkeypatch):
+    conn = _StubConn(
+        {
+            "A": ("ok", Y, T),
+            "B": ("raise", RuntimeError("Connection reset by peer")),
+            "C": ("raise", RuntimeError("%TREE-W-NODATA, no data available")),
+        }
+    )
+    out = {c.name: c for c in _fetch(monkeypatch, conn, ["A", "B", "C"], batch_size=3)}
+    assert out["A"].ok
+    # transport failure → retried next pull (issue #63); MDSplus no-data → cached
+    assert out["B"].error_kind == "transport"
+    assert out["C"].error_kind == "no_data"
+
+
+def test_whole_batch_failure_falls_back_to_per_channel(monkeypatch):
+    conn = _StubConn({p: ("ok", Y, T) for p in ["A", "B", "C"]}, fail_batch=True)
+    out = _fetch(monkeypatch, conn, ["A", "B", "C"])
+    assert [c.name for c in out] == ["A", "B", "C"]
+    assert all(c.ok for c in out)
+    # one failed getMany attempt, then everything through the 2-round-trip path —
+    # and the fallback sticks for later batches (no second getMany attempt)
+    assert conn.getmany_execs == 1
+    assert conn.per_channel_calls == 3
+
+
+def test_per_channel_mode_never_tries_getmany(monkeypatch):
+    conn = _StubConn({"A": ("ok", Y, T)})
+    out = _fetch(monkeypatch, conn, ["A"], per_channel=True)
+    assert out[0].ok
+    assert conn.getmany_execs == 0
+    assert conn.per_channel_calls == 1
+
+
+def test_progress_reaches_one_across_workers(monkeypatch):
+    fracs = []
+    conn = _StubConn({p: ("ok", Y, T) for p in "ABCDEF"})
+    _fetch(
+        monkeypatch,
+        conn,
+        list("ABCDEF"),
+        workers=3,
+        batch_size=1,
+        progress=lambda f, m: fracs.append(f),
+    )
+    assert fracs and max(fracs) == pytest.approx(1.0)
+
+
+def test_transport_failure_end_to_end_is_not_cached_as_missing(monkeypatch, tmp_path):
+    """The #63 chain through the real writer: a batch whose connection died must
+    leave the channel retryable — in the run report but NOT in channels_missing."""
+    import h5py
+
+    conn = _StubConn({"A": ("ok", Y, T), "B": ("raise", ConnectionResetError("reset"))})
+    channels = _fetch(monkeypatch, conn, ["A", "B"])
+    out = str(tmp_path / "shot.h5")
+    toksearch._write_h5(
+        out,
+        184927,
+        "rotating",
+        "mdsthin",
+        channels,
+        compression="lzf",
+        tmin=1000.0,
+        tmax=2000.0,
+        stride=1,
+    )
+    with h5py.File(out, "r") as h5:
+        missing = {x.decode() for x in h5.attrs["channels_missing"]}
+        fetched = {x.decode() for x in h5.attrs["channels_fetched"]}
+    assert fetched == {"A"}
+    assert missing == set()  # B is retryable, not poisoned into the cache
diff --git a/tests/test_mode_shape.py b/tests/test_mode_shape.py
index 9d73086..a671476 100644
--- a/tests/test_mode_shape.py
+++ b/tests/test_mode_shape.py
@@ -152,9 +152,10 @@ class TestMacNSpectrum:
     def test_peaks_at_true_mode(self):
         phi = np.array([0.0, 33.0, 66.0, 120.0, 200.0, 300.0])
         n_true = 3
-        z = shape_vector(-n_true * phi, np.ones_like(phi))  # phase = -n·φ
+        # phase = +n·φ: the conj(sig)·ref cross-phase ramp a cos(nφ − ωt) mode measures
+        z = shape_vector(n_true * phi, np.ones_like(phi))
         ns, macs, best = mac_n_spectrum(phi, z)
-        assert abs(best) == n_true
+        assert best == n_true  # signed
         assert macs.max() > 0.99
         assert ns.shape == macs.shape
 
@@ -259,7 +260,7 @@ def test_recovers_mode_number_during_mode(self):
 
         sigs, phi, t = self._array()
         tr = track_mode_shape(sigs, phi, t, frequency=8000.0, n_slices=20)
-        assert abs(int(np.median(tr.n_by_time[tr.t_ms < 25.0]))) == 2
+        assert int(np.median(tr.n_by_time[tr.t_ms < 25.0])) == 2  # signed
 
 
 class TestGPScaleInvariance:
diff --git a/tests/test_nodes.py b/tests/test_nodes.py
index ece9701..50677da 100644
--- a/tests/test_nodes.py
+++ b/tests/test_nodes.py
@@ -14,13 +14,18 @@
 
 
 def _first_shot():
-    # Prefer a DIII-D shot: these tests assume full geometry + QS arrays, which the
-    # synthetic KSTAR shot (device=kstar, toroidal-only, no r/z) does not yet carry.
-    ms = nodes.machines()
+    # The full DIII-D synthetic shot: these tests assume full geometry + QS arrays.
+    # (This used to take the first non-KSTAR machine, but list_shots() sorts shot ids
+    # as STRINGS, so the NSTX shot '204718' sorted before '990000' and the whole
+    # module silently ran against NSTX — making the DIII-D family-mixing guard
+    # tautological and permanently skipping the contour test.)
+    ms = {m["id"] for m in nodes.machines()}
     if not ms:
         pytest.skip("no fetched HDF5 in the data dir")
-    non_kstar = [m for m in ms if (m.get("device") or "").lower() != "kstar"]
-    return (non_kstar or ms)[0]["id"]
+    from .conftest import SYNTH_SHOT
+
+    assert str(SYNTH_SHOT) in ms, "synthetic DIII-D fixture shot missing from data dir"
+    return str(SYNTH_SHOT)
 
 
 def test_machines_shape():
@@ -88,10 +93,9 @@ def test_spectrogram_denoise_off_matches_default():
 
 def test_contour_node():
     shot = _first_shot()
-    try:
-        n = nodes.build_node(shot, "contour")
-    except Exception as e:  # noqa: BLE001 — shot may lack the MPID toroidal array
-        pytest.skip(f"no MPID toroidal array in this shot: {e}")
+    # The DIII-D fixture guarantees the MPID toroidal array — a builder crash here
+    # must be RED, not a skip (this test skipped for weeks while pointed at NSTX).
+    n = nodes.build_node(shot, "contour")
     assert n["kind"] == "contour"
     assert len(n["z"]) == len(n["y"]) and len(n["z"][0]) == len(n["x"])
 
@@ -133,10 +137,7 @@ def test_toroidal_array_single_family():
 
 def test_poloidal_shape_node():
     shot = _first_shot()
-    try:
-        n = nodes.build_node(shot, "poloidal_shape")
-    except Exception as e:  # noqa: BLE001 — shot may lack the MPID poloidal array
-        pytest.skip(f"no poloidal array in this shot: {e}")
+    n = nodes.build_node(shot, "poloidal_shape")
     assert n["kind"] == "line"
     assert {s["name"] for s in n["series"]} >= {"Re", "Im"}
     assert all("markers" in s for s in n["series"])
@@ -191,7 +192,10 @@ def test_fit_quality_node_has_finite_k():
     shot = _first_shot()
     n = nodes.build_node(shot, "fit_quality")
     assert n["kind"] == "metrics"
-    assert n["fields"]
+    k_fields = [f for f in n["fields"] if str(f["label"]).startswith("K")]
+    assert k_fields, n["fields"]
+    for f in k_fields:
+        assert float(f["value"]) == float(f["value"]) and float(f["value"]) < 1e18
 
 
 def test_sensor_map_rz_threads_shot_to_wall_loader(synthetic_shot, monkeypatch):
@@ -230,10 +234,7 @@ def test_real_theta_has_full_poloidal_coverage():
 
 def test_mode_pattern_node():
     shot = _first_shot()
-    try:
-        n = nodes.build_node(shot, "mode_pattern")
-    except Exception as e:  # noqa: BLE001 — shot may lack the poloidal array
-        pytest.skip(f"no poloidal array in this shot: {e}")
+    n = nodes.build_node(shot, "mode_pattern")
     assert n["kind"] == "contour"
     assert len(n["z"]) == len(n["y"]) and len(n["z"][0]) == len(n["x"])  # [θ][φ]
 
@@ -241,11 +242,6 @@ def test_mode_pattern_node():
 def test_elongation_theta_star_threads_into_poloidal_nodes(monkeypatch):
     """With κ available the poloidal axis is the corrected θ*; absent κ it's geometric."""
     shot = _first_shot()
-    try:
-        nodes.build_node(shot, "mode_pattern")  # needs the poloidal array
-    except Exception as e:  # noqa: BLE001
-        pytest.skip(f"no poloidal array in this shot: {e}")
-
     # κ absent → geometric θ, honest "no κ" label
     monkeypatch.setattr(nodes, "_kappa_at", lambda *a, **k: None)
     mp = nodes.build_node(shot, "mode_pattern", {"time": 3000})
@@ -277,10 +273,7 @@ def test_toroidal_stripes_node():
 
 def test_poloidal_phase_fit_node():
     shot = _first_shot()
-    try:
-        n = nodes.build_node(shot, "poloidal_phase_fit", {"time": 3000})
-    except Exception as e:  # noqa: BLE001 — shot may lack the poloidal array
-        pytest.skip(f"no poloidal array in this shot: {e}")
+    n = nodes.build_node(shot, "poloidal_phase_fit", {"time": 3000})
     assert n["kind"] == "scatter2d" and n["points"]
     assert "m_fit" in n["meta"]
 
@@ -313,3 +306,101 @@ def test_quality_for_k_thresholds():
     assert contracts.quality_for_k(15) == "warn"
     assert contracts.quality_for_k(25) == "bad"
     assert contracts.quality_for_k(float("nan")) == "bad"
+
+
+# ---------------------------------------------------------------------------
+# NaN samples in fetched channels must serialize as JSON null, not 500
+# ---------------------------------------------------------------------------
+
+
+def test_raw_trace_with_nan_samples_serializes(synthetic_shot, monkeypatch):
+    """starlette serializes with allow_nan=False, so a NaN gap in a raw channel
+    used to turn the whole node response into an opaque 500. Non-finite samples
+    must come out as JSON null (a Plotly gap)."""
+    import json
+
+    from magnetics.data import h5source
+
+    orig = h5source.load_channel
+
+    def nan_load(shot, name):
+        t, d = orig(shot, name)
+        d = np.asarray(d, dtype=float).copy()
+        d[::7] = np.nan
+        return t, d
+
+    monkeypatch.setattr(nodes.h5source, "load_channel", nan_load)
+    try:
+        node = nodes.build_node(synthetic_shot, "raw_trace")
+        json.dumps(node, allow_nan=False)  # must not raise
+        assert any(v is None for v in node["series"][0]["y"])
+        assert any(v is not None for v in node["series"][0]["y"])
+    finally:
+        nodes.refresh()  # don't leak NaN-poisoned caches into other tests
+
+
+def test_json_finite_helper():
+    from magnetics.core.contracts import json_finite
+
+    assert json_finite([1.0, np.nan, np.inf, -np.inf, 2.5]) == [1.0, None, None, None, 2.5]
+    assert json_finite([[1.0, np.nan], [3.0, 4.0]]) == [[1.0, None], [3.0, 4.0]]
+
+
+def test_quality_for_k_matches_ts_on_nonfinite():
+    """contract.ts uses !Number.isFinite(K) -> "bad"; the Python twin must agree
+    for NaN and BOTH infinities (K = -inf used to return "good")."""
+    assert contracts.quality_for_k(float("nan")) == "bad"
+    assert contracts.quality_for_k(float("inf")) == "bad"
+    assert contracts.quality_for_k(float("-inf")) == "bad"
+    assert contracts.quality_for_k(5.0) == "good"
+    assert contracts.quality_for_k(15.0) == "warn"
+    assert contracts.quality_for_k(25.0) == "bad"
+
+
+def test_refresh_clears_real_theta_cache(synthetic_shot):
+    """Regression: refresh() (fired after every fetch job) cleared every lru_cache
+    except _real_theta, so a re-pull with more probes kept serving the stale
+    channel->theta map and poloidal nodes 422'd until server restart."""
+    nodes._real_theta(str(synthetic_shot))
+    assert nodes._real_theta.cache_info().currsize > 0
+    nodes.refresh()
+    assert nodes._real_theta.cache_info().currsize == 0
+
+
+def test_channel_usage_tags_plasma_and_qs_channels(synthetic_shot):
+    """Regression: channel_usage marked kappa/ip/bt (and the QS fit array) as
+    "unused / droppable" although the theta* correction and the QS helicity
+    consume them — trimming a pull per that endpoint silently degraded physics."""
+    usage = nodes.channel_usage(synthetic_shot)
+    unused = set(usage["unused"])
+    present = {u["name"] for u in usage["used"]}
+    for nm in ("kappa", "ip", "bt"):
+        assert nm not in unused, f"{nm} reported droppable but analyses consume it"
+        assert nm in present
+    # the QS midplane array must carry a QS role, not sit in unused
+    roles = {u["name"]: u["roles"] for u in usage["used"]}
+    assert any("QS fit array" in r for rs in roles.values() for r in rs)
+
+
+def test_mode_number_band_follows_fmax(synthetic_shot, monkeypatch):
+    """Issue #85: the n-map's compute band was pinned at 0-50 kHz, so raising the
+    GUI's f_max never showed data above 50 kHz no matter the request. The requested
+    fmax must reach the array STFT as its band ceiling, quantized to 50-kHz steps
+    (so crops within a step still reuse the cached STFT)."""
+    captured = {}
+    orig = nodes.spectral.array_shape_spectrum
+
+    def spy(*a, **k):
+        captured["fmax"] = k.get("fmax")
+        return orig(*a, **k)
+
+    monkeypatch.setattr(nodes.spectral, "array_shape_spectrum", spy)
+    try:
+        nodes.refresh()
+        nodes.build_node(synthetic_shot, "mode_number", {"fmax": "120"})
+        assert captured["fmax"] == 150_000.0  # ceil(120/50) * 50 kHz
+        nodes.refresh()
+        nodes.build_node(synthetic_shot, "mode_number", {"fmax": "40"})
+        assert captured["fmax"] == 50_000.0  # within the default step
+    finally:
+        nodes.refresh()  # don't leak spy-built cache entries
diff --git a/tests/test_qs_bridge.py b/tests/test_qs_bridge.py
index f66fb19..d175934 100644
--- a/tests/test_qs_bridge.py
+++ b/tests/test_qs_bridge.py
@@ -58,8 +58,11 @@ def test_fit_quality_statuses_use_contract_vocabulary(fit_ds):
 
 def test_amplitude_sigma_is_finite(fit_ds):
     node = qs_bridge.fit_to_amplitude_node(fit_ds)
-    sigma = np.asarray(node["meta"]["sigma"], dtype=float)
-    assert np.all(np.isfinite(sigma))
+    for series in node["series"]:
+        lower = np.asarray(series["lower"], dtype=float)
+        upper = np.asarray(series["upper"], dtype=float)
+        assert np.all(np.isfinite(lower)) and np.all(np.isfinite(upper))
+        assert np.all(upper >= np.asarray(series["y"], dtype=float))
 
 
 def test_svd_energy_node_is_monotonic_and_bounded(fit_ds):
@@ -82,9 +85,14 @@ def test_svd_condition_node_is_finite_and_positive(fit_ds):
 def test_sigma_override_changes_amplitude_uncertainty(synthetic_shot):
     default_fit = nodes._prep_qs_ds(synthetic_shot, {}).fit
     overridden_fit = nodes._prep_qs_ds(synthetic_shot, {"sigma": "1.0"}).fit
-    default_sigma = qs_bridge.fit_to_amplitude_node(default_fit)["meta"]["sigma"]
-    overridden_sigma = qs_bridge.fit_to_amplitude_node(overridden_fit)["meta"]["sigma"]
-    assert np.mean(overridden_sigma) > np.mean(default_sigma)
+
+    def band_width(fit):
+        node = qs_bridge.fit_to_amplitude_node(fit)
+        return np.mean(
+            [np.mean(np.asarray(s["upper"]) - np.asarray(s["lower"])) for s in node["series"]]
+        )
+
+    assert band_width(overridden_fit) > band_width(default_fit)
 
 
 def test_fit_basis_param_reaches_fit(synthetic_shot):
@@ -130,3 +138,45 @@ def test_reconstruction_uses_minus_i_sign_convention():
         f"reconstructed peak at φ={peak_phi}°, expected ~{delta}° (−i convention); "
         f"a peak near {360 - delta}° means the exp(+i…) sign bug is back"
     )
+
+
+def test_sigma_override_does_not_corrupt_fit_signal(synthetic_shot):
+    """Regression: fit_signal was de-normalized by the dataset signal_sigma while the
+    design matrix and RHS were normalized by the override sigma, scaling fit_signal
+    by signal_sigma/override and corrupting residual and chi_sq (the coefficients
+    stayed correct). With a uniform sigma the weighted LS solution is
+    override-invariant: fit_signal must be identical, and chi_sq must scale by
+    exactly (sigma_default / sigma_override)**2."""
+    default = nodes._prep_qs_ds(synthetic_shot, {}).fit
+    override = 1e-4
+    overridden = nodes._prep_qs_ds(synthetic_shot, {"sigma": str(override)}).fit
+
+    np.testing.assert_allclose(
+        overridden["fit_coeffs"].values, default["fit_coeffs"].values, rtol=1e-9
+    )
+    np.testing.assert_allclose(
+        overridden["fit_signal"].values, default["fit_signal"].values, rtol=1e-9
+    )
+    sig0 = float(np.nanmean(default["signal_sigma"].values))
+    np.testing.assert_allclose(
+        overridden["chi_sq"].values,
+        default["chi_sq"].values * (sig0 / override) ** 2,
+        rtol=1e-6,
+    )
+
+
+def test_default_fit_cond_is_inversion_cutoff_not_trust_threshold(synthetic_shot):
+    """The default cutoff must be OMFIT SLCONTOUR's 1e3 inversion cutoff (1/rcond)
+    at every layer. When it was 10 (the GUI's K-trust threshold), any fit with K in
+    (10, 1e3) silently zeroed basis directions the reference fit keeps — and K(eff)
+    read <= 10 by construction, so the quality panel looked healthiest exactly when
+    the regularization was distorting the fit."""
+    import inspect
+
+    from magnetics.core import qs_fit
+
+    # core default (what a direct fit() call regularizes with)
+    assert inspect.signature(qs_fit.fit).parameters["fit_cond"].default == 1e3
+    # service default (what a GUI request without fit_cond regularizes with)
+    run = nodes._prep_qs_ds(synthetic_shot, {})
+    assert run.fit.attrs["fit_condition"] == 1e3
diff --git a/tests/test_qs_prep_branches.py b/tests/test_qs_prep_branches.py
new file mode 100644
index 0000000..db833a0
--- /dev/null
+++ b/tests/test_qs_prep_branches.py
@@ -0,0 +1,120 @@
+"""qs_prep branches that the end-to-end pipeline tests never reach.
+
+The GUI offers four detrend types, but only `baseline` ever ran under test; the
+`linear`/`endpoints` estimators and the DIII-D 2019 ESLD wiring-swap branch were
+complete dead zones. All offline: `_detrend` is exercised directly on hand-built
+Datasets, the swap through `prepare()` on a minimal ShotData.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+import xarray as xr
+
+from magnetics.core.qs_io_data import ShotData
+from magnetics.core.qs_prep import _detrend, prepare
+
+TIME = np.linspace(0.0, 1.0, 201)  # s
+
+
+def _ds(*signals):
+    names = [f"c{i}" for i in range(len(signals))]
+    return xr.Dataset(
+        {"signal": (("channel", "time"), np.array(signals, dtype=float))},
+        coords={"channel": names, "time": TIME},
+    )
+
+
+def _quiet(*a):
+    pass
+
+
+class TestDetrend:
+    def test_none_is_a_no_op(self):
+        ds = _ds(np.sin(2 * np.pi * 5 * TIME) + 3.0)
+        before = ds["signal"].values.copy()
+        _detrend(ds, ds["channel"].values, "none", (0.0, 0.2), _quiet)
+        np.testing.assert_array_equal(ds["signal"].values, before)
+
+    def test_baseline_removes_in_band_mean(self):
+        ds = _ds(np.sin(2 * np.pi * 50 * TIME) + 7.5)
+        _detrend(ds, ds["channel"].values, "baseline", (0.0, 0.2), _quiet)
+        band = (TIME >= 0.0) & (TIME <= 0.2)
+        assert np.mean(ds["signal"].values[0, band]) == pytest.approx(0.0, abs=1e-12)
+
+    def test_linear_removes_a_pure_ramp_everywhere(self):
+        # The trend is fit inside the band but subtracted over the WHOLE record —
+        # a pure ramp must vanish everywhere, not just in the band.
+        ds = _ds(4.0 * TIME + 2.0)
+        _detrend(ds, ds["channel"].values, "linear", (0.1, 0.4), _quiet)
+        np.testing.assert_allclose(ds["signal"].values[0], 0.0, atol=1e-10)
+
+    def test_linear_preserves_oscillation_on_top_of_ramp(self):
+        osc = np.sin(2 * np.pi * 40 * TIME)
+        ds = _ds(osc + 4.0 * TIME + 2.0)
+        _detrend(ds, ds["channel"].values, "linear", (0.0, 1.0), _quiet)
+        # residual is the oscillation (the LSQ line through a whole number of
+        # periods absorbs none of it to good approximation)
+        assert np.corrcoef(ds["signal"].values[0], osc)[0, 1] > 0.999
+
+    def test_endpoints_zeroes_the_band_endpoints(self):
+        rng = np.random.default_rng(7)
+        ds = _ds(rng.standard_normal(TIME.size) + 3.0 * TIME)
+        band = (0.1, 0.9)
+        sel = (TIME >= band[0]) & (TIME <= band[1])
+        first, last = np.flatnonzero(sel)[[0, -1]]
+        _detrend(ds, ds["channel"].values, "endpoints", band, _quiet)
+        # the line through the two band endpoints is removed → both go to ~0
+        assert ds["signal"].values[0, first] == pytest.approx(0.0, abs=1e-10)
+        assert ds["signal"].values[0, last] == pytest.approx(0.0, abs=1e-10)
+
+    def test_per_channel_independence(self):
+        ds = _ds(np.full(TIME.size, 5.0), np.full(TIME.size, -2.0))
+        _detrend(ds, ds["channel"].values, "baseline", (0.0, 1.0), _quiet)
+        np.testing.assert_allclose(ds["signal"].values, 0.0, atol=1e-12)
+
+    def test_unknown_type_raises(self):
+        ds = _ds(TIME)
+        with pytest.raises(ValueError, match="detrend_type"):
+            _detrend(ds, ds["channel"].values, "quadratic", (0.0, 1.0), _quiet)
+
+
+# ---------------------------------------------------------------------------
+# DIII-D 2019 wiring mix-up: ESLD66M079 ⇄ ESLD66M319 swap for shot > 177705
+# ---------------------------------------------------------------------------
+
+
+def _shotdata(shot, ch_a="ESLD66M079", ch_b="ESLD66M319"):
+    t = np.linspace(2.85, 3.05, 2001)  # spans prepare()'s default (2.9, 3.0) trim
+    raw = xr.Dataset(
+        {
+            "signal": (
+                ("channel", "time"),
+                np.vstack([np.full(t.size, 1.0), np.full(t.size, 2.0)]),
+            ),
+            "signal_sigma": (("channel",), np.array([2e-5, 2e-5])),
+        },
+        coords={"channel": [ch_a, ch_b], "time": t},
+    )
+    plasma = xr.Dataset(coords={"time": t})
+    return ShotData(shot=shot, device="DIII-D", raw=raw, plasma=plasma)
+
+
+def test_esld_2019_swap_applies_after_the_wiring_fix_shot():
+    sd = _shotdata(180000)  # > 177705
+    prepared, _ = prepare(sd, channel_filter=".*", detrend_type="none", verbose=False)
+    # channel labeled 079 must carry 319's samples and vice versa
+    assert float(prepared["signal"].sel(channel="ESLD66M079").values[0]) == pytest.approx(2.0)
+    assert float(prepared["signal"].sel(channel="ESLD66M319").values[0]) == pytest.approx(1.0)
+
+
+def test_esld_swap_skipped_for_pre_fix_shots_and_partial_pairs():
+    sd = _shotdata(170000)  # < 177705: wiring was correct, no swap
+    prepared, _ = prepare(sd, channel_filter=".*", detrend_type="none", verbose=False)
+    assert float(prepared["signal"].sel(channel="ESLD66M079").values[0]) == pytest.approx(1.0)
+
+    # only one of the pair present → swap must not run (guard on both names)
+    sd = _shotdata(180000, ch_b="MPID66M067")
+    prepared, _ = prepare(sd, channel_filter=".*", detrend_type="none", verbose=False)
+    assert float(prepared["signal"].sel(channel="ESLD66M079").values[0]) == pytest.approx(1.0)
diff --git a/tests/test_remote_orchestration.py b/tests/test_remote_orchestration.py
new file mode 100644
index 0000000..d7ad09f
--- /dev/null
+++ b/tests/test_remote_orchestration.py
@@ -0,0 +1,151 @@
+"""run_remote's SSH/rsync command construction, with subprocess faked.
+
+The cluster round-trip can't run offline, but every command it builds is pure
+logic — and it has already burned us once (the remote result clobbering the
+local shot file). These tests script `subprocess.run` and assert the exact
+orchestration: master connect + jump selection, per-user /tmp staging, the
+remote fetch argv, the rsync copy-back target, cleanup, and teardown.
+"""
+
+from __future__ import annotations
+
+import subprocess
+from types import SimpleNamespace
+
+import pytest
+
+from magnetics.data.fetch import remote
+
+
+class _Runner:
+    """Record every subprocess.run; script per-command return codes/stdout."""
+
+    def __init__(self, fail_on=None):
+        self.calls: list[list[str]] = []
+        self.fail_on = fail_on  # substring of the joined argv that should fail
+
+    def __call__(self, cmd, env=None, **kw):
+        self.calls.append(list(cmd))
+        joined = " ".join(str(c) for c in cmd)
+        rc = 1 if (self.fail_on and self.fail_on in joined) else 0
+        stdout = "gauser\n" if "id -un" in joined else ""
+        return SimpleNamespace(returncode=rc, stdout=stdout, stderr="")
+
+    def joined(self):
+        return [" ".join(str(c) for c in call) for call in self.calls]
+
+
+@pytest.fixture()
+def runner(monkeypatch):
+    r = _Runner()
+    monkeypatch.setattr(subprocess, "run", r)
+    monkeypatch.setattr(remote, "on_site_network", lambda device: False)  # laptop
+    return r
+
+
+def test_off_site_pull_builds_the_full_round_trip(runner, tmp_path):
+    out = remote.run_remote(
+        184927,
+        "both",
+        username="myuser",
+        tmin=1000,
+        tmax=2000,
+        decimate=1,
+        device="diiid",
+        local_out_dir=str(tmp_path),
+    )
+    cmds = runner.joined()
+
+    # 1) one authenticated master with the site gateway as ProxyJump
+    master = cmds[0]
+    assert "ControlMaster=auto" in master and "myuser@omega.gat.com" in master
+    assert "-J myuser@cybele.gat.com:2039" in master
+
+    # 2) per-user /tmp stage (shared login node) + code sync via the master
+    assert any("mkdir -p" in c and "/tmp/magnetics_out_myuser" in c for c in cmds)
+    rsync_up = next(c for c in cmds if c.startswith("rsync") and "magnetics" in c)
+    assert "myuser@omega.gat.com:~/magnetics_fetch/" in rsync_up
+
+    # 3) the remote fetch runs THIS package with the toksearch backend and the
+    #    caller's window, staging into the per-user /tmp file
+    fetch = next(c for c in cmds if "-m magnetics.data.fetch.toksearch" in c)
+    for frag in (
+        "--backend toksearch",
+        "--shot 184927",
+        "--analysis both",
+        "--tmin 1000",
+        "--tmax 2000",
+        "--out /tmp/magnetics_out_myuser/shot_184927.h5",
+    ):
+        assert frag in fetch, f"{frag!r} missing from: {fetch}"
+    assert "--decimate" not in fetch  # decimate=1 is the default, don't forward it
+
+    # 4) the result is rsync'd back to the REQUESTED local dir, the stage removed,
+    #    and the ControlMaster torn down
+    rsync_back = next(c for c in cmds if c.startswith("rsync") and str(tmp_path) in c)
+    assert "myuser@omega.gat.com:/tmp/magnetics_out_myuser/shot_184927.h5" in rsync_back
+    assert any("rm -f" in c and "shot_184927.h5" in c for c in cmds)
+    assert "-O exit" in cmds[-1]
+    assert out == str(tmp_path / "shot_184927.h5")
+
+
+def test_custom_pointnames_replace_analysis_selection(runner, tmp_path):
+    remote.run_remote(
+        1, "both", username="u", raw_pointnames=["betan", "li"], local_out_dir=str(tmp_path)
+    )
+    fetch = next(c for c in runner.joined() if "-m magnetics.data.fetch.toksearch" in c)
+    assert "--pointnames betan,li" in fetch
+    assert "--analysis" not in fetch and "--sensor-set" not in fetch
+
+
+def test_sensor_set_overrides_analysis(runner, tmp_path):
+    remote.run_remote(1, "rotating", username="u", sensor_set="Bp All", local_out_dir=str(tmp_path))
+    fetch = next(c for c in runner.joined() if "-m magnetics.data.fetch.toksearch" in c)
+    assert (
+        "--sensor-set 'Bp All'" in fetch
+        or '--sensor-set "Bp All"' in fetch
+        or "--sensor-set Bp All" in fetch
+    )
+    assert "--analysis" not in fetch
+
+
+def test_on_site_drops_the_jump(monkeypatch, tmp_path):
+    r = _Runner()
+    monkeypatch.setattr(subprocess, "run", r)
+    monkeypatch.setattr(remote, "on_site_network", lambda device: True)
+    remote.run_remote(1, "both", username="u", local_out_dir=str(tmp_path))
+    assert "-J" not in r.joined()[0]
+
+
+def test_explicit_empty_jump_uses_alias_proxyjump(runner, tmp_path):
+    # jump="" means "my ssh-config alias carries its own ProxyJump" — forcing -J
+    # on top would double-jump.
+    remote.run_remote(1, "both", host="myalias", jump="", local_out_dir=str(tmp_path))
+    assert "-J" not in runner.joined()[0]
+
+
+def test_unknown_user_is_resolved_on_the_cluster(runner, tmp_path):
+    # no username: an ssh-config alias supplies it, so the stage-dir owner is
+    # asked from the cluster (`id -un`, works with $USER unset in command shells)
+    remote.run_remote(1, "both", local_out_dir=str(tmp_path))
+    cmds = runner.joined()
+    assert any("id -un" in c for c in cmds)
+    assert any("/tmp/magnetics_out_gauser" in c for c in cmds)
+
+
+def test_failed_master_connect_exits_before_any_sync(monkeypatch, tmp_path):
+    r = _Runner(fail_on="true")  # the master probe command ends in `true`
+    monkeypatch.setattr(subprocess, "run", r)
+    monkeypatch.setattr(remote, "on_site_network", lambda device: False)
+    with pytest.raises(SystemExit):
+        remote.run_remote(1, "both", username="u", local_out_dir=str(tmp_path))
+    # nothing was synced or executed remotely after the failed connect
+    assert not any("rsync" in c for c in r.joined())
+
+
+def test_unsafe_remote_dir_rejected_before_any_ssh(monkeypatch, tmp_path):
+    r = _Runner()
+    monkeypatch.setattr(subprocess, "run", r)
+    with pytest.raises(ValueError, match="unsafe remote_dir"):
+        remote.run_remote(1, "both", remote_dir="~/x; rm -rf /", local_out_dir=str(tmp_path))
+    assert r.calls == []
diff --git a/tests/test_slcontour_fit.py b/tests/test_slcontour_fit.py
index 0f404d7..705d0e2 100644
--- a/tests/test_slcontour_fit.py
+++ b/tests/test_slcontour_fit.py
@@ -87,3 +87,178 @@ def test_fit_reports_covariance_diagonal_sigmas():
 def test_bad_basis_raises():
     with pytest.raises(ValueError):
         form_basis_function(1, 0, 0.0, 1.0, 0.0, 1.0, fit_basis="not-a-basis")
+
+
+# ---------------------------------------------------------------------------
+# Gaussian bases — value-pinned against direct numerical integration.
+# These branches had NO tests, and this function is exactly where the port
+# deliberately diverges from the OMFIT source (its `)(` typo and integral
+# scale bugs) — so the intended semantics must be pinned, not inherited.
+# ---------------------------------------------------------------------------
+
+
+class TestGaussianBases:
+    X1, X2 = np.array([10.0]), np.array([30.0])  # sensor phi extent (deg)
+    Y1, Y2 = np.array([-5.0]), np.array([15.0])  # sensor theta extent (deg)
+    NE, ME = 25.0, 30.0  # RBF widths (deg)
+    CX, CY = 40.0, 10.0  # RBF centre (deg) — n/m carry the CENTRE for gaussians
+
+    def test_gaussian_point_is_rbf_at_sensor_centre(self):
+        v = form_basis_function(
+            self.CX,
+            self.CY,
+            self.X1,
+            self.X2,
+            self.Y1,
+            self.Y2,
+            "gaussian-point",
+            nepsilon=self.NE,
+            mepsilon=self.ME,
+        )
+        xc, yc = 20.0, 5.0  # midpoints of the extents
+        expect = np.exp(-(((self.CX - xc) / self.NE) ** 2 + ((self.CY - yc) / self.ME) ** 2))
+        assert v.dtype.kind == "f"
+        assert v[0] == pytest.approx(expect, rel=1e-12)
+
+    def test_gaussian_point_uniform_direction_drops_term(self):
+        # mepsilon=inf → ridge in phi only (theta term contributes nothing)
+        v = form_basis_function(
+            self.CX,
+            999.0,
+            self.X1,
+            self.X2,
+            self.Y1,
+            self.Y2,
+            "gaussian-point",
+            nepsilon=self.NE,
+            mepsilon=np.inf,
+        )
+        assert v[0] == pytest.approx(np.exp(-(((self.CX - 20.0) / self.NE) ** 2)), rel=1e-12)
+
+    def test_gaussian_point_periodic_copies_wrap(self):
+        # centre at 350°, sensor centred at 20°: the +360° copy (cx=710? no — cx-360=-10)
+        # sits 30° away and dominates the direct 330° separation.
+        base = form_basis_function(
+            350.0,
+            self.CY,
+            self.X1,
+            self.X2,
+            self.Y1,
+            self.Y2,
+            "gaussian-point",
+            ncycle=0,
+            nepsilon=self.NE,
+            mepsilon=self.ME,
+        )
+        wrapped = form_basis_function(
+            350.0,
+            self.CY,
+            self.X1,
+            self.X2,
+            self.Y1,
+            self.Y2,
+            "gaussian-point",
+            ncycle=1,
+            nepsilon=self.NE,
+            mepsilon=self.ME,
+        )
+        manual = sum(
+            np.exp(
+                -((((350.0 + k * 360) - 20.0) / self.NE) ** 2 + ((self.CY - 5.0) / self.ME) ** 2)
+            )
+            for k in (-1, 0, 1)
+        )
+        assert wrapped[0] == pytest.approx(manual, rel=1e-12)
+        assert wrapped[0] > base[0] * 10  # the wrap matters physically
+
+    def test_gaussian_integral_2d_matches_dblquad(self):
+        from scipy.integrate import dblquad
+
+        v = form_basis_function(
+            self.CX,
+            self.CY,
+            self.X1,
+            self.X2,
+            self.Y1,
+            self.Y2,
+            "gaussian-integral",
+            nepsilon=self.NE,
+            mepsilon=self.ME,
+        )
+        num, _ = dblquad(
+            lambda y, x: np.exp(-(((x - self.CX) / self.NE) ** 2 + ((y - self.CY) / self.ME) ** 2)),
+            float(self.X1[0]),
+            float(self.X2[0]),
+            float(self.Y1[0]),
+            float(self.Y2[0]),
+        )
+        # the basis IS the plain area integral of the RBF over the sensor extent
+        assert v[0] == pytest.approx(num, rel=1e-10)
+
+    @pytest.mark.parametrize("uniform", ["theta", "phi"])
+    def test_gaussian_integral_ridge_matches_quad(self, uniform):
+        from scipy.integrate import quad
+
+        if uniform == "theta":  # mepsilon=inf → 1-D line integral along phi
+            v = form_basis_function(
+                self.CX,
+                0.0,
+                self.X1,
+                self.X2,
+                self.Y1,
+                self.Y2,
+                "gaussian-integral",
+                nepsilon=self.NE,
+                mepsilon=np.inf,
+            )
+            num, _ = quad(
+                lambda x: np.exp(-(((x - self.CX) / self.NE) ** 2)),
+                float(self.X1[0]),
+                float(self.X2[0]),
+            )
+        else:  # nepsilon=inf → 1-D line integral along theta
+            v = form_basis_function(
+                0.0,
+                self.CY,
+                self.X1,
+                self.X2,
+                self.Y1,
+                self.Y2,
+                "gaussian-integral",
+                nepsilon=np.inf,
+                mepsilon=self.ME,
+            )
+            num, _ = quad(
+                lambda y: np.exp(-(((y - self.CY) / self.ME) ** 2)),
+                float(self.Y1[0]),
+                float(self.Y2[0]),
+            )
+        assert v[0] == pytest.approx(num, rel=1e-10)
+
+    def test_sinusoidal_integral_mixed_mode_is_exact_area_average(self):
+        """The (n≠0, m≠0) sinusoidal-integral value must be the exact area-average
+        of e^{i(nφ+mθ)} over the sensor extent. This is the site where the port
+        DELIBERATELY diverges from OMFIT (whose expression is off by −180/π); the
+        divergence is intentional and this pins the correct semantics."""
+        from scipy.integrate import dblquad
+
+        n_mode, m_mode = 2, 1
+        v = form_basis_function(
+            n_mode, m_mode, self.X1, self.X2, self.Y1, self.Y2, "sinusoidal-integral"
+        )
+        area = np.deg2rad(self.X2[0] - self.X1[0]) * np.deg2rad(self.Y2[0] - self.Y1[0])
+
+        def integrand(re):
+            def f(y, x):
+                val = np.exp(1j * (n_mode * np.deg2rad(x) + m_mode * np.deg2rad(y)))
+                return val.real if re else val.imag
+
+            num, _ = dblquad(
+                f, float(self.X1[0]), float(self.X2[0]), float(self.Y1[0]), float(self.Y2[0])
+            )
+            return num
+
+        # ∫∫ over degrees × (rad/deg)² = ∫∫ over radians; divide by the rad² area
+        avg = (integrand(True) + 1j * integrand(False)) * np.deg2rad(1.0) ** 2 / area
+        assert v[0].real == pytest.approx(avg.real, rel=1e-9)
+        assert v[0].imag == pytest.approx(avg.imag, rel=1e-9)
diff --git a/tests/test_spectral.py b/tests/test_spectral.py
index a2b810b..7b1535b 100644
--- a/tests/test_spectral.py
+++ b/tests/test_spectral.py
@@ -101,7 +101,7 @@ def test_recovers_synthetic_mode_number(self, synthetic_n2):
         d = synthetic_n2
         result = cross_spectrum(d["sig1"], d["sig2"], d["fs"], delta_phi=d["delta_phi"])
         peak_idx = np.argmax(result.power)
-        assert abs(result.mode_number[peak_idx]) == d["n_true"]
+        assert result.mode_number[peak_idx] == d["n_true"]  # signed
 
     def test_peak_frequency_near_expected(self, synthetic_n2):
         d = synthetic_n2
@@ -113,7 +113,7 @@ def test_rms_by_mode_peaks_at_true_n(self, synthetic_n2):
         d = synthetic_n2
         result = cross_spectrum(d["sig1"], d["sig2"], d["fs"], delta_phi=d["delta_phi"])
         peak_mode_idx = np.argmax(result.rms_by_mode)
-        assert abs(result.mode_indices[peak_mode_idx]) == d["n_true"]
+        assert result.mode_indices[peak_mode_idx] == d["n_true"]  # signed
 
     def test_coherence_high_for_correlated_signals(self, synthetic_n2):
         d = synthetic_n2
@@ -168,11 +168,11 @@ def test_recovers_mode_across_all_slices(self, synthetic_n2):
         )
         for i in range(len(result.time)):
             peak_idx = np.argmax(result.power[i])
-            assert abs(result.mode_number[i, peak_idx]) == d["n_true"]
+            assert result.mode_number[i, peak_idx] == d["n_true"]  # signed
 
     def test_on_real_data(self, shot_174446):
         d = shot_174446
-        delta_phi = d["phi_307"] - d["phi_340"]
+        delta_phi = d["phi_340"] - d["phi_307"]  # φ2 − φ1, service convention
         result = compute_spectrogram(
             d["time_s"],
             d["sig_307"],
@@ -290,7 +290,7 @@ def test_preserves_dominant_mode_peak(self, synthetic_n2):
         dn = denoise_spectrogram(spec, coherence_min=0.5, power_floor_k=None)
         peak = np.unravel_index(np.argmax(spec.power), spec.power.shape)
         assert dn.power[peak] > 0
-        assert abs(dn.mode_number[peak]) == d["n_true"]
+        assert dn.mode_number[peak] == d["n_true"]  # signed
 
     def test_power_floor_removes_more_than_coherence_alone(self, synthetic_n2):
         spec = self._spec(synthetic_n2)
@@ -300,7 +300,7 @@ def test_power_floor_removes_more_than_coherence_alone(self, synthetic_n2):
 
     def test_on_real_data(self, shot_174446):
         d = shot_174446
-        delta_phi = d["phi_307"] - d["phi_340"]
+        delta_phi = d["phi_340"] - d["phi_307"]  # φ2 − φ1, service convention
         spec = compute_spectrogram(
             d["time_s"], d["sig_307"], d["sig_340"], delta_phi, slice_duration=0.004
         )
@@ -513,7 +513,7 @@ def test_reference_probe_flagged_nan(self):
     def test_fit_reports_confidence_and_sigma(self):
         mode, n_true = self._modes()
         fit = fit_toroidal_mode(mode)
-        assert abs(fit.n) == n_true
+        assert fit.n == n_true  # signed
         assert fit.phase_sigma is not None and fit.phase_sigma > 0
         assert fit.n_confidence is not None
         assert 0.0 < fit.n_confidence <= 1.0
@@ -656,7 +656,7 @@ def test_recovers_high_mode_number_at_peak(self, n_true):
         fi = int(np.argmin(np.abs(ms.freq_band - 8000.0)))  # the injected mode bin
         it = ms.time.size // 2
         n_fit = fit_toroidal_mode(mode_from_spectrum(spec, phi, t[it], 8000.0)).n
-        assert abs(ms.mode_number[it, fi]) == n_true  # right |n| (not aliased)
+        assert ms.mode_number[it, fi] == n_true  # right SIGNED n (not aliased)
         assert ms.mode_number[it, fi] == n_fit  # agrees with the phase fit
         assert ms.quality[it, fi] > 0.9  # clean single-n → high q
 
@@ -688,3 +688,168 @@ def test_smooth_mode_spectrogram(self):
         assert sm.quality.shape == ms.quality.shape
         assert np.all((sm.quality >= 0.0) & (sm.quality <= 1.0 + 1e-6))
         np.testing.assert_array_equal(sm.mode_number, ms.mode_number)
+
+
+# -----------------------------------------------------------------------
+# wide-pair aliasing regression (probes 320° apart the long way = -40° wrapped)
+# -----------------------------------------------------------------------
+
+
+class TestWidePairAliasing:
+    """Regression for the production DIII-D pair MPI66M020D/340D (φ = 19.52°/339.7°):
+    Δφ = φ2 − φ1 = 320.2° the long way, −39.8° wrapped. Unwrapped, round(phase/Δφ)
+    reported n = 0 for every real mode. n is SIGNED here — sig ~ sin(ωt − nφ) with
+    Δφ = φ2 − φ1 (the service convention) must recover +n."""
+
+    PHI1, PHI2 = 19.52, 339.7  # real MPI66M020D / MPI66M340D toroidal angles
+
+    @staticmethod
+    def _pair(n, phi1, phi2, fs=50_000, f_mode=3_000.0, duration=0.1):
+        t = np.linspace(0, duration, int(fs * duration), endpoint=False)
+        w = 2 * np.pi * f_mode
+        sig1 = np.sin(w * t - np.deg2rad(n * phi1))
+        sig2 = np.sin(w * t - np.deg2rad(n * phi2))
+        return t, sig1, sig2, float(fs)
+
+    @pytest.mark.parametrize("n_true", [1, 2, 3, 4, -2])
+    def test_cross_spectrum_recovers_signed_n(self, n_true):
+        _, sig1, sig2, fs = self._pair(n_true, self.PHI1, self.PHI2)
+        result = cross_spectrum(sig1, sig2, fs, delta_phi=self.PHI2 - self.PHI1)
+        peak = np.argmax(result.power)
+        assert result.mode_number[peak] == n_true
+
+    @pytest.mark.parametrize("n_true", [1, 3])
+    def test_spectrogram_recovers_signed_n(self, n_true):
+        t, sig1, sig2, fs = self._pair(n_true, self.PHI1, self.PHI2)
+        res = compute_spectrogram(t, sig1, sig2, delta_phi=self.PHI2 - self.PHI1)
+        it, fi = np.unravel_index(np.argmax(res.power), res.power.shape)
+        assert res.mode_number[it, fi] == n_true
+
+    def test_long_and_short_way_separations_agree(self):
+        _, sig1, sig2, fs = self._pair(2, self.PHI1, self.PHI2)
+        long_way = cross_spectrum(sig1, sig2, fs, delta_phi=self.PHI2 - self.PHI1)
+        short_way = cross_spectrum(sig1, sig2, fs, delta_phi=self.PHI2 - self.PHI1 - 360.0)
+        np.testing.assert_array_equal(long_way.mode_number, short_way.mode_number)
+        np.testing.assert_array_equal(long_way.mode_indices, short_way.mode_indices)
+
+    def test_full_turn_separation_raises(self):
+        _, sig1, sig2, fs = self._pair(1, 0.0, 360.0)
+        with pytest.raises(ValueError, match="mod 360"):
+            cross_spectrum(sig1, sig2, fs, delta_phi=360.0)
+
+
+class TestWrapAngleDeg:
+    def test_values(self):
+        from magnetics.core.spectral import wrap_angle_deg
+
+        assert wrap_angle_deg(320.18) == pytest.approx(-39.82)
+        assert wrap_angle_deg(-200.0) == pytest.approx(160.0)
+        assert wrap_angle_deg(33.0) == pytest.approx(33.0)
+        assert wrap_angle_deg(180.0) == 180.0
+        assert wrap_angle_deg(-180.0) == 180.0
+        assert wrap_angle_deg(540.0) == 180.0
+        assert wrap_angle_deg(720.0) == 0.0
+
+
+# -----------------------------------------------------------------------
+# cross-family signed-n agreement — every estimator, one physical mode
+# -----------------------------------------------------------------------
+
+
+class TestSignedNConventionAgreement:
+    """One physical mode b ∝ cos(nφ − ωt) (n = +2, rotating toward +φ) must be
+    reported with the SAME signed n by every estimator family: the 2-point
+    cross-spectrum/spectrogram, the phase-ramp fit, the array projection, and the
+    shape MAC. Regression for the families disagreeing in sign (docstrings claimed
+    agreement while the fit family returned −n)."""
+
+    N_TRUE = 2
+    F_MODE = 8_000.0
+
+    @pytest.fixture()
+    def array_mode(self):
+        fs = 200_000
+        t = np.linspace(0, 0.05, int(fs * 0.05), endpoint=False)
+        phi = np.linspace(0.0, 330.0, 12)
+        rng = np.random.default_rng(11)
+        sigs = np.vstack(
+            [
+                np.sin(2 * np.pi * self.F_MODE * t - np.deg2rad(self.N_TRUE * p))
+                + 0.02 * rng.standard_normal(t.size)
+                for p in phi
+            ]
+        )
+        return sigs, phi, t, float(fs)
+
+    def test_all_families_agree_on_signed_n(self, array_mode):
+        from magnetics.core.mode_shape import mac_n_spectrum, shape_vector
+        from magnetics.core.spectral import (
+            array_mode_spectrogram,
+            array_shape_spectrum,
+        )
+
+        sigs, phi, t, fs = array_mode
+
+        # family A1: 2-point cross-spectrum (adjacent probes, Δφ = φ2 − φ1 = 30°)
+        two_pt = cross_spectrum(sigs[0], sigs[1], fs, delta_phi=phi[1] - phi[0])
+        n_2pt = int(two_pt.mode_number[np.argmax(two_pt.power)])
+
+        # family A2: 2-point spectrogram
+        spec2 = compute_spectrogram(t, sigs[0], sigs[1], delta_phi=phi[1] - phi[0])
+        it, fi = np.unravel_index(np.argmax(spec2.power), spec2.power.shape)
+        n_specgram = int(spec2.mode_number[it, fi])
+
+        # family B1: phase-ramp fit on the full array
+        mode = extract_mode_at_frequency(sigs, phi, t, frequency=self.F_MODE)
+        n_fit = fit_toroidal_mode(mode).n
+
+        # family B2: array harmonic projection
+        aspec = array_shape_spectrum(sigs, t)
+        ms = array_mode_spectrogram(aspec, phi)
+        fb = int(np.argmin(np.abs(ms.freq_band - self.F_MODE)))
+        n_array = int(ms.mode_number[ms.time.size // 2, fb])
+
+        # family B3: shape MAC
+        _, _, n_mac = mac_n_spectrum(phi, shape_vector(mode.phase, mode.amplitude))
+
+        assert n_2pt == n_specgram == n_fit == n_array == n_mac == self.N_TRUE, (
+            f"2pt={n_2pt} specgram={n_specgram} fit={n_fit} array={n_array} mac={n_mac}"
+        )
+
+
+# -----------------------------------------------------------------------
+# spectrogram power normalization — must match csd's density scaling
+# -----------------------------------------------------------------------
+
+
+class TestSpectrogramPowerNormalization:
+    """Regression: compute_spectrogram used raw |conj(S2)·S1| with no window-power
+    or fs scaling yet applied cross_spectrum's density-RMS formula, so its
+    'rms amplitude' was off by a data-dependent factor of thousands."""
+
+    def _pair(self, amp=1.0, fs=50_000, f_mode=3_000.0, duration=0.2):
+        t = np.linspace(0, duration, int(fs * duration), endpoint=False)
+        w = 2 * np.pi * f_mode
+        n, phi1, phi2 = 2, 30.0, 63.0
+        sig1 = amp * np.sin(w * t - np.deg2rad(n * phi1))
+        sig2 = amp * np.sin(w * t - np.deg2rad(n * phi2))
+        return t, sig1, sig2, float(fs), phi2 - phi1
+
+    def test_rms_by_mode_recovers_physical_rms(self):
+        # a unit-amplitude sine has RMS 1/sqrt(2); both estimators must report it
+        t, sig1, sig2, fs, dphi = self._pair(amp=1.0)
+        two_pt = cross_spectrum(sig1, sig2, fs, delta_phi=dphi, nperseg=1024)
+        spec = compute_spectrogram(t, sig1, sig2, dphi, slice_duration=0.02)
+        rms_2pt = float(np.max(two_pt.rms_by_mode))
+        rms_spec = float(np.max(spec.rms_by_mode))
+        assert rms_2pt == pytest.approx(2**-0.5, rel=0.05)
+        assert rms_spec == pytest.approx(2**-0.5, rel=0.05)
+
+    def test_rms_scales_linearly_with_amplitude(self):
+        _, s1a, s2a, fs, dphi = self._pair(amp=1.0)
+        t, s1b, s2b, _, _ = self._pair(amp=3.0)
+        a = compute_spectrogram(t, s1a, s2a, dphi, slice_duration=0.02)
+        b = compute_spectrogram(t, s1b, s2b, dphi, slice_duration=0.02)
+        assert float(np.max(b.rms_by_mode)) == pytest.approx(
+            3.0 * float(np.max(a.rms_by_mode)), rel=1e-3
+        )
diff --git a/tests/test_toksearch_writer.py b/tests/test_toksearch_writer.py
index 461948b..668a4c9 100644
--- a/tests/test_toksearch_writer.py
+++ b/tests/test_toksearch_writer.py
@@ -182,3 +182,153 @@ def compute_serial(self):
     assert by_name["PT"].ok is True
     assert by_name["kappa"].ok is False
     assert by_name["kappa"].error == "no tree server configured"
+
+
+# ---------------------------------------------------------------------------
+# _merge_h5_files — the remote-backend result must fold into an existing file
+# ---------------------------------------------------------------------------
+
+
+def test_merge_h5_files_preserves_existing_channels(tmp_path):
+    """Regression: a remote pull used to rsync the cluster file OVER the local shot
+    file, so a one-signal custom pull (e.g. 'betan') destroyed every previously
+    fetched channel. The staged file must merge in, replacing on name conflicts and
+    keeping everything else."""
+    import h5py
+
+    from magnetics.data.fetch.toksearch import _merge_h5_files
+
+    t = np.linspace(0.0, 10.0, 100)
+    a = Channel("MPI66M307D", t, np.ones(100, np.float32), ok=True)
+    b = Channel("MPI66M340D", t, np.full(100, 2.0, np.float32), ok=True)
+    dst, _ = _write(tmp_path, [a, b])
+
+    stage = tmp_path / "stage"
+    stage.mkdir()
+    src = str(stage / "shot.h5")
+    t2 = np.linspace(0.0, 10.0, 50)
+    new = Channel("betan", t2, np.full(50, 1.5, np.float32), ok=True)
+    refetch = Channel("MPI66M340D", t, np.full(100, 9.0, np.float32), ok=True)
+    gone = Channel("ip", ok=False, error="no data")
+    _write_h5(
+        src,
+        123,
+        "custom",
+        "toksearch",
+        [new, refetch, gone],
+        compression="lzf",
+        tmin=None,
+        tmax=None,
+        stride=1,
+    )
+
+    _merge_h5_files(src, dst)
+
+    with h5py.File(dst, "r") as h5:
+        # pre-existing channel untouched, re-fetched channel replaced, new one added
+        assert set(h5) >= {"MPI66M307D", "MPI66M340D", "betan"}
+        np.testing.assert_array_equal(h5["MPI66M307D/data"][()], 1.0)
+        np.testing.assert_array_equal(h5["MPI66M340D/data"][()], 9.0)
+        np.testing.assert_array_equal(h5["betan/data"][()], 1.5)
+        np.testing.assert_allclose(h5["betan/time"][()], t2)
+        fetched = {x.decode() for x in h5.attrs["channels_fetched"]}
+        missing = {x.decode() for x in h5.attrs["channels_missing"]}
+        assert {"MPI66M307D", "MPI66M340D", "betan"} <= fetched
+        assert missing == {"ip"}
+        # both selection labels recorded
+        assert "custom" in str(h5.attrs["analysis"])
+
+
+def test_fetch_shot_remote_merges_instead_of_clobbering(tmp_path, monkeypatch):
+    """End-to-end dispatch regression: fetch_shot(backend='remote') with an existing
+    compatible shot file must stage the cluster result and merge it — the local file
+    keeps its earlier channels."""
+    import h5py
+
+    # existing local file with one 'earlier pull' channel, default window/stride
+    out = tmp_path / "shot_184927.h5"
+    t = np.linspace(0.0, 10.0, 100)
+    _write_h5(
+        str(out),
+        184927,
+        "both",
+        "toksearch",
+        [Channel("MPI66M307D", t, np.ones(100, np.float32), ok=True)],
+        compression="lzf",
+        tmin=None,
+        tmax=None,
+        stride=1,
+    )
+
+    def fake_run_remote(shot, analysis="both", *, local_out_dir=None, progress=None, **kw):
+        # the cluster always writes a FRESH file with only the requested selection
+        path = str(f"{local_out_dir}/shot_{shot}.h5")
+        import os
+
+        os.makedirs(local_out_dir, exist_ok=True)
+        t2 = np.linspace(0.0, 10.0, 50)
+        _write_h5(
+            path,
+            int(shot),
+            "custom",
+            "toksearch",
+            [Channel("betan", t2, np.full(50, 1.5, np.float32), ok=True)],
+            compression="lzf",
+            tmin=None,
+            tmax=None,
+            stride=1,
+        )
+        return path
+
+    from magnetics.data.fetch import remote as remote_mod
+
+    monkeypatch.setattr(remote_mod, "run_remote", fake_run_remote)
+
+    result = toksearch.fetch_shot(184927, backend="remote", raw_pointnames=["betan"], out=str(out))
+
+    assert result == str(out)
+    with h5py.File(out, "r") as h5:
+        assert "MPI66M307D" in h5, "remote pull clobbered the existing channel"
+        assert "betan" in h5
+    assert not (tmp_path / "_remote_stage" / "shot_184927.h5").exists()
+
+
+# ---------------------------------------------------------------------------
+# issue #63 — transport failures must not poison the incremental cache
+# ---------------------------------------------------------------------------
+
+
+def test_transport_failures_are_not_cached_as_missing(tmp_path):
+    """A channel that failed because the tunnel died must be retried on the next
+    incremental pull: it may appear in the run's report, but never in the file's
+    channels_missing attr. Genuine no-data channels ARE cached."""
+    import h5py
+
+    from magnetics.data.fetch.toksearch import _existing_channels
+
+    t = np.linspace(0.0, 10.0, 100)
+    got = Channel("GOOD", t, np.ones(100, np.float32), ok=True)
+    dead = Channel("DROPPED", ok=False, error="Connection reset by peer", error_kind="transport")
+    empty = Channel("EMPTY", ok=False, error="no data")  # default kind: no_data
+
+    out, (ok_ch, missing) = _write(tmp_path, [got, dead, empty])
+
+    # the run report still surfaces both failures...
+    assert {c.name for c in missing} == {"DROPPED", "EMPTY"}
+    with h5py.File(out, "r") as h5:
+        cached = {x.decode() for x in h5.attrs["channels_missing"]}
+    # ...but only the genuine no-data channel is cached as missing
+    assert cached == {"EMPTY"}
+    # so the next same-window pull retries DROPPED and skips GOOD + EMPTY
+    assert _existing_channels(out, None, None, 1) == {"GOOD", "EMPTY"}
+
+
+def test_error_kind_classifier():
+    from magnetics.data.fetch.toksearch import _error_kind
+
+    assert _error_kind(ConnectionResetError("reset")) == "transport"
+    assert _error_kind(TimeoutError()) == "transport"
+    assert _error_kind(OSError(32, "Broken pipe")) == "transport"
+    assert _error_kind(RuntimeError("ssh tunnel exited")) == "transport"
+    assert _error_kind(RuntimeError("%TREE-W-NODATA, no data available")) == "no_data"
+    assert _error_kind(RuntimeError("TreeNNF: node not found")) == "no_data"
diff --git a/uv.lock b/uv.lock
index 7c6e96c..472a7d1 100644
--- a/uv.lock
+++ b/uv.lock
@@ -11,7 +11,7 @@ resolution-markers = [
 ]
 
 [options]
-exclude-newer = "2026-07-01T05:08:45.379638Z"
+exclude-newer = "2026-07-03T05:43:21.277622Z"
 exclude-newer-span = "P1W"
 
 [[package]]
@@ -1176,6 +1176,7 @@ dependencies = [
     { name = "matplotlib" },
     { name = "mdsthin" },
     { name = "numpy" },
+    { name = "pexpect" },
     { name = "platformdirs" },
     { name = "scipy" },
     { name = "uvicorn", extra = ["standard"] },
@@ -1202,6 +1203,7 @@ requires-dist = [
     { name = "matplotlib", specifier = ">=3.11.0" },
     { name = "mdsthin", specifier = ">=1.0" },
     { name = "numpy", specifier = ">=2.4.6" },
+    { name = "pexpect", specifier = ">=4.9" },
     { name = "platformdirs", specifier = ">=4" },
     { name = "scipy", specifier = ">=1.18.0" },
     { name = "uvicorn", extras = ["standard"], specifier = ">=0.49.0" },
@@ -1610,7 +1612,7 @@ name = "pexpect"
 version = "4.9.0"
 source = { registry = "https://pypi.org/simple" }
 dependencies = [
-    { name = "ptyprocess", marker = "sys_platform != 'emscripten' and sys_platform != 'win32'" },
+    { name = "ptyprocess" },
 ]
 sdist = { url = "https://files.pythonhosted.org/packages/42/92/cc564bf6381ff43ce1f4d06852fc19a2f11d180f23dc32d9588bee2f149d/pexpect-4.9.0.tar.gz", hash = "sha256:ee7d41123f3c9911050ea2c2dac107568dc43b2d3b0c7557a33212c398ead30f", size = 166450, upload-time = "2023-11-25T09:07:26.339Z" }
 wheels = [