# Probes and S-Parameters

Channel: release; package: 1.8.0; source: 7d1e97331656dd74ad0f5471653ac0c5effdff70

Page: https://remilab.ai/rfx/versions/v1.8.0/guide/probes-sparams/
Source: https://github.com/bk-squared/rfx/blob/7d1e97331656dd74ad0f5471653ac0c5effdff70/docs/public/guide/probes-sparams.mdx

rfx provides three probe types for extracting field data, plus several S-parameter extraction methods for different port models.

---

## Point probe

Records a single field component at one Yee cell over time. The time series is stored in `result.time_series`.

Setup used by the snippets below:

```python
import numpy as np
from rfx import GaussianPulse, Simulation

sim = Simulation(
    freq_max=5e9,
    domain=(0.05, 0.05, 0.05),
    dx=2.5e-3,
    boundary="cpml",
    cpml_layers=6,
)
sim.add_source(
    (0.015, 0.025, 0.025),
    "ez",
    waveform=GaussianPulse(f0=2.5e9, bandwidth=0.8),
    amplitude_kind="current",
)

# A two-port S-matrix cube for the host-side report checks below.
quality_sim = Simulation(
    freq_max=5e9,
    domain=(0.05, 0.04, 0.02),
    dx=2.5e-3,
    boundary="cpml",
    cpml_layers=6,
)
quality_sim.add_port((0.01, 0.02, 0.01), "ez", impedance=50.0)
quality_sim.add_port((0.04, 0.02, 0.01), "ez", impedance=50.0)
s_matrix = quality_sim.run(n_steps=600, compute_s_params=True).s_params
```

```python
sim.add_probe(
    position=(0.025, 0.025, 0.010),  # (x, y, z) in meters
    component="ez",                   # ex / ey / ez / hx / hy / hz
)
```

After running:

```python
result = sim.run(n_steps=4000)
ts = result.time_series          # shape (n_steps, n_probes)
ez_vs_time = ts[:, 0]            # first probe, all timesteps
```

Multiple probes are indexed in the order they were registered.

---

## Vector probe

Records all six field components (Ex, Ey, Ez, Hx, Hy, Hz) at the same position.

```python
sim.add_vector_probe((0.025, 0.025, 0.010))
```

The six columns appended by this call correspond to Ex, Ey, Ez, Hx, Hy, Hz in
that order. They are columns 0–5 only when no probe was registered earlier.

---

## DFT plane probe

Accumulates selected frequencies directly on a 2-D cross-section during the
time loop. This uses less memory than saving full snapshots and avoids
FFT-bin/interpolation constraints, but it does not remove finite-window leakage
or incomplete-settling error. Repeat with a longer run before interpreting a
narrow feature.

```python
import jax.numpy as jnp

sim.add_dft_plane_probe(
    axis="z",
    coordinate=0.010,           # z-plane at 10 mm
    component="ez",
    freqs=jnp.linspace(1e9, 4e9, 50),
    name="ez_cross_section",
)
```

Access the result:

```python
result = sim.run(n_steps=4000)
plane = result.dft_planes["ez_cross_section"]
# plane.accumulator : shape (n_freqs, nx, ny) complex  (in-plane dims depend on axis)
# plane.freqs       : shape (n_freqs,)
```

!!! tip
    Omit `freqs` and set `n_freqs` to let rfx choose a frequency array automatically
    from `freq_max / 10` to `freq_max`.

### Exactly-zero flux diagnostics

During normal Python execution outside JAX tracing, `flux_spectrum(...)` checks
float32 results that are exactly zero at every frequency even though their DFT
accumulators are nonzero. If a float64 recomputation of the same Poynting sum is
nonzero, rfx warns that the per-cell `E × H*` products underflowed and were
flushed to zero. The returned float32 array is not rewritten. Enable x64 for
that computation, increase the source amplitude, or recompute from the healthy
accumulators in float64; do not interpret the exact zeros as a physical
transmission null.

---

## S-parameter extraction

Use the calculator that matches the port family:

