Probes and S-Parameters

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.

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:

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.

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.

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:

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

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():

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.

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:

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.

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

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

Microstrip-line S-matrix

MSL ports use their own calculator and result object:

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:

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 for the waveguide-port setup workflow. Use compute_waveguide_s_matrix(...) for the full multi-port matrix:

wg = sim.compute_waveguide_s_matrix(num_periods=30, normalize=True)
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 dx_profile / dy_profile mesh, normalize=False is unsupported; use normalize="flux" for strong-reflector S11. See Migration and 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 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:

sp = result.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).

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

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:

from rfx import harminv
import numpy as np

dt = 1e-12  # timestep in seconds
signal = np.loadtxt("probe_data.txt")
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:

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).

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
network = read_touchstone_full("device_v2.s4p")
print(network.reference)