def field_transport_contributions(values, *, axis=0, cancellation_rtol=1e-8):
    """Signed field totals and fractions; cancellation is explicit, never clipped.

    Parameters
    ----------
    values : array-like
        Signed spectral contributions or aligned time samples [quantity].
    axis : int, optional
        Axis containing spectral bins or additive field contributions [-].
    cancellation_rtol : float, optional
        Numerical threshold for total versus absolute summed magnitude [-].

    Returns
    -------
    dict
        Signed total, unbounded signed fractions and cancellation mask [quantity; fractions -].

    Applicability
    -------------
    Machine-independent. Requires the declared spectral convention, retained axes and units.
    These are numerical reductions of solver data, not independent evidence of physical validity.

    Notes
    -----
    Signed field totals and fractions; cancellation is explicit, never clipped.

    Provenance
    ----------
    .. [1] Canonical mathematical definition in #2028; uses the retained-source and weighting contract of #2025.
    """
    values = np.asarray(values, float)
    total = np.sum(values, axis=axis)
    magnitude = np.sum(np.abs(values), axis=axis)
    cancelled = np.abs(total) <= cancellation_rtol * magnitude
    good = np.isfinite(total) & (total != 0) & ~cancelled
    fractions = np.divide(values, np.expand_dims(np.where(good, total, np.nan), axis))
    return {
        "total": total,
        "fractions": fractions,
        "cancellation": cancelled,
        "definition": "signed contribution / signed total",
    }

def prepare_diagnostic_traces(diagnostics, frame_times):
    """Prepare causal, gap-broken traces; renderers never interpolate or reduce.

    Parameters
    ----------
    diagnostics : mapping
        Labelled traces with time, values, unit and optional maximum gap [declared trace units].
    frame_times : array-like
        Simulation timestamps of the selected field frames [trace time units].

    Returns
    -------
    dict
        Causal per-frame traces, broken gaps and fixed display limits [input trace units].

    Applicability
    -------------
    Machine-independent. Requires the declared physical definitions, source identity and units.
    Matching display conventions alone does not establish physical equivalence or convergence.

    Notes
    -----
    Prepare causal, gap-broken traces; renderers never interpolate or reduce.

    Provenance
    ----------
    .. [issue] #2030 four-axis comparison and time-selection contract; #2028 diagnostic ownership.
    """
    frame_times = np.asarray(frame_times, float)
    if frame_times.ndim != 1 or not np.isfinite(frame_times).all():
        raise ValueError("frame times must be finite and one-dimensional")
    result = {}
    for label, trace in diagnostics.items():
        t = np.asarray(trace["time"], float)
        y = np.asarray(trace["values"], float)
        if (
            t.ndim != 1
            or y.shape != t.shape
            or not len(t)
            or not np.isfinite(t).all()
            or np.any(np.diff(t) <= 0)
        ):
            raise ValueError(f"invalid diagnostic trace {label}")
        gap = trace.get(
            "max_gap", 3 * float(np.median(np.diff(t))) if len(t) > 1 else np.inf
        )
        gap = float(gap)
        if np.isnan(gap) or gap <= 0:
            raise ValueError("diagnostic max_gap must be positive and not NaN")
        frames = []
        for instant in frame_times:
            n = int(np.searchsorted(t, instant, side="right"))
            x = t[:n]
            v = y[:n].copy()
            # Insert separators; retain the first observed sample after each gap.
            if n > 1:
                breaks = np.flatnonzero(np.diff(x) > gap) + 1
                x = np.insert(x, breaks, np.nan)
                v = np.insert(v, breaks, np.nan)
            frames.append((x, v))
        finite = y[np.isfinite(y)]
        low, high = (
            (float(finite.min()), float(finite.max())) if finite.size else (0.0, 1.0)
        )
        pad = max((high - low) * 0.05, abs(high) * 0.01, 1e-12)
        result[label] = {
            "frames": frames,
            "limits": (low - pad, high + pad),
            "unit": trace.get("unit", "unspecified"),
            "max_gap": gap,
        }
    return result

