def equilibrium_safety_factor(equilibrium) -> tuple[np.ndarray, np.ndarray]:
    """The magnitude of the safety factor on normalized flux, as the equilibrium knows it.

    Parameters
    ----------
    equilibrium : EquilibriumData or source
        An equilibrium record, or anything
        :func:`vaft.process.equilibrium.as_equilibrium` accepts [-].

    Returns
    -------
    tuple of numpy.ndarray
        ``(psi_n, abs(q))``, sorted by increasing ``psi_n`` [-].

    Raises
    ------
    ValueError
        The record carries neither a q profile nor the flux map and poloidal
        current function to compute one, or its flux storage family is unknown.

    Convention
    ----------
    The record's own profile when it has one; otherwise q is computed from the
    flux map and ``F`` by :func:`calculate_q_profile_from_psi` on 97 interior
    levels, which is the case for an analytic Solov'ev export. Only the
    magnitude is returned: resonance is ``|q| = m/n``, and the direction is the
    helicity's job.

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

    Limitations
    -----------
    A computed profile starts at ``psi_n = 0.02`` and ends at ``0.98``, so a
    resonance closer to the axis or the boundary than that is not found.

    Provenance
    ----------
    .. [1] :func:`vaft.process.equilibrium.calculate_q_profile_from_psi`, the
       contour-integral safety factor it falls back to.
    """
    eq = _coerce(equilibrium)
    span = eq.psi_boundary - eq.psi_axis
    if eq.q is not None and eq.psi_1d is not None and eq.q.size == eq.psi_1d.size:
        psi_n = (eq.psi_1d - eq.psi_axis) / span
        q = np.abs(eq.q)
    else:
        if eq.psi is None or eq.f is None or eq.psi_1d is None:
            raise ValueError(
                "the equilibrium has no q profile and lacks the flux map or F to compute one"
            )
        psi_n = np.linspace(0.02, 0.98, 97)
        q = np.abs(np.asarray(calculate_q_profile_from_psi(
            eq.psi, eq.r, eq.z, (eq.psi_1d, eq.f), eq.psi_axis, eq.psi_boundary, psi_n,
            axis_rz=eq.magnetic_axis, cocos=_cocos_for_q(eq),
        ), dtype=float))
    order = np.argsort(psi_n)
    psi_n, q = psi_n[order], q[order]
    keep = np.isfinite(q)
    return psi_n[keep], q[keep]

def resolve_rational_surface(equilibrium, m: int, n: int, *, index: int = 0) -> tuple[float, float]:
    """Where the equilibrium's own safety factor equals ``m/n``.

    Parameters
    ----------
    equilibrium : EquilibriumData or source
        An equilibrium record, or anything
        :func:`vaft.process.equilibrium.as_equilibrium` accepts [-].
    m : int
        Poloidal mode number [-].
    n : int
        Toroidal mode number [-].
    index : int, optional
        Which crossing when ``q`` crosses ``m/n`` more than once, innermost
        first [-].

    Returns
    -------
    tuple of float
        ``(psi_n_s, q_s)``: the normalized flux of the resonant surface and the
        profile's value there, which equals ``m/n`` to the root tolerance [-].

    Raises
    ------
    ValueError
        ``q = m/n`` does not occur in the profile, or ``index`` asks for a
        crossing that does not exist.

    Convention
    ----------
    Resonance is ``|q| = m/n``. The crossings are bracketed by
    :func:`find_rational_surfaces` and each is refined by a root solve on a
    cubic spline through the profile, so the answer is not limited to the
    linear interpolation between samples.

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

    Limitations
    -----------
    Only as good as the profile the equilibrium carries; a reconstruction with
    a poorly constrained core q places a core resonance equally poorly.

    Provenance
    ----------
    .. [1] :func:`vaft.process.equilibrium.find_rational_surfaces` for the
       bracketing.
    """
    from scipy.interpolate import CubicSpline
    from scipy.optimize import brentq

    psi_n, q = equilibrium_safety_factor(equilibrium)
    target = float(m) / float(n)
    found = find_rational_surfaces(psi_n, q, n, m_range=(m, m))
    coarse = np.asarray(found["psi_n_rational"], dtype=float)
    coarse = coarse[np.asarray(found["m"]) == m]
    if coarse.size == 0:
        raise ValueError(
            f"q = {m}/{n} = {target:.4g} does not occur in this equilibrium: |q| spans "
            f"[{np.min(q):.4g}, {np.max(q):.4g}] over psi_N in [{psi_n[0]:.3g}, {psi_n[-1]:.3g}]"
        )
    coarse = np.sort(coarse)
    if not 0 <= index < coarse.size:
        raise ValueError(f"q = {m}/{n} crosses {coarse.size} time(s); index {index} does not exist")
    guess = float(coarse[index])
    unique, first = np.unique(psi_n, return_index=True)
    # anti-alias: not time-domain -- q over normalized flux.
    spline = CubicSpline(unique, q[first])
    j = int(np.clip(np.searchsorted(unique, guess), 1, unique.size - 1))
    lo, hi = unique[j - 1], unique[j]
    f_lo, f_hi = spline(lo) - target, spline(hi) - target
    if f_lo * f_hi > 0.0:
        return guess, float(spline(guess))
    root = float(brentq(lambda s: float(spline(s)) - target, lo, hi, xtol=1e-12))
    return root, float(spline(root))