diff --git a/research/duetsch_re100_kc5/capture.py b/research/duetsch_re100_kc5/capture.py
new file mode 100644
index 0000000..f6a6266
--- /dev/null
+++ b/research/duetsch_re100_kc5/capture.py
@@ -0,0 +1,52 @@
+"""Research-only field capture helpers for the Dütsch driver (numpy only).
+
+Snapshots are taken at the first step of the scored phase bins in declared
+cycles: a near-body crop at full resolution and a fluid-weighted block-mean
+downsample of the whole domain; plus scalar trends at every bin read.
+"""
+from __future__ import annotations
+
+import numpy as np
+
+CROP_X, CROP_Y, DOWN = 8.0, 6.0, 8      # crop half-widths in D; downsample factor
+STRIP = 2.0                             # inlet/outlet strip width in D
+
+
+def crop_box(xm, ym, D, n, body_extent):
+    """Inclusive integer crop (x0, x1, y0, y1); fail closed if it leaves the domain
+    interior or does not contain the body's swept extent (x_lo, x_hi, y_lo, y_hi)."""
+    x0, x1 = int(round(xm - CROP_X * D)), int(round(xm + CROP_X * D))
+    y0, y1 = int(round(ym - CROP_Y * D)), int(round(ym + CROP_Y * D))
+    if x0 < 1 or y0 < 1 or x1 > n - 2 or y1 > n - 2:
+        raise ValueError("capture crop leaves the domain interior")
+    bx0, bx1, by0, by1 = body_extent
+    if not (x0 <= bx0 and bx1 <= x1 and y0 <= by0 and by1 <= y1):
+        raise ValueError("capture crop clips the body's swept extent")
+    return x0, x1, y0, y1
+
+
+def downsample(field, fluid, f=DOWN):
+    """Fluid-weighted f x f block means; NaN where a block has no fluid cell."""
+    ny, nx = field.shape
+    h, w = ny // f, nx // f
+    v = np.where(fluid, field, 0.0)[:h * f, :w * f].reshape(h, f, w, f).sum(axis=(1, 3))
+    c = fluid[:h * f, :w * f].reshape(h, f, w, f).sum(axis=(1, 3))
+    with np.errstate(invalid="ignore", divide="ignore"):
+        return np.where(c > 0, v / np.maximum(c, 1), np.nan).astype(np.float32)
+
+
+def fluid_mean(field, fluid, xs=slice(None), ys=slice(None)):
+    m = fluid[ys, xs]
+    return float(field[ys, xs][m].mean()) if m.any() else float("nan")
+
+
+def trend_row(s, k, ph, rho, ux, uy, fluid, D, n, crop):
+    """(s, cycle, phase, <rho>, <ux>, <uy>, <rho>_crop, inlet <rho>, <ux>, outlet <rho>, <ux>)."""
+    w = max(1, int(round(STRIP * D)))
+    x0, x1, y0, y1 = crop
+    cx, cy = slice(x0, x1 + 1), slice(y0, y1 + 1)
+    inlet, outlet, inner = slice(1, 1 + w), slice(n - 1 - w, n - 1), slice(1, n - 1)
+    return (s, k, ph, fluid_mean(rho, fluid), fluid_mean(ux, fluid), fluid_mean(uy, fluid),
+            fluid_mean(rho, fluid, cx, cy),
+            fluid_mean(rho, fluid, inlet, inner), fluid_mean(ux, fluid, inlet, inner),
+            fluid_mean(rho, fluid, outlet, inner), fluid_mean(ux, fluid, outlet, inner))
diff --git a/research/duetsch_re100_kc5/run.py b/research/duetsch_re100_kc5/run.py
index 9e9bc15..1278de7 100644
--- a/research/duetsch_re100_kc5/run.py
+++ b/research/duetsch_re100_kc5/run.py
@@ -49,6 +49,8 @@ def main():
     ap.add_argument("--T", type=int, default=None, help="period in steps (multiple of 360)")
     ap.add_argument("--cycles", type=int, default=p.LAST_SCORED, help="run to the end of this cycle")
     ap.add_argument("--out", required=True)
+    ap.add_argument("--dump-fields", default="",
+                    help="research capture: comma-separated paper cycles, e.g. 11,21,30,31")
     args = ap.parse_args()
     D, T = args.D, args.T or 126 * args.D
     if T % 360:
@@ -69,6 +71,13 @@ def main():
         body = sim.add_body("circle", (xm + A, ym), R,
                             translation=Translation((-2.0 * A, 0.0, 0.0, 0.0), 1.01 * u_max))
         sim.initialize()
+        dump = {int(c) for c in args.dump_fields.split(",") if c}
+        if dump:
+            import capture
+            if max(dump) > args.cycles or min(dump) <= p.DISCARD:
+                raise SystemExit("--dump-fields cycles must lie in the scored range of this run")
+            crop = capture.crop_box(xm, ym, D, n, (xm - A - R, xm + A + R, ym - R, ym + R))
+            trend = []
         steps = int(round((0.25 + args.cycles) * T))
         force = np.zeros((steps, 2)); diag = []; samples = []
         for m in range(steps):