def radial_correlation_from_spectrum(kx, spectrum, *, box_length, samples=513):
    """Cosine autocorrelation 1/e crossing and integral scale; unresolved is a bound.

    Parameters
    ----------
    kx : array-like
        Radial Fourier wavenumbers [inverse perpendicular length].
    spectrum : array-like
        Nonnegative radial power spectrum on kx [intensity].
    box_length : float
        Positive periodic radial box extent [perpendicular length].
    samples : int, optional
        Number of separation samples across half the box; default 513 [-].

    Returns
    -------
    dict
        1/e crossing, lower bound if unresolved, and first-zero integral scale [perpendicular length].

    Applicability
    -------------
    Machine-independent. Requires the declared spectral convention, retained axes and units.
    These are numerical reductions of solver data, not independent evidence of physical validity.

    Notes
    -----
    Cosine autocorrelation 1/e crossing and integral scale; unresolved is a bound.

    Array-only counterpart of the CGYRO locality reader, with its existing
    half-box search and first-zero integral definition preserved.

    Provenance
    ----------
    .. [1] Existing CGYRO locality algorithm (#1354/#1484), spectrum-based cosine autocorrelation; #2028.
    """
    kx, power = np.asarray(kx, float), np.asarray(spectrum, float)
    length = float(box_length)
    if (
        kx.ndim != 1
        or power.shape != kx.shape
        or not np.isfinite(kx).all()
        or not np.isfinite(power).all()
        or np.any(power < 0)
        or not np.isfinite(length)
        or length <= 0
        or samples < 2
    ):
        raise ValueError("finite kx/nonnegative spectrum and positive box required")
    if not np.any(power > 0):
        return {"reason": "zero turbulent amplitude in the window"}
    dx = np.linspace(0, 0.5 * length, samples)
    correlation = (power[None, :] * np.cos(np.outer(dx, kx))).sum(axis=1)
    correlation /= correlation[0]
    below = np.flatnonzero(correlation <= np.exp(-1.0))
    crossing = None
    if below.size:
        j = int(below[0])
        x0, x1, c0, c1 = dx[j - 1], dx[j], correlation[j - 1], correlation[j]
        crossing = float(x0 + (np.exp(-1.0) - c0) * (x1 - x0) / (c1 - c0))
    zero = np.flatnonzero(correlation <= 0)
    stop = int(zero[0]) if zero.size else dx.size
    integrate = getattr(np, "trapezoid", None) or np.trapz
    return {
        "l_corr": crossing,
        "l_corr_lower_bound": None if crossing is not None else 0.5 * length,
        "l_integral": float(integrate(correlation[:stop], dx[:stop])),
    }

def rebin_contributions(values, edges, new_edges):
    """Conserve signed flux using piecewise-uniform density within audited bins.

    Parameters
    ----------
    values : array-like
        Signed spectral contributions or aligned time samples [quantity].
    edges : array-like, optional
        Audited increasing spectral bin edges [wavenumber].
    new_edges : array-like
        Destination bin edges spanning exactly the source interval [wavenumber].

    Returns
    -------
    numpy.ndarray
        Conservative destination-bin contributions under piecewise-uniform density [quantity].

    Applicability
    -------------
    Machine-independent. Requires the declared spectral convention, retained axes and units.
    These are numerical reductions of solver data, not independent evidence of physical validity.

    Notes
    -----
    Conserve signed flux using piecewise-uniform density within audited bins.

    The destination must cover exactly the same interval; no extrapolation.

    Provenance
    ----------
    .. [1] Canonical mathematical definition in #2028; uses the retained-source and weighting contract of #2025.
    """
    values, edges, new = map(lambda x: np.asarray(x, float), (values, edges, new_edges))
    if (
        values.shape[-1] + 1 != len(edges)
        or len(new) < 2
        or not np.isfinite(edges).all()
        or not np.isfinite(new).all()
        or np.any(np.diff(edges) <= 0)
        or np.any(np.diff(new) <= 0)
        or not np.allclose([edges[0], edges[-1]], [new[0], new[-1]], rtol=0, atol=1e-12)
    ):
        raise ValueError("audited increasing edges with identical coverage required")
    overlap = np.maximum(
        0,
        np.minimum(new[1:, None], edges[None, 1:])
        - np.maximum(new[:-1, None], edges[None, :-1]),
    )
    return (values / np.diff(edges)) @ overlap.T

