def amplification_ratio(
    psi_norm,
    values,
    reference_psi_norm,
    reference_values,
    *,
    window: ResonantWindow,
    statistic: str = "rms",
) -> float:
    """How much larger one run's resonant response is than another's.

    Parameters
    ----------
    psi_norm : array_like
        Rational-surface coordinate of the case [-].
    values : array_like
        The case's resonant column [any].
    reference_psi_norm : array_like
        Rational-surface coordinate of the reference, usually the vacuum run
        [-].
    reference_values : array_like
        The reference's resonant column, in the same unit as ``values`` [any].
    window : ResonantWindow
        The interval both are reduced over [-].
    statistic : str, optional
        One of :data:`RESONANT_STATISTICS` [n/a].

    Returns
    -------
    float
        Case over reference; ``nan`` when either window is empty or the
        reference reduces to zero [-].

    Processing steps
    ----------------
    1. Reduce each run over ``window`` independently, so the two need not
       share a rational-surface set.
    2. Divide.

    Limitations
    -----------
    Both runs are reduced on their *own* surfaces. Two equilibria resonate at
    different places, and interpolating one onto the other's surfaces would
    invent a resonance where there is none.

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

    Provenance
    ----------
    .. [legacy] ``gpec_analysis.calculate_amplification_ratio``, which fell
       back to the vacuum column silently when the plasma one was absent --
       returning a ratio of one and calling it no amplification. Not carried
       over: an absent column raises here.
    """
    case = reduce_resonant(psi_norm, values, window=window, statistic=statistic)
    reference = reduce_resonant(
        reference_psi_norm, reference_values, window=window, statistic=statistic
    )
    if not np.isfinite(case) or not np.isfinite(reference) or reference == 0.0:
        return float("nan")
    return float(case / reference)

def chirikov(psi_norm, width, *, definition: str = "surface") -> np.ndarray:
    """The island-overlap parameter, by either of the two definitions in use.

    Parameters
    ----------
    psi_norm : array_like
        Rational-surface positions [-].
    width : array_like
        Full island width per surface [-].
    definition : str, optional
        ``"surface"`` centres it on each surface, ``width / d_nn``, and
        returns one value per surface. ``"pair"`` takes each adjacent pair,
        ``(w_inner + w_outer) / 2 / spacing``, and returns one value per
        *pair* -- one fewer than the surfaces [n/a].

    Returns
    -------
    ndarray
        For ``"surface"``, one value per surface in the caller's order. For
        ``"pair"``, one per adjacent pair **sorted outward** -- which two
        islands are neighbours is a property of the sorted sequence, and
        :class:`IslandPairs` carries the indices that map each pair back [-].

    Raises
    ------
    ValueError
        ``definition`` is neither, or two surfaces share a position.

    Defaults
    --------
    ``definition="surface"`` is a legacy compatibility value: it is the form
    GPEC itself reports as ``K_isl``, so it is what a caller comparing with a
    file will want. It also returns a different array length from ``"pair"``.

    Convention
    ----------
    The two definitions answer different questions and disagree by design:
    the surface form asks whether one island reaches its nearest neighbour,
    the pair form whether two particular islands reach each other. GPEC's own
    ``K_isl`` is the surface form.

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

    Provenance
    ----------
    .. [legacy] ``gpec_multimode_metrics.compute_surface_chirikov`` and
       ``compute_pair_chirikov``.
    """
    if definition not in ("surface", "pair"):
        raise ValueError(f"definition must be 'surface' or 'pair', not {definition!r}")
    if definition == "surface":
        return penetration_ratio(psi_norm, width)
    return island_pairs(psi_norm, width).chirikov