| Primitive | Calculator | Result | Port model / scope |
|---|---|---|---|
| `add_port(..., extent=None)` | `run(compute_s_params=True)` | `Result.s_params`, `Result.freqs` | Single-cell lumped port; S normalized to each port impedance |
| `add_port(..., extent=...)` | `run(compute_s_params=True)` | `Result.s_params`, `Result.freqs` | Multi-cell wire port; V/I DFT integrals |
| lumped/wire `add_port(...)` in AD objectives | `forward(port_s11_freqs=...)` | per-port S11 in `ForwardResult.s_params` | Differentiable S11 through `forward()` |
| `add_msl_port(...)` | `compute_msl_s_matrix()` | `MSLSMatrixResult.S`, `.freqs`, `.Z0`, `.beta` | Microstrip line; uniform mesh is externally compared; nonuniform `laplace`/`uniform` modes have internal checks only; `eigenmode` is unsupported |
| `add_waveguide_port(...)` | `compute_waveguide_s_matrix()` | `WaveguideSMatrixResult.s_params`, `.freqs`, metadata | Rectangular metal waveguide; modal V/I or power-flux |
| `add_coaxial_port(...)` | `compute_coaxial_line_reflection(...)` | `CoaxialLineReflectionResult.s11`, `.freqs`, diagnostics | float32 precision; nonperiodic 3D second-order uniform Yee grid; CPML on all six boundary faces with positive thickness on both z faces; `cpml_axes="z"`; exactly one `face="top"` coaxial port; see the [API constraints](https://remilab.ai/rfx/versions/v1.8.0/api/sources-ports/) |

Sources, TFSF, probes, DFT plane probes, and flux monitors are non-port
observables: they record fields but do not define a port impedance or
S-matrix reference. Pick the calculator by the port family you registered.

For host-side report checks on any S-matrix cube, use
`network_quality_metrics()`:

```python
from rfx import network_quality_metrics

quality = network_quality_metrics(s_matrix)
print(quality["passivity_excess"], quality["reciprocity_error"])
```

This helper reports finite-data, passivity, reciprocity, max-magnitude, and
column-power checks. These are host-side report diagnostics; they do not
alter any solver path and are not a substitute for physically validating the
S-parameters themselves.

Automatic checks differ by calculator during normal Python execution:

- lumped/wire `run()` and uniform single-device `forward()` warn on non-finite entries or an
  individual `|S| > 1.1`; they do not test total outgoing column power;
- MSL and waveguide full-matrix calculators check individual entries and
  column power with calculator-specific tolerances;
- `compute_coaxial_line_reflection()` does not run that shared passivity guard,
  so inspect `status`, `recurrence_residual`, `fit_residual`, and `s11`
  explicitly. It is a self-contained line calculation and rejects separately
  registered geometry, circuit elements, monitors, and termination helpers;
  see the [complete argument rules](https://remilab.ai/rfx/versions/v1.8.0/api/sources-ports/).

Settling is a separate guard from passivity. `compute_msl_s_matrix`,
`compute_waveguide_s_matrix`, `compute_mixed_s_matrix`,
`compute_coaxial_two_port`, and `compute_coax_msl_transition` compare each
drive's `settling_db` with the shared -40 dB threshold and warn when a fixed
record ends too early. The warning does not modify `S`, and no warning is not
proof of mesh convergence or external accuracy.
`compute_coaxial_line_reflection()` is not one of those five full-matrix lanes;
use its returned recurrence and fit diagnostics.

Warnings leave every computed value intact. If a passive model exceeds its
documented limit, treat that frequency bin as suspect and check settling time,
normalization, the reference plane, and the extraction setup. Silence only
means that an implemented check did not fire; it does not certify passivity.

### Lumped port S-matrix

When one or more lumped ports are added with `add_port()`, pass `compute_s_params=True`:

```python
# Ports drive the run themselves, so this simulation carries no soft source;
# the probe records the ring-down the Harminv section reads.
sim = Simulation(
    freq_max=5e9,
    domain=(0.05, 0.05, 0.05),
    dx=2.5e-3,
    boundary="cpml",
    cpml_layers=6,
)
sim.add_probe((0.025, 0.025, 0.025), "ez")

sim.add_port((0.01, 0.02, 0.01), "ez", impedance=50)
sim.add_port((0.04, 0.02, 0.01), "ez", impedance=50)

result = sim.run(n_steps=3000, compute_s_params=True)
s_matrix = result.s_params           # (2, 2, n_freqs) complex
freqs     = result.freqs             # (n_freqs,) Hz

import numpy as np
s11_db = 20 * np.log10(np.abs(s_matrix[0, 0, :]))
s21_db = 20 * np.log10(np.abs(s_matrix[1, 0, :]))
```

The S-matrix is normalized to each port's `impedance`.

### Wire port S-matrix (JIT-DFT)

Wire ports use voltage/current DFT integrals accumulated during the time loop. The extraction is performed inside a JIT-compiled scan, making it faster than post-processing FFTs for long simulations.

```python
# A mixed lumped+wire port set is unsupported, so the wire port gets its own
# simulation.
wire_sim = Simulation(
    freq_max=5e9,
    domain=(0.05, 0.05, 0.02),
    dx=2.5e-3,
    boundary="cpml",
    cpml_layers=6,
)

wire_sim.add_port(
    (0.01, 0.02, 0.0),
    "ez",
    impedance=50,
    extent=0.0016,               # wire length
)

wire_result = wire_sim.run(n_steps=5000, compute_s_params=True)
s11 = wire_result.s_params[0, 0, :]  # complex, shape (n_freqs,)
```

`add_port(..., reference_plane_cells=N)` (wire ports, multi-drive S-matrix
runs only) reads V/I at two planes N and 2N cells from the port along the
line. The build rejects a plane that leaves the declared domain or sits within
the absorber's own thickness of a domain face, and a plane zone that reaches
another port — including any rung of a microstrip port's de-embedding probe
ladder, at the auto-resolved offset and spacing the run actually probes. Both
crossing checks compare positions along the line axis only, so a transversely
separated port on a parallel trace also trips them.

### Microstrip-line S-matrix

MSL ports use their own calculator and result object:

```python
from rfx import Box

# A microstrip board of its own: ground, substrate, trace, and two MSL ports.
msl_sim = Simulation(
    freq_max=6e9,
    domain=(0.020, 0.010, 0.004),
    dx=0.25e-3,
    boundary="cpml",
    cpml_layers=3,
)
msl_sim.add_material("substrate", eps_r=3.2, sigma=0.01)
msl_sim.add(Box((0.75e-3, 0.75e-3, 0.75e-3), (19.25e-3, 9.25e-3, 1.25e-3)), material="pec")
msl_sim.add(Box((0.75e-3, 0.75e-3, 1.25e-3), (19.25e-3, 9.25e-3, 2.25e-3)), material="substrate")
msl_sim.add(Box((0.75e-3, 4.50e-3, 2.25e-3), (19.25e-3, 5.50e-3, 2.75e-3)), material="pec")
for x_feed, feed_direction, port_name in ((0.004, "+x", "in"), (0.016, "-x", "out")):
    msl_sim.add_msl_port(
        (x_feed, 0.005, 1.25e-3),
        width=1.0e-3,
        height=1.0e-3,
        direction=feed_direction,
        eps_r_sub=3.2,
        name=port_name,
    )

msl = msl_sim.compute_msl_s_matrix(n_freqs=101, num_periods=40)
S = msl.S          # (n_ports, n_ports, n_freqs)
freqs = msl.freqs  # (n_freqs,) Hz
Z0 = msl.Z0        # extracted characteristic impedance
reliable = msl.reliable  # (port, frequency), or None while tracing
                         # `port` indexes the PROBE PLANE, not the driven port
```

`num_periods=40` is the docstring default; values below that may under-settle
the fields. Do not expect `run(compute_s_params=True)` to include MSL ports in
`Result.s_params`.

#### Low-signal MSL bins

During normal Python execution outside JAX tracing,
`MSLSMatrixResult.reliable[p, k]` is `False` if, **in any one of the drives**,
both measured `|V|` and `|I|` on port `p`'s plane at frequency index `k` fall
below 10% of that individual `(driven, port)` record's own band median. The
mask is screened per record and then reduced over the drive axis, so `p`
identifies the probe plane to investigate — not which drive was running. This
condition can occur near a strong-reflector measurement null, but the mask does
not diagnose the cause.

The S-matrix is solved from the wave amplitudes of **all** drives at once
(`S = B·A⁻¹`), so a collapsed wave pair at any one port contaminates the
entire frequency slice `S[:, :, k]` — not just that port's column. A false
entry at any port therefore means the whole S-matrix at that frequency should
not be plotted, fitted, or optimized as physical data. The values remain
present for diagnosis.

`reliable[p, k]` is `False` when port `p`'s plane collapsed at bin `k` in at
least one drive: every driven/port record the solve consumes is covered, so the
index does attribute which plane to investigate.

A true entry certifies that the low-signal threshold did not fire — never that
the result is accurate.

The threshold is relative to each record's own band median, so a port sitting in
a deep stopband is not flagged wholesale — but individual deep bins are. Live
extractor runs on rfx's two filter geometries flagged 2 bins of 100 on the
microstrip notch fixture and 12 of 120 on the Sheen low-pass leg — and the notch
fixture's two are the notch centre itself, 3.6273 GHz, which its committed
fixture records at −30.66 dB. The two counts are not recomputable from the
committed fixtures, which store S magnitudes only with no V/I dump; checking them
means re-running the extractor and reading `reliable`. That is not a false alarm: at a
−30 dB notch the passive port's wave split really is low-signal, so the extractor
cannot certify the depth there. It does mean a filter measurement loses its most
interesting bin to the screen; read the depth from `.S_raw` or a flux monitor
with that caveat rather than treating the screen as free.

The filter below is therefore the right per-bin screen — it keeps exactly the
bins where no plane the solve reads had collapsed:

```python
import numpy as np

if msl.reliable is not None:
    usable_bins = np.all(msl.reliable, axis=0)
    S_usable = msl.S[..., usable_bins]
    freqs_usable = msl.freqs[usable_bins]
```

### Waveguide S-parameters

See [Sources & Ports](https://remilab.ai/rfx/versions/v1.8.0/guide/sources-ports/) for the waveguide-port setup workflow.
Use `compute_waveguide_s_matrix(...)` for the full multi-port matrix:

```python
# Waveguide ports cannot share a simulation with lumped ports.
wg_sim = Simulation(
    freq_max=11e9,
    domain=(0.10, 0.023, 0.010),
    dx=1e-3,
    boundary="cpml",
    cpml_layers=8,
)
for x_face, face_direction, port_name in ((0.01, "+x", "port1"), (0.09, "-x", "port2")):
    wg_sim.add_waveguide_port(
        x_position=x_face,
        y_range=(0.0, 0.023),
        z_range=(0.0, 0.010),
        direction=face_direction,
        f0=9e9,
        bandwidth=0.5,
        name=port_name,
    )

wg = wg_sim.compute_waveguide_s_matrix(num_periods=40, normalize=True)  # 40 = docstring default; 30 under-settles (-35.8 dB vs the -40 dB bar)
S = wg.s_params      # (n_ports, n_ports, n_freqs)
freqs = wg.freqs
```

`normalize=True` is suitable for transmission measurements. On a uniform mesh,
use `normalize=False` or `normalize="flux"` for strong-reflector S11 (for
example, a PEC short) to avoid the reference-run subtraction artifact. On a
non-uniform mesh (`dx_profile` / `dy_profile` / `dz_profile`), `normalize=False`
is unsupported;
use `normalize="flux"` for strong-reflector S11. See
[Migration](https://remilab.ai/rfx/versions/v1.8.0/guide/migration/) and
[Support Boundaries](https://remilab.ai/rfx/versions/v1.8.0/api/support-boundaries/).

For memory-heavy waveguide AD runs, pass `checkpoint_segments=K` to segment
the reverse-mode tape (peak backward memory drops from O(n_steps) toward
O(√n_steps) near K≈√n_steps). On the uniform path, `K` must divide the timestep
count exactly. On a non-uniform (`dx_profile` / `dy_profile` / `dz_profile`) path, `K` is a
target segment count and need not divide the timestep count: rfx selects the
divisor chunk size nearest to `n_steps / K` and applies checkpointing to the
device run. This non-uniform AD path requires `normalize="flux"` with an
`eps_override` or `sigma_override` design variable; checkpointing does not
expand the RF validation scope.

`run(...)` can also expose per-port waveguide diagnostics:

```python
# run() supports a single waveguide port, so drive one port here.
diag_sim = Simulation(
    freq_max=11e9,
    domain=(0.10, 0.023, 0.010),
    dx=1e-3,
    boundary="cpml",
    cpml_layers=8,
)
diag_sim.add_waveguide_port(
    x_position=0.01,
    y_range=(0.0, 0.023),
    z_range=(0.0, 0.010),
    f0=9e9,
    bandwidth=0.5,
    name="port1",
)
wg_run = diag_sim.run(n_steps=3000)

sp = wg_run.waveguide_sparams["port1"]

print(sp.freqs.shape)              # (n_freqs,)
print(sp.s11.shape)                # (n_freqs,) complex
print(sp.calibration_preset)
print(f"Reference plane: {sp.reference_plane*1e3:.1f} mm")
print(f"Probe plane:     {sp.probe_plane*1e3:.1f} mm")
```

`WaveguideSParamResult` fields:

| Field | Type | Description |
|---|---|---|
| `freqs` | ndarray | Frequency array (Hz) |
| `s11` | complex ndarray | Reflection coefficient at port |
| `s21` | complex ndarray | Transmission to measurement plane |
| `calibration_preset` | str | `"measured"`, `"source_to_probe"`, or `"explicit"` |
| `source_plane` | float | Actual source injection plane (m) |
| `reference_plane` | float | Reported S11 reference plane (m) |
| `probe_plane` | float | Reported S21 measurement plane (m) |

---

## Resonance extraction — Harminv

`result.find_resonances()` runs harmonic inversion (the Matrix Pencil Method, Sarkar & Pereira 1995) on the probe ring-down signal to extract complex resonant frequencies. It returns a list of `HarminvMode` objects sorted by amplitude (strongest first).

```python
modes = result.find_resonances(
    freq_range=(1e9, 5e9),   # search band (Hz)
    probe_idx=0,             # which probe to analyse (default 0)
)

for m in modes:
    print(f"f={m.freq/1e9:.4f} GHz  Q={m.Q:.0f}  amp={m.amplitude:.3e}")
```

`HarminvMode` fields:

| Field | Description |
|---|---|
| `freq` | Resonant frequency (Hz) |
| `decay` | Exponential decay rate (1/s) |
| `Q` | Quality factor = pi*f/decay |
| `amplitude` | Mode amplitude (magnitude) |
| `phase` | Phase (radians) |
| `error` | Relative fit-error estimate for the mode |

### Parameters

```python
modes = result.find_resonances(
    freq_range=(1e9, 4e9),
    probe_idx=0,
    source_decay_time=None,  # heuristic from the requested frequency range; see note
    bandpass=None,           # auto: True for CPML (removes DC artifacts), False for PEC
)
```

!!! note
    Harminv requires the signal to have decayed to ring-down before the analysis
    window starts. `source_decay_time` is subtracted from the beginning of the
    time series. When it is `None`, rfx does **not** inspect the registered
    source: it uses the midpoint of `freq_range` as `f_center`, a fixed
    `bandwidth=0.8`, and `2 * (3 / (pi * f_center * bandwidth))`. Pass the
    actual source decay time when its center frequency, bandwidth, waveform, or
    turn-off differs from that heuristic.

### Standalone Harminv

You can call Harminv directly on any time series:

```python
from rfx import harminv
import numpy as np

dt = result.dt                              # timestep in seconds
signal = np.asarray(result.time_series[:, 0])
modes = harminv(signal, dt, f_min=1e9, f_max=5e9)
```

Standalone `harminv(...)` uses `decimate="auto"` by default. When the sample
rate is more than eight times `f_max`, it applies staged, anti-aliased FIR
decimation before the matrix-pencil solve while preserving the record duration
and using the effective timestep. This removes redundant oversampling from a
calculation whose cost grows steeply with sample count. Pass `decimate=False`
when exact full-rate reproduction is required.

---

## Plotting helpers

Each helper takes arrays (or the final `state` + `grid`) rather than the
`Result` object, and returns a matplotlib `Figure`:

```python
from rfx import plot_s_params, plot_time_series, plot_field_slice

# Magnitude (dB) of every S_ij vs frequency
plot_s_params(result.s_params, result.freqs)

# All probe time series on one axes
plot_time_series(result.time_series, result.dt)

# 2-D Ez slice through the final field state
plot_field_slice(result.state, result.grid, component="ez", axis="z", index=30)
```

---

## Touchstone I/O

Touchstone export uses rfx's canonical S-matrix shape `(n_ports, n_ports, n_freqs)`. Reading returns `(s_params, freqs, z0)`.

```python
from rfx import write_touchstone, read_touchstone, read_touchstone_full

# Write 2-port S-params to a .s2p file
write_touchstone("patch_antenna.s2p", result.s_params, result.freqs, z0=50.0)

# Read back (legacy tuple API)
s_matrix, freqs, z0 = read_touchstone("patch_antenna.s2p")

# Metadata-aware Touchstone 2.0 read path preserves per-port references
write_touchstone("device_v2.s2p", result.s_params, result.freqs,
                 version="2.0", layout="standard", port_z0=[50.0, 50.0])
network = read_touchstone_full("device_v2.s2p")
print(network.reference)
```