def saturated_intensity(ql_weights, potential_intensity):
    """TGLF moment intensity = QL intensity weight * saturated phi intensity.

    Parameters
    ----------
    ql_weights : array-like
        Native TGLF quasilinear moment weights (species, ky, mode, moment) [native weight].
    potential_intensity : array-like
        Saturated potential intensity (ky, mode) [native phi squared].

    Returns
    -------
    numpy.ndarray
        Saturated moment intensity on the original four axes [native moment intensity].

    Applicability
    -------------
    Machine-independent. Requires the declared spectral convention, retained axes and units.
    These are numerical reductions of solver data, not independent evidence of physical validity.

    Notes
    -----
    TGLF moment intensity = QL intensity weight * saturated phi intensity.

    Weights: (species, ky, mode, moment); phi: (ky, mode). This is not a
    transport-flux reconstruction: transport also needs its audited quadrature.

    Provenance
    ----------
    .. [1] GACODE b493397 tglf_inout.f90 field/intensity writers; native SAT0/SAT2 fixture identities; #2025/#2028.
    """
    q, p = np.asarray(ql_weights, float), np.asarray(potential_intensity, float)
    if q.ndim != 4 or p.shape != q.shape[1:3]:
        raise ValueError("QL weights and saturated phi intensity axes disagree")
    return q * p[None, :, :, None]

def select_time(times, requested=None, *, tolerance=None):
    """Nearest finite time and evidence; no implicit scientific tolerance.

    Parameters
    ----------
    times : array-like
        Stored finite, distinct experimental timestamps [s].
    requested : float, optional
        Requested experimental time, or first stored time when omitted [s].
    tolerance : float, optional
        Explicit maximum time mismatch; absence leaves scientific alignment unverified [s].

    Returns
    -------
    dict
        Stored index, requested/selected time, difference and alignment status [s; index -].

    Applicability
    -------------
    Machine-independent. Requires the declared physical definitions, source identity and units.
    Matching display conventions alone does not establish physical equivalence or convergence.

    Notes
    -----
    Nearest finite time and evidence; no implicit scientific tolerance.

    Provenance
    ----------
    .. [issue] #2030 four-axis comparison and time-selection contract; #2028 diagnostic ownership.
    """
    t = np.asarray(times, float)
    if (
        t.ndim != 1
        or not len(t)
        or not np.isfinite(t).all()
        or len(np.unique(t)) != len(t)
    ):
        raise ValueError("unavailable or ambiguous time stamps")
    if tolerance is not None and (not np.isfinite(tolerance) or tolerance < 0):
        raise ValueError("time tolerance must be finite and nonnegative")
    target = float(t[0] if requested is None else requested)
    if not np.isfinite(target):
        raise ValueError("requested time must be finite")
    distances = np.abs(t - target)
    i = int(np.argmin(distances))
    delta = float(t[i] - target)
    if (
        np.count_nonzero(
            np.isclose(
                distances,
                distances[i],
                rtol=0,
                atol=4 * np.finfo(float).eps * max(1, abs(target)),
            )
        )
        > 1
    ):
        raise ValueError("ambiguous equidistant time slices")
    if tolerance is not None and abs(delta) > tolerance:
        raise ValueError("selected experimental time outside tolerance")
    return {
        "index": i,
        "requested_time": target,
        "selected_time": float(t[i]),
        "delta_time": delta,
        "tolerance": tolerance,
        "status": "matched" if tolerance is not None else "unverified",
    }

def spectral_integral(values, *, representation="contribution", edges=None, axis=-1):
    """Integrate signed bin contributions or densities with explicit bin edges.

    Parameters
    ----------
    values : array-like
        Signed spectral contributions or aligned time samples [quantity].
    representation : str
        Contribution means bin-integrated values; density requires edges [-].
    edges : array-like, optional
        Audited increasing spectral bin edges [wavenumber].
    axis : int, optional
        Axis containing spectral bins or additive field contributions [-].

    Returns
    -------
    numpy.ndarray
        Signed integral with the selected spectral axis removed [quantity].

    Applicability
    -------------
    Machine-independent. Requires the declared spectral convention, retained axes and units.
    These are numerical reductions of solver data, not independent evidence of physical validity.

    Notes
    -----
    Integrate signed bin contributions or densities with explicit bin edges.

    Provenance
    ----------
    .. [1] Canonical mathematical definition in #2028; uses the retained-source and weighting contract of #2025.
    """
    values = np.asarray(values, float)
    if representation == "contribution":
        return np.sum(values, axis=axis)
    if representation != "density" or edges is None:
        raise ValueError("density requires audited bin edges")
    edges = np.asarray(edges, float)
    if (
        edges.size != values.shape[axis] + 1
        or not np.isfinite(edges).all()
        or np.any(np.diff(edges) <= 0)
    ):
        raise ValueError("increasing finite edges must bracket every bin")
    shape = [1] * values.ndim
    shape[axis] = len(edges) - 1
    return np.sum(values * np.diff(edges).reshape(shape), axis=axis)