def critical_island_width(psi_norm) -> np.ndarray:
    """The island width at which a surface's islands reach its neighbour's.

    Purely geometric: an island overlaps its neighbour when it is as wide as
    the distance to that neighbour, so the critical width *is* the
    nearest-neighbour spacing.

    Parameters
    ----------
    psi_norm : array_like
        Rational-surface positions [-].

    Returns
    -------
    ndarray
        Critical width per surface, in the same units as ``psi_norm`` [-].

    Convention
    ----------
    Widths are full widths, in normalized poloidal flux, matching what GPEC
    writes as ``w_isl``. The identity this rests on -- ``w_isl / w_crit``
    equalling GPEC's own ``K_isl`` -- holds because GPEC forms its Chirikov
    parameter from a half width over a half distance, and the halves cancel.

    Limitations
    -----------
    Undefined for a single surface, which has no neighbour: ``nan``.

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

    Provenance
    ----------
    .. [gpec] Reproduces GPEC's own ``K_isl`` exactly when divided into
       ``w_isl``: zero relative difference on every surface of the DIII-D
       n=1, n=2 and n=3 ideal runs.
    .. [ideal] GPEC leaves its own ``w_isl_v_crit`` at zero on an ideal run,
       which is why this is computed rather than read.
    """
    return nearest_surface_spacing(psi_norm)

def dominant_mode(coupling, *, response_area=None, control_area=None) -> DominantMode:
    """The leading direction of a coupling matrix, weighted by surface area.

    Parameters
    ----------
    coupling : array_like
        Response-by-control coupling, complex [n/a].
    response_area : array_like, optional
        Area per response row. Unweighted when omitted [m^2].
    control_area : array_like, optional
        Area per control column. Unweighted when omitted [m^2].

    Returns
    -------
    DominantMode
        The leading singular value, and the control and response vectors in
        the unweighted basis [n/a].

    Raises
    ------
    ValueError
        The matrix is not 2-D, an area does not match its dimension, or an
        area is not finite and positive.

    Convention
    ----------
    Rows are weighted by ``sqrt(area)`` and columns divided by it, so the
    decomposition is of ``W_r M W_c**-1`` and the vectors are mapped back
    through the weights. The returned control vector is the right singular
    vector itself -- ``vh[0].conj()`` -- and a projection onto it is
    ``np.vdot(control_vector, x)``, which conjugates it; that is not the
    ``vh.conj() @ x`` defect :func:`edge_overlap_metric` once had.

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

    Processing steps
    ----------------
    1. Check the areas **before** taking their square root.
    2. Weight rows and columns, decompose, map the leading pair back.

    Provenance
    ----------
    .. [legacy] ``resistive_gpec_metrics.py::dominant_mode``. Its guard,
       ``np.any(np.sqrt(area) <= 0)``, is defeated by a negative area:
       ``sqrt`` of one is ``nan``, ``nan <= 0`` is ``False``, and the error
       surfaced later as ``SVD did not converge`` rather than naming the area.
    """
    matrix = np.asarray(coupling, dtype=complex)
    if matrix.ndim != 2:
        raise ValueError(f"the coupling must be 2-D, not {matrix.ndim}-D")
    if not np.all(np.isfinite(matrix)):
        raise ValueError("the coupling holds a non-finite value")

    def weights(area, size, label):
        if area is None:
            return np.ones(size)
        values = np.atleast_1d(np.asarray(area, dtype=float))
        if values.shape != (size,):
            raise ValueError(f"{values.size} {label} areas against {size} {label}s")
        # Checked on the area, not on its root: sqrt of a negative is nan,
        # and nan <= 0 is False.
        if not np.all(np.isfinite(values)) or np.any(values <= 0.0):
            raise ValueError(f"every {label} area must be finite and positive")
        return np.sqrt(values)

    rows = weights(response_area, matrix.shape[0], "response")
    columns = weights(control_area, matrix.shape[1], "control")
    weighted = rows[:, None] * matrix / columns[None, :]
    left, singular, right = np.linalg.svd(weighted, full_matrices=False)
    return DominantMode(
        singular_value=float(singular[0]),
        control_vector=right[0].conj() / columns,
        response_vector=left[:, 0] / rows,
    )