@@ -83,6 +92,20 @@ def main():
                 fields = sim.get_macroscopic()
                 fluid = (sim.get_flags().reshape(n, n) & 1) != 0
                 ux, uy = (np.asarray(fields[k]).reshape(n, n) for k in ("ux", "uy"))
+                if dump:  # research capture only; the default path is unchanged
+                    rho = np.asarray(fields["rho"]).reshape(n, n)
+                    trend.append(capture.trend_row(s, p.cycle(s, T), ph, rho, ux, uy, fluid, D, n, crop))
+                    if p.cycle(s, T) in dump and p.phase(s, T) == ph:  # first step of the bin
+                        x0c, x1c, y0c, y1c = crop
+                        cy_, cx_ = slice(y0c, y1c + 1), slice(x0c, x1c + 1)
+                        np.savez_compressed(
+                            os.path.join(args.out, f"fields_c{p.cycle(s, T)}_p{int(ph)}.npz"),
+                            s=s, cycle=p.cycle(s, T), phase=ph, crop=np.asarray(crop),
+                            force=force[m], centre=np.asarray((xm + p.position(s, T, A), ym)),
+                            ux=ux[cy_, cx_].astype(np.float32), uy=uy[cy_, cx_].astype(np.float32),
+                            rho=rho[cy_, cx_].astype(np.float32), fluid=fluid[cy_, cx_].astype(np.uint8),
+                            down_ux=capture.downsample(ux, fluid), down_uy=capture.downsample(uy, fluid),
+                            down_rho=capture.downsample(rho, fluid))
                 centre = (xm + p.position(s, T, A), ym)  # analytic centre after the step
                 for xs in p.STATIONS:
                     for x2 in p.X2:
@@ -92,6 +115,8 @@ def main():
                         samples.append((s, p.cycle(s, T), ph, xs, x2,
                                         math.nan if u is None else u, math.nan if v is None else v,
                                         wsum, json.dumps(wu)))
+        if dump:
+            np.save(os.path.join(args.out, "trend.npy"), np.asarray(trend, float))
         np.savez(os.path.join(args.out, "raw.npz"), force=force, diag=np.asarray(diag),
                  samples=np.asarray([r[:8] for r in samples], float))
         with open(os.path.join(args.out, "case.json"), "w") as fh:
diff --git a/tests/unit/test_duetsch_protocol.py b/tests/unit/test_duetsch_protocol.py
index 16c0259..cdde4fe 100644
--- a/tests/unit/test_duetsch_protocol.py
+++ b/tests/unit/test_duetsch_protocol.py
@@ -150,3 +150,42 @@ def test_horizon_excludes_future_coverage_and_events(monkeypatch):
     entry = json.loads(out[0])["cases"][0]
     assert sorted(int(k) for k in entry["coverage"]) == [3, 4]   # nothing from cycles 5, 6
     assert entry["no_donor"] == 0.0
+
+
+
+def _capture():
+    spec = importlib.util.spec_from_file_location(
+        "duetsch_capture", Path(__file__).resolve().parents[2] / "research/duetsch_re100_kc5/capture.py")
+    mod = importlib.util.module_from_spec(spec)
+    spec.loader.exec_module(mod)
+    return mod
+
+
+def test_capture_crop_downsample_and_trends():
+    c = _capture()
+    n, D = 2400, 40
+    xm = ym = 0.5 * (n - 1)
+    A, R = p.amplitude(D), 20.0
+    box = c.crop_box(xm, ym, D, n, (xm - A - R, xm + A + R, ym - R, ym + R))
+    assert box == (880, 1520, 960, 1440) and box[1] - box[0] + 1 == 641 and box[3] - box[2] + 1 == 481
+    with pytest.raises(ValueError):
+        c.crop_box(xm, ym, D, 600, (0, 1, 0, 1))                    # leaves the domain
+    with pytest.raises(ValueError):
+        c.crop_box(xm, ym, D, n, (xm - 9 * D, xm, ym, ym))           # clips the swept body
+    field = np.arange(64.0).reshape(8, 8)
+    fluid = np.ones((8, 8), bool)
+    fluid[:4, :4] = False
+    down = c.downsample(field, fluid, 4)
+    assert np.isnan(down[0, 0]) and down[1, 1] == field[4:, 4:].mean()
+    rho = np.ones((n, n)); rho[:, 1:81] = 1.01
+    row = c.trend_row(10, 3, 180.0, rho, np.zeros((n, n)), np.zeros((n, n)),
+                      np.ones((n, n), bool), D, n, box)
+    assert row[7] == pytest.approx(1.01) and row[9] == pytest.approx(1.0)
+
+
+def test_dump_steps_are_the_first_step_of_each_bin():
+    T = 5040
+    for k in (11, 21, 30, 31):
+        starts = [s for s in range(T // 4 + (k - 1) * T, T // 4 + k * T)
+                  if p.bin_phase(s, T) is not None and p.phase(s, T) == p.bin_phase(s, T)]
+        assert len(starts) == 3 and all(p.cycle(s, T) == k for s in starts)