def time_window_average(time, values, *, window, max_gap=None):
    """Average aligned scalar or multi-channel histories with the shared gap policy.

    Parameters
    ----------
    time : array-like
        Strictly increasing sample times [time].
    values : array-like
        Histories with time on the final axis [quantity].
    window : tuple of float
        Selected integration bounds in the input time unit [time].
    max_gap : float, optional
        Largest integrated interval; default three median sample spacings [time].

    Returns
    -------
    numpy.ndarray
        Means with the time axis removed; unavailable integrals are NaN [quantity].

    Applicability
    -------------
    Machine-independent. Uses observed adjacent samples only, with the same
    no-extrapolation and missing-interval policy as reduce_time_series.

    Provenance
    ----------
    .. [1] Shared piecewise-linear reduction contract in #2028.
    """
    t, y = np.asarray(time, float), np.asarray(values, float)
    if y.ndim < 1 or y.shape[-1] != t.size:
        raise ValueError("histories must end in the simulation time axis")
    means = [
        reduce_time_series(t, row, window=window, max_gap=max_gap, blocks=2)["mean"]
        for row in y.reshape(-1, len(t))
    ]
    return np.asarray([np.nan if v is None else v for v in means]).reshape(y.shape[:-1])

def zonal_fraction(phi, ky, *, weights=None, ky_tolerance=1e-12, weighting=None):
    """I_Z/(I_Z+I_T); all-zero intensity is undefined, not zero.

    Parameters
    ----------
    phi : array-like
        Complex potential coefficients (kx, ky, time) [native potential].
    ky : array-like
        Binormal wavenumbers, including zero [inverse perpendicular length].
    weights : array-like, optional
        Nonnegative intensity weights broadcast over kx and ky [-].
    ky_tolerance : float, optional
        Numerical threshold identifying the zero binormal wavenumber [ky unit].
    weighting : str, optional
        Explicit interpretation of the supplied intensity weights [-].

    Returns
    -------
    numpy.ndarray
        Zonal fraction at each time; all-zero intensity is NaN [-].

    Applicability
    -------------
    Machine-independent. Requires the declared spectral convention, retained axes and units.
    These are numerical reductions of solver data, not independent evidence of physical validity.

    Notes
    -----
    I_Z/(I_Z+I_T); all-zero intensity is undefined, not zero.

    Provenance
    ----------
    .. [1] Canonical mathematical definition in #2028; uses the retained-source and weighting contract of #2025.
    """
    result = zonal_turbulence_intensity(
        phi, ky, weights=weights, ky_tolerance=ky_tolerance, weighting=weighting
    )
    return result.fraction

def zonal_shear_proxy(phi, kx, ky, *, ky_tolerance=1e-12):
    r"""Spectral proxies for the zonal-flow energy and the zonal shear.

    $E_Z \propto \sum_{k_x} k_x^2|\phi_{k_x,0}|^2$ and
    $S_Z^2 \propto \sum_{k_x} k_x^4|\phi_{k_x,0}|^2$, from
    $v_{E,Z} \propto k_x\phi_Z$ and $\partial_x v_{E,Z} \propto k_x^2\phi_Z$.

    Parameters
    ----------
    phi : array-like
        Potential Fourier components, shape ``(n_kx, n_ky, n_time)`` [solver units].
    kx : array-like
        Radial wavenumber of each ``kx`` row, length ``n_kx`` [1/length].
    ky : array-like
        Binormal wavenumber of each ``ky`` column, length ``n_ky`` [1/length].
    ky_tolerance : float, optional
        Largest ``|ky|`` counted as zonal [1/length].

    Returns
    -------
    ZonalShearProxy
        ``flow`` and ``shear`` per time sample and the ``normalization`` label
        [kx^2 phi^2, kx^4 phi^2].

    Raises
    ------
    ValueError
        Shapes do not match, inputs are not finite, or no zonal ``ky`` exists.

    Defaults
    --------
    ``ky_tolerance=1e-12`` is a numerical convenience (solvers store the zonal
    wavenumber as exact zero).

    Convention
    ----------
    Spectral proxies in the units of $\phi$ and $k_x$: physical flow energy
    and shearing rate need the solver's potential and perpendicular-length
    normalisation and geometric factors, which this routine does not apply.

    Applicability
    -------------
    Machine-independent.

    Provenance
    ----------
    .. [1] P. H. Diamond, S.-I. Itoh, K. Itoh and T. S. Hahm,
           Plasma Phys. Control. Fusion 47 (2005) R35.
    .. [issue] #1820 Sec. 7.
    """
    cube = _spectral_cube(phi)
    zonal_mask, _ = _ky_masks(ky, cube.shape[1], ky_tolerance)
    k = np.asarray(kx, dtype=float).ravel()
    if k.size != cube.shape[0] or not np.all(np.isfinite(k)):
        raise ValueError(f"kx must be finite with {cube.shape[0]} entries")
    power = (np.abs(cube[:, zonal_mask]) ** 2).sum(axis=1)       # (n_kx, n_time)
    flow = (k[:, None] ** 2 * power).sum(axis=0)
    shear = (k[:, None] ** 4 * power).sum(axis=0)
    return ZonalShearProxy(flow=flow, shear=shear,
                           normalization="spectral proxy: sum kx^p |phi_{kx,0}|^2 in solver units")