def helical_phase_sweep(
    tables: Mapping[int, Mapping[str, object]],
    q: float,
    column: str,
    *,
    angles,
    weights: Mapping[int, complex] | None = None,
    integrality_tolerance: float = 1e-6,
) -> np.ndarray:
    """The composite drive at one surface, as a function of the helical angle.

    Parameters
    ----------
    tables : mapping of int to mapping
        One resonant table per toroidal mode number [n/a].
    q : float
        The rational surface to sweep at [-].
    column : str
        The column to sum [n/a].
    angles : array_like
        Helical angles to evaluate at [rad].
    weights : mapping of int to complex, optional
        Per-mode complex weight, applied before the sum [n/a].
    integrality_tolerance : float, optional
        Passed to :func:`align_surfaces_by_q` [-].

    Returns
    -------
    ndarray
        The real amplitude at each angle, in the order given. Units follow
        ``column`` [n/a].

    Raises
    ------
    ValueError
        As :func:`composite_drive_at_q`, or an angle is not finite.

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

    Convention
    ----------
    The angle is the helical one, ``u = q theta - zeta``, in which the mode
    ``n`` advances as ``exp(i n u)``; :func:`composite_drive_at_q` derives
    it. Its origin is the input's, and shifting that origin translates the
    whole curve rather than changing its shape.

    Processing steps
    ----------------
    1. Align, weight and collect the contributions, as
       :func:`composite_drive_at_q` does.
    2. Evaluate ``sum_n Re(z_n exp(i n u))`` at each requested angle.

    Limitations
    -----------
    A sweep resolves a series whose highest harmonic is the largest
    resonating ``n``; sampling it more coarsely than a few points per period
    of that harmonic will miss the peak, and nothing here checks the caller's
    spacing against it.

    Provenance
    ----------
    .. [legacy] ``gpec_multimode_metrics.phase_sweep_composite``, which
       rotated one nominated mode through a default three angles -- 0, 90 and
       180 degrees, so a sweep could not tell ``+90`` from ``-90`` -- and
       returned ``nan`` rows rather than an error when no surface matched.
       Rotating one mode and sweeping the shared angle differ once more than
       two modes resonate.
    """
    requested = np.atleast_1d(np.asarray(angles, dtype=float))
    if not np.all(np.isfinite(requested)):
        raise ValueError("every angle must be finite")
    modes, values = _contributions_at_q(
        tables, q, column, weights=weights, integrality_tolerance=integrality_tolerance
    )
    orders = np.asarray(modes, dtype=float)
    return (values[None, :] * np.exp(1j * np.outer(requested, orders))).real.sum(axis=1)

def island_pairs(psi_norm, width) -> IslandPairs:
    """Every adjacent pair of islands, with the gap between their separatrices.

    Parameters
    ----------
    psi_norm : array_like
        Rational-surface positions [-].
    width : array_like
        Full island width per surface [-].

    Returns
    -------
    IslandPairs
        One entry per adjacent pair, sorted outward; empty for fewer than two
        surfaces [-].

    Raises
    ------
    ValueError
        ``width`` is not one per surface, two surfaces share a position, or a
        position or width is not finite or a width is negative.

    Processing steps
    ----------------
    1. Sort outward.
    2. For each pair, the gap is the spacing less the two half widths.
    3. A gap at or below zero -- within rounding -- is contact.

    Defaults
    --------
    The contact tolerance is a numerical convenience, 1e-12 of the spacing.

    Convention
    ----------
    Residual gaps below that fraction of the spacing count as contact rather
    than separation, so two separatrices that meet exactly are not reported
    as apart -- and "round-off" means the same thing at every radius, which a
    fixed absolute tolerance does not.

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

    Provenance
    ----------
    .. [legacy] ``gpec_multimode_metrics.compute_pair_chirikov``.
    """
    ordered, ordered_width, order = _island_arrays(psi_norm, width)
    if ordered.size < 2:
        empty = np.empty(0, dtype=float)
        index = np.empty(0, dtype=int)
        return IslandPairs(index, index, empty, empty, empty, empty, empty,
                           empty, empty, np.empty(0, dtype=bool))
    spacing = np.diff(ordered)
    half_sum = 0.5 * (ordered_width[:-1] + ordered_width[1:])
    gap = spacing - half_sum
    return IslandPairs(
        inner_index=order[:-1],
        outer_index=order[1:],
        inner_psi_norm=ordered[:-1],
        outer_psi_norm=ordered[1:],
        spacing=spacing,
        inner_right=ordered[:-1] + 0.5 * ordered_width[:-1],
        outer_left=ordered[1:] - 0.5 * ordered_width[1:],
        gap=gap,
        chirikov=half_sum / spacing,
        # Scaled to the spacing, so "round-off" means the same thing whether
        # the surfaces are a thousandth apart or a whole unit. A fixed atol --
        # what the code this replaces used, with an rtol against zero that
        # could never contribute -- is a different criterion at every radius.
        overlaps=gap <= _CONTACT_TOLERANCE * spacing,
    )

