feat(compliance): USBM RI8507/OSMRE compliance chart + reference doc

sfm/compliance.py renders the velocity-vs-frequency blasting compliance chart
Blastware draws on its Event Report:
- limit_at()/limit_curve() — the RI8507 Fig B-1 / 30 CFR 816.67 curve as data
  (Drywall 0.75 + plaster 0.50 lines): 0.030in low-freq bound, plateau, 0.008in
  rising diagonal to a 2.0 in/s cap at ~40 Hz, drawn continuous.
- channel_compliance_points() — the per-cycle (freq, peak-velocity) scatter by
  the zero-crossing method (matches Blastware; cloud ceiling = channel PPV).
- draw_compliance_chart() — matplotlib rendering (both lines + scatter, BW tick
  scales + channel markers).

Verified against 7 BE12844 Blastware reports. docs/ri8507_compliance_curve.md
captures the curve construction, the SHM basis, and the scatter method.

Not yet wired into report_pdf.py — that placeholder is the next step.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01YDXjZCr4RqT2U3QvMDhgzf

docs/ri8507_compliance_curve.md

@@ -0,0 +1,135 @@
+# USBM RI8507 / OSMRE Blasting Compliance Curve — Reference
+
+Reference for the **velocity-vs-frequency blasting compliance chart** Blastware
+draws on its Event Report ("USBM RI8507 And OSMRE"), and how seismo-relay
+reproduces it. Implemented in [`sfm/compliance.py`](../sfm/compliance.py); the
+spectral (FFT) side lives in [`waveform_fft.py`](../waveform_fft.py).
+
+Reverse-engineered 2026-09-14 against 7 BE12844 (MiniMate Plus) events, each
+with a Blastware Event Report + FFT Report as ground truth. Curve values from
+USBM RI8507 Appendix B and 30 CFR 816.67.
+
+---
+
+## What it is
+
+Two closely-related sources for the same limit curve:
+
+- **USBM RI8507** — Bureau of Mines *Report of Investigations 8507* (Siskind
+  et al., 1980), *"Structure Response and Damage Produced by Ground Vibration
+  From Surface Mine Blasting."* The curve is **Figure B-1**, Appendix B
+  ("Alternative Blasting Level Criteria"), p.73–74.
+- **OSMRE / OSM** — the Office of Surface Mining Reclamation and Enforcement
+  codified it as **30 CFR 816.67, Figure 1**. "CFR" = the U.S. Code of Federal
+  Regulations. Same curve, regulatory force.
+
+The chart plots each geophone channel's significant vibration cycles as
+`(frequency, peak velocity)` points against this limit. A point **below** the
+line passes; **above** fails.
+
+---
+
+## The limit curve
+
+A structure has a resonance band (~4–12 Hz for whole structures) where it is
+most vulnerable, so the safe velocity is **lower** at those frequencies and
+**higher** away from them. The curve captures this by alternating two kinds of
+bound:
+
+- **Constant-velocity** segments — a flat horizontal line at a fixed PPV.
+- **Constant-displacement** segments — a fixed peak *displacement* `d`. For
+  simple harmonic motion, peak velocity `v = 2πf·d`, so on a velocity-vs-
+  frequency **log-log** plot this is a straight line of slope +1 (velocity rises
+  with frequency). This is why the low- and high-frequency bounds are sloped.
+
+### Two lines — structure type
+
+RI8507 gives two lines for two interior-wall constructions (Table 13, p.67):
+
+| line | construction | plateau PPV |
+|---|---|---|
+| **Drywall** (solid) | modern gypsum wallboard | **0.75 in/s** |
+| **Plaster** (dashed) | older plaster on wood lath | **0.50 in/s** |
+
+Plaster-on-lath is more damage-prone, hence the lower limit. You apply **one**
+line depending on the monitored structure.
+
+### The four segments (Figure B-1, p.74)
+
+Going low → high frequency, each line is:
+
+1. **Ultimate low-frequency bound** — constant displacement **0.030 in**
+   (`v = 2πf·0.030`). Only relevant below ~4 Hz.
+2. **Plateau** — constant velocity **0.75** (Drywall) / **0.50** (plaster) in/s.
+3. **Rising diagonal** — constant displacement **0.008 in** (`v = 2πf·0.008`),
+   climbing from the plateau up to the high-frequency cap.
+4. **High-frequency cap** — constant velocity **2.0 in/s** above ~40 Hz.
+
+The segments are drawn **continuous**: each bound is used over the frequency
+range where it is the binding (lowest) limit, and consecutive bounds meet where
+they are equal — so there are no vertical steps. Transition frequencies come
+straight from the values (`f = V / (2π·d)`):
+
+| transition | formula | Drywall | Plaster |
+|---|---|---|---|
+| 0.030 in → plateau | `V_mid / (2π·0.030)` | 3.98 Hz | 2.65 Hz |
+| plateau → 0.008 in | `V_mid / (2π·0.008)` | 14.92 Hz | 9.95 Hz |
+| 0.008 in → 2.0 in/s | `2.0 / (2π·0.008)` | 39.79 Hz | 39.79 Hz |
+
+Because both lines share the same **0.008 in** rising diagonal, above ~15 Hz
+they lie on the *same* line (both reach 2.0 in/s at ~40 Hz) — RI8507's literal
+construction merges them there. Blastware renders the dashed line as a separate
+parallel diagonal, but that is cosmetic: above ~15 Hz both structure types carry
+the identical limit, so compliance is unaffected.
+
+> ⚠ RI8507's *Table 13* is a simpler two-range criterion with a **sharp
+> discontinuity at 40 Hz** (flat plateau, then a jump to 2.0). Figure B-1 is the
+> **smoothed** version that adds the 0.008 in transition — that is the one drawn
+> on reports and implemented here.
+
+---
+
+## The compliance scatter (the points)
+
+The cloud is **not** the FFT spectrum. It is a per-cycle, time-domain measure by
+the **zero-crossing method** (`channel_compliance_points`):
+
+- Split the channel's waveform at its zero crossings.
+- Each half-cycle contributes one point: **frequency** `= 1 / (2 · half-period)`
+  (from the samples between the two crossings), **velocity** `= peak |amplitude|`
+  in that half-cycle.
+
+This yields ~90–110 points per channel, and — by construction — each channel's
+**highest** point equals that channel's PPV. Verified against Blastware: the
+cloud shape, density, and ceiling all match.
+
+### Why not the FFT?
+
+A broadband blast spreads its energy across many FFT bins, so no single bin
+reaches the time-domain peak — the FFT amplitudes come out ~10× below the
+compliance-chart velocities. The compliance chart is a *per-cycle peak* view;
+the **FFT** is a separate analysis (Blastware's *FFT Report*), reproduced by
+[`waveform_fft.py`](../waveform_fft.py) and used for the dominant-frequency
+readout and the #10 FFT view — not for this scatter.
+
+---
+
+## Implementation
+
+- `sfm/compliance.py`
+  - `limit_at(freq, curve)` — the limit PPV at a frequency (`curve` = `"Drywall"`
+    or `"Plaster"`); curves are data in `_CURVES`, so more standards can be added.
+  - `channel_compliance_points(samples, sps)` — the zero-crossing scatter.
+  - `draw_compliance_chart(ax, channels, sps)` — matplotlib rendering (both
+    limit lines + per-channel scatter, Blastware's tick scales and channel
+    markers: Tran `+` red, Vert `×` green, Long `o` blue).
+- Tests: `tests/test_compliance.py`.
+
+---
+
+## Sources
+
+- USBM **RI8507** (Siskind, Stagg, Kopp, Dowding, 1980), Appendix B / Figure B-1,
+  p.73–74; Table 13, p.67. (`ref-stuff/usbm-ri8507-ground_vibration.pdf`.)
+- **30 CFR 816.67**, "Use of explosives: Control of adverse effects," Figure 1 —
+  <https://www.ecfr.gov/current/title-30/chapter-VII/subchapter-K/part-816/section-816.67>

sfm/compliance.py

@@ -0,0 +1,132 @@
+"""USBM RI8507 / OSMRE blasting compliance chart.
+
+Renders the velocity-vs-frequency compliance scatter Blastware draws on its Event
+Report: each channel's significant waveform cycles as ``(frequency, peak
+velocity)`` points on log-log axes against the regulatory limit curve(s). A point
+below the curve passes; above fails.
+
+Two pieces, kept separate so both can be reused/extended:
+  * ``limit_at`` / ``limit_curve`` — the regulatory limit curve(s), as data.
+  * ``channel_compliance_points`` — the per-cycle (freq, velocity) scatter, by
+    the zero-crossing method (matches Blastware: each channel's cloud tops out
+    at that channel's PPV).
+
+Limit curves (USBM RI8507 Figure B-1 / OSM 30 CFR 816.67), drawn CONTINUOUS — a
+constant-displacement bound (sloped, ``v = 2πf·d``) meets a constant-velocity
+plateau at the frequency where they're equal, so there are no vertical steps
+(matching how Blastware draws it). Two lines:
+  * **Drywall** (modern gypsum board) — 0.75 in/s plateau (solid).
+  * **Plaster** on wood lath (older homes) — 0.50 in/s plateau (dashed).
+Both use a 0.030 in low-frequency displacement bound and rise through a 0.010 in
+displacement bound to a 2.0 in/s high-frequency plateau. Values from USBM RI8507
+(Appendix B) / 30 CFR 816.67; ⚠ confirm the exact shape against a Blastware
+report before trusting for compliance.
+"""
+from __future__ import annotations
+
+import math
+from typing import Dict, Sequence, Tuple
+
+import numpy as np
+from matplotlib.ticker import FixedLocator, NullLocator
+
+# curve name → (low-freq "ultimate" displacement in, mid velocity plateau in/s,
+#               high-freq displacement in, high-freq velocity plateau in/s).
+# RI8507 Fig B-1 (p.74): ultimate max displacement 0.030 in (< ~4 Hz), plateau
+# 0.75 (Drywall) / 0.50 (plaster), rising diagonal at 0.008 in displacement up to
+# a 2.0 in/s plateau reached at ~40 Hz.
+_CURVES: Dict[str, Tuple[float, float, float, float]] = {
+    "Drywall": (0.030, 0.75, 0.008, 2.00),
+    "Plaster": (0.030, 0.50, 0.008, 2.00),
+}
+# how each curve is stroked on the chart
+_CURVE_STYLE = {"Drywall": {"ls": "-", "lw": 1.0}, "Plaster": {"ls": "--", "lw": 0.9}}
+
+STANDARDS = tuple(_CURVES)
+
+# Blastware's channel markers/colours on the compliance chart.
+_CHANNEL_STYLE = {
+    "Tran": ("+", "#d62728"),   # red +
+    "Vert": ("x", "#2ca02c"),   # green x
+    "Long": ("o", "#1f77b4"),   # blue o
+}
+
+
+def limit_at(freq_hz: float, curve: str = "Drywall") -> float:
+    """Max allowed PPV (in/s) at ``freq_hz`` for ``curve`` (continuous)."""
+    d_low, v_mid, d_high, v_high = _CURVES[curve]
+    f = max(freq_hz, 1.0)
+    f_a = v_mid / (2.0 * math.pi * d_low)    # disp_low → vel_mid
+    f_b = v_mid / (2.0 * math.pi * d_high)   # vel_mid → disp_high
+    f_c = v_high / (2.0 * math.pi * d_high)  # disp_high → vel_high
+    if f <= f_a:
+        return 2.0 * math.pi * f * d_low
+    if f <= f_b:
+        return v_mid
+    if f <= f_c:
+        return 2.0 * math.pi * f * d_high
+    return v_high
+
+
+def limit_curve(curve: str = "Drywall", fmin: float = 1.0, fmax: float = 100.0, n: int = 400):
+    """(freqs, limits) sampled across the band for plotting one curve."""
+    freqs = np.logspace(np.log10(fmin), np.log10(fmax), n)
+    return freqs, np.array([limit_at(f, curve) for f in freqs])
+
+
+def channel_compliance_points(
+    samples: Sequence[float], sps: float, fmin: float = 1.0, fmax: float = 100.0,
+    vmin: float = 0.0,
+) -> Tuple[np.ndarray, np.ndarray]:
+    """Per-cycle (frequency, peak velocity) scatter for one channel.
+
+    Zero-crossing method: split the trace at sign changes; each half-cycle
+    contributes one point at ``(1/(2·half_period), max|amplitude|)``. Matches
+    Blastware — the cloud's ceiling is the channel PPV. ``samples`` must be in the
+    velocity unit you want plotted (in/s). Points outside ``[fmin, fmax]`` or at
+    or below ``vmin`` are dropped.
+    """
+    x = np.asarray(samples, dtype=float)
+    if x.size < 3:
+        return np.empty(0), np.empty(0)
+    zc = np.where(np.diff(np.signbit(x)))[0]
+    freqs, vels = [], []
+    for a, b in zip(zc[:-1], zc[1:]):
+        half_period = (b - a) / sps
+        if half_period <= 0:
+            continue
+        freqs.append(1.0 / (2.0 * half_period))
+        vels.append(float(np.abs(x[a:b + 1]).max()))
+    f = np.array(freqs)
+    v = np.array(vels)
+    keep = (f >= fmin) & (f <= fmax) & (v > vmin)
+    return f[keep], v[keep]
+
+
+def draw_compliance_chart(ax, channels: Dict[str, Sequence[float]], sps: float) -> None:
+    """Draw the compliance chart (both limit curves + per-channel scatter)."""
+    for name, style in _CURVE_STYLE.items():
+        cf, cv = limit_curve(name)
+        ax.plot(cf, cv, color="#333", zorder=3, **style)
+
+    for ch, (marker, color) in _CHANNEL_STYLE.items():
+        samples = channels.get(ch)
+        if samples is None or len(samples) == 0:
+            continue
+        f, v = channel_compliance_points(samples, sps)
+        ax.scatter(f, v, marker=marker, s=12, c=color, linewidths=0.7, zorder=4, label=ch)
+
+    ax.set_xscale("log")
+    ax.set_yscale("log")
+    ax.set_xlim(1, 100)
+    ax.set_ylim(0.0394, 10)
+    xt = [1, 2, 5, 10, 20, 50, 100]
+    yt = [0.0394, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10]
+    ax.xaxis.set_major_locator(FixedLocator(xt)); ax.xaxis.set_minor_locator(NullLocator())
+    ax.yaxis.set_major_locator(FixedLocator(yt)); ax.yaxis.set_minor_locator(NullLocator())
+    ax.set_xticklabels([str(v) for v in xt])
+    ax.set_yticklabels([("%g" % v) for v in yt])
+    ax.set_xlabel("Frequency (Hz)", fontsize=7)
+    ax.set_ylabel("Velocity (in/s)", fontsize=7)
+    ax.tick_params(labelsize=6)
+    ax.grid(True, which="both", ls=":", lw=0.4, color="#ccc")

tests/test_compliance.py

@@ -0,0 +1,35 @@
+"""USBM/OSMRE compliance curve + scatter logic (sfm.compliance).
+Rendering is verified visually against Blastware reports."""
+import math
+
+import numpy as np
+
+from sfm.compliance import limit_at, channel_compliance_points
+
+
+def test_osmre_velocity_segments():
+    assert abs(limit_at(6.0) - 0.75) < 1e-9      # 3.5–12 Hz flat
+    assert abs(limit_at(50.0) - 2.00) < 1e-9     # 30–100 Hz flat
+
+
+def test_displacement_segments():
+    assert abs(limit_at(2.0) - 2 * math.pi * 2.0 * 0.030) < 1e-9    # low-freq 0.030 in
+    assert abs(limit_at(20.0) - 2 * math.pi * 20.0 * 0.008) < 1e-9  # rising diagonal 0.008 in
+
+
+def test_limit_clamps_below_1hz():
+    assert limit_at(0.1) == limit_at(1.0)
+
+
+def test_scatter_ceiling_is_ppv_at_dominant_freq():
+    # ~27 Hz blast-like trace whose energy peaks mid-record (inside full cycles,
+    # as a real event does): the scatter cloud's ceiling is the trace PPV and the
+    # top point sits near the dominant frequency.
+    sps, n = 1024.0, 3328
+    t = np.arange(n) / sps
+    env = np.exp(-((t - 1.5) ** 2) / (2 * 0.3 ** 2))
+    x = 0.9 * env * np.sin(2 * np.pi * 27.0 * t)
+    f, v = channel_compliance_points(x, sps)
+    assert len(f) > 20
+    assert v.max() >= 0.99 * np.abs(x).max()
+    assert 20.0 < f[int(np.argmax(v))] < 35.0