def zonal_turbulence_intensity(phi, ky, *, weights=None, ky_tolerance=1e-12, weighting=None):
    r"""Split fluctuation intensity into its zonal ($k_y = 0$) and non-zonal parts.

    $I_\mathrm{turb}(t) = \sum_{k_x, k_y\ne0} w\,|\phi_{k_x,k_y}(t)|^2$,
    $I_\mathrm{zonal}(t) = \sum_{k_x} w\,|\phi_{k_x,0}(t)|^2$ and
    $R_{ZT} = I_\mathrm{zonal}/I_\mathrm{turb}$.

    Parameters
    ----------
    phi : array-like
        Potential Fourier components, shape ``(n_kx, n_ky, n_time)``, real or
        complex; average or select any poloidal coordinate beforehand [solver units].
    ky : array-like
        Binormal wavenumber of each ``ky`` column; the zonal column is the one
        with ``|ky| <= ky_tolerance`` [1/length].
    weights : array-like, optional
        Spectral weight broadcastable to ``(n_kx, n_ky)``; ``None`` for 1 [caller units].
    ky_tolerance : float, optional
        Largest ``|ky|`` counted as zonal [1/length].
    weighting : str, optional
        Name of the weighting, recorded in the result [-].

    Returns
    -------
    ZonalTurbulenceIntensity
        ``turbulence``, ``zonal`` and ``ratio`` per time sample (``ratio`` is
        NaN where the turbulence intensity is zero) and the ``weighting`` label
        [weight * phi^2; ratio -].

    Raises
    ------
    ValueError
        ``phi`` is not 3-D or not finite, ``ky`` does not match its ky axis, no
        ``ky`` is zonal within the tolerance, or a weight is negative.

    Defaults
    --------
    ``ky_tolerance=1e-12`` is a numerical convenience: solvers store the zonal
    wavenumber as exact zero.  Unit weights (``weights=None``) are the
    conventional potential-intensity proxy.

    Convention
    ----------
    Unweighted $|\phi|^2$ is a potential-intensity proxy, labelled as such; it
    becomes an energy only with the metric weights the caller supplies (e.g.
    $k_\perp^2$ for an $E\times B$ kinetic-energy proxy).  Solvers that store
    only $k_y \ge 0$ count each non-zonal mode once.

    Applicability
    -------------
    Machine-independent.

    Provenance
    ----------
    .. [1] S. Kobayashi, O. D. Gurcan and P. H. Diamond, Phys. Plasmas 22 (2015) 090702.
    .. [issue] #1820 Sec. 6.
    """
    cube = _spectral_cube(phi)
    zonal_mask, turb_mask = _ky_masks(ky, cube.shape[1], ky_tolerance)
    power = np.abs(cube) ** 2
    if weights is not None:
        w = np.asarray(weights, dtype=float)
        if not np.all(np.isfinite(w)) or np.any(w < 0):
            raise ValueError("weights must be finite and non-negative")
        power = power * np.broadcast_to(w, cube.shape[:2])[..., None]
    turbulence = power[:, turb_mask].sum(axis=(0, 1))
    zonal = power[:, zonal_mask].sum(axis=(0, 1))
    with np.errstate(divide="ignore", invalid="ignore"):
        ratio = np.where(turbulence > 0, zonal / np.where(turbulence > 0, turbulence, 1.0), np.nan)
    label = weighting or ("unweighted |phi|^2 (potential-intensity proxy)" if weights is None
                          else "caller-weighted |phi|^2")
    return ZonalTurbulenceIntensity(turbulence=turbulence, zonal=zonal, ratio=ratio, weighting=label)