def jump_width(lundquist, model: JumpWidthModel) -> np.ndarray:
    """The resistive-layer width a Lundquist number implies.

    Parameters
    ----------
    lundquist : array_like
        Lundquist number per rational surface [-].
    model : JumpWidthModel
        The width law. **Required**: the coefficients are one paper's fit,
        not a property of the physics, so there is no default to fall back
        on [n/a].

    Returns
    -------
    ndarray
        Full layer width per surface, in normalized poloidal flux [-].

    Raises
    ------
    ValueError
        A Lundquist number is not finite and positive.

    Convention
    ----------
    ``w = floor + scale * S**(-1/3)``. The ``S**(-1/3)`` scaling is the
    resistive-inertial layer's; the floor keeps the width finite as ``S``
    grows, where the pure scaling would shrink the layer below the grid.

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

    Processing steps
    ----------------
    1. Check the Lundquist numbers are finite and positive.
    2. Apply the model's law.

    Provenance
    ----------
    .. [legacy] ``resistive_gpec_metrics.py::jump_width_from_lundquist``,
       which carried the paper's coefficients as keyword defaults -- so a
       caller who never chose a model got one anyway.
    """
    if not isinstance(model, JumpWidthModel):
        raise TypeError(
            f"model must be a JumpWidthModel, not {type(model).__name__}; "
            "PAPER_JUMP_WIDTH reproduces the legacy scan"
        )
    values = np.atleast_1d(np.asarray(lundquist, dtype=float))
    if not np.all(np.isfinite(values)) or np.any(values <= 0.0):
        raise ValueError("every Lundquist number must be finite and positive")
    return model.floor + model.scale * np.power(values, -1.0 / 3.0)

def lundquist_number(tau_resistive, tau_alfven) -> np.ndarray:
    """The Lundquist number, the resistive over the Alfven time.

    Parameters
    ----------
    tau_resistive : array_like
        Resistive diffusion time per rational surface [s].
    tau_alfven : array_like
        Alfven time per rational surface [s].

    Returns
    -------
    ndarray
        ``S = tau_r / tau_a`` per surface [-].

    Raises
    ------
    ValueError
        The two disagree in shape, or a time is not finite and positive.

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

    Processing steps
    ----------------
    1. Check both times are finite, positive and the same shape.
    2. Divide.

    Provenance
    ----------
    .. [legacy] ``library/resistive_gpec_metrics.py::lundquist_number``
       (branch ``codex/transp``).
    """
    resistive = np.atleast_1d(np.asarray(tau_resistive, dtype=float))
    alfven = np.atleast_1d(np.asarray(tau_alfven, dtype=float))
    if resistive.shape != alfven.shape:
        raise ValueError(f"{resistive.shape} resistive times against {alfven.shape} Alfven times")
    for label, values in (("resistive", resistive), ("Alfven", alfven)):
        if not np.all(np.isfinite(values)) or np.any(values <= 0.0):
            raise ValueError(f"every {label} time must be finite and positive")
    return resistive / alfven

def nearest_surface_spacing(psi_norm) -> np.ndarray:
    """Distance from each rational surface to its nearest neighbour.

    Parameters
    ----------
    psi_norm : array_like
        Rational-surface positions, in any order [-].

    Returns
    -------
    ndarray
        One distance per surface, in the input's order; ``nan`` for a single
        surface, which has no neighbour [-].

    Raises
    ------
    ValueError
        Two surfaces share a position.

    Processing steps
    ----------------
    1. Sort the surfaces outward.
    2. Take the smaller of the two gaps either side; an end surface has one.
    3. Put the result back in the caller's order.

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

    Provenance
    ----------
    .. [legacy] ``gpec_analysis_metrics._nearest_neighbour_spacing`` and
       ``gpec_multimode_metrics.compute_surface_chirikov``, which computed the
       same spacing twice and disagreed about repeated positions: one
       tolerated them, the other refused. This refuses.
    """
    ordered, _, order = _island_arrays(psi_norm, np.zeros_like(np.asarray(psi_norm, float)))
    result = np.full(ordered.size, np.nan)
    if ordered.size < 2:
        return result
    gaps = np.diff(ordered)
    result[0], result[-1] = gaps[0], gaps[-1]
    if ordered.size > 2:
        result[1:-1] = np.minimum(gaps[:-1], gaps[1:])
    unsorted = np.empty_like(result)
    unsorted[order] = result
    return unsorted

def penetration_ratio(psi_norm, width) -> np.ndarray:
    """Island width as a fraction of the width at which islands would overlap.

    Parameters
    ----------
    psi_norm : array_like
        Rational-surface positions [-].
    width : array_like
        Full island width per surface, in the same units [-].

    Returns
    -------
    ndarray
        ``width / critical_island_width``; one at overlap [-].

    Raises
    ------
    ValueError
        ``width`` is not one per surface, two surfaces share a position, or a
        position or width is not finite or a width is negative.

    Convention
    ----------
    A ratio, not a verdict: the value one is where islands touch, and
    comparing against it is the caller's to do.

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

    Provenance
    ----------
    .. [legacy] ``gpec_analysis_metrics.compute_critical_island_width``.
    """
    _island_arrays(psi_norm, width)  # refuses a tie or a shape mismatch
    with np.errstate(divide="ignore", invalid="ignore"):
        return np.asarray(width, dtype=float) / critical_island_width(psi_norm)

def rms_resonant_field(psi_norm, phi_res, *, window: ResonantWindow) -> float:
    """Root-mean-square resonant field over a radial window.

    Parameters
    ----------
    psi_norm : array_like
        Normalized poloidal flux of each rational surface [-].
    phi_res : array_like
        Resonant flux per surface, complex as GPEC writes it [T].
    window : ResonantWindow
        The interval to reduce over [-].

    Returns
    -------
    float
        RMS magnitude, in tesla; ``nan`` when the window holds no surface [T].

    Convention
    ----------
    Tesla, not gauss. GPEC normalizes this flux by the surface area, so it
    carries field units already; the code this replaces multiplied by 1e4
    inside the reduction, which put a unit conversion where nothing said so.

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

    Provenance
    ----------
    .. [legacy] ``gpec_analysis_metrics.compute_rms_resonant_field``, whose
       ``to_gauss=True`` default is the factor not applied here.
    """
    return reduce_resonant(psi_norm, phi_res, window=window, statistic="rms")

def shielded_field(delta, self_inductance) -> np.ndarray:
    """The field a finite-width layer's shielding current produces.

    Parameters
    ----------
    delta : array_like
        The derivative jump per surface, from :func:`finite_width_delta`,
        complex [T].
    self_inductance : array_like
        The surfaces' self-inductance matrix, square and matching ``delta``.
        **Required** [n/a].

    Returns
    -------
    ndarray
        ``L @ delta``, complex [T].

    Raises
    ------
    ValueError
        The matrix is not square, does not match ``delta``, or holds a
        non-finite value.

    Convention
    ----------
    The legacy returned ``delta`` itself when no inductance was supplied --
    a derivative jump relabelled as a shielded field, with nothing in the
    result to say the multiplication had not happened. There is no such
    fallback here: without the inductance there is no shielded field.

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

    Processing steps
    ----------------
    1. Check the matrix matches the surfaces.
    2. Multiply.

    Provenance
    ----------
    .. [legacy] ``resistive_gpec_metrics.py::shielded_field_from_delta``.
    """
    jump = np.atleast_1d(np.asarray(delta, dtype=complex))
    matrix = np.asarray(self_inductance, dtype=complex)
    if matrix.shape != (jump.size, jump.size):
        raise ValueError(
            f"the inductance is {matrix.shape} but there are {jump.size} surfaces; "
            "it must be square and one row per surface"
        )
    if not np.all(np.isfinite(matrix)) or not np.all(np.isfinite(jump)):
        raise ValueError("the inductance or the jump holds a non-finite value")
    return matrix @ jump