def compute_br_bz_phi(
    r_obs: np.ndarray,
    z_obs: np.ndarray,
    r_src: float,
    z_src: float,
    shift: float = 0.01
    ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Field and flux of one unit-current ring at a set of observation points.

    Parameters
    ----------
    r_obs : array_like
        Major radius of the observation points [m].
    z_obs : array_like
        Height of the observation points [m].
    r_src : float
        Major radius of the source ring [m].
    z_src : float
        Height of the source ring [m].
    shift : float, optional
        Offset used to step around the source singularity [m].

    Returns
    -------
    tuple of np.ndarray
        The radial and vertical field per unit current in tesla per ampere, and
        the poloidal flux per unit current in weber per ampere [-].

    Convention
    ----------
    Per unit current, so the caller multiplies by an actual current. An
    observation point sitting on the source has no finite response; rather than
    excluding it, the response is evaluated at two points either side and
    averaged, which keeps the returned array the same shape as the input and gives
    a finite value whose error is second order in the offset.

    Defaults
    --------
    ``shift = 0.01`` is a numerical convenience. The averaging is applied only
    within a third of it from the source, so a point comfortably away from the
    ring is untouched.

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

    Limitations
    -----------
    Free space: no iron, no shielding, no image currents. The source is an ideal
    filament, so the response very close to a real conductor of finite section is
    not physical, which is what the offset is papering over rather than solving.

    Provenance
    ----------
    .. [1] The axisymmetric Green's functions in :mod:`vaft.formula.green`.
    """
    distances = calculate_distance(r_obs, r_src, z_obs, z_src)

    condition = distances < (shift / 3.0)

    br1_shifted, bz1_shifted = green_br_bz(r_obs + shift, z_obs, r_src, z_src)
    br2_shifted, bz2_shifted = green_br_bz(r_obs - shift, z_obs, r_src, z_src)
    phi1_shifted = green_r(r_obs + shift, z_obs, r_src, z_src)
    phi2_shifted = green_r(r_obs - shift, z_obs, r_src, z_src)

    br_direct, bz_direct = green_br_bz(r_obs, z_obs, r_src, z_src)
    phi_direct = green_r(r_obs, z_obs, r_src, z_src)

    br_final = np.where(condition, (br1_shifted + br2_shifted) / 2.0, br_direct)
    bz_final = np.where(condition, (bz1_shifted + bz2_shifted) / 2.0, bz_direct)
    phi_final = np.where(condition, (phi1_shifted + phi2_shifted) / 2.0, phi_direct)

    return br_final, bz_final, phi_final

def compute_mutual_passive_active(
    passive_loop_geometry: List[Tuple[str, float, float, float]],
    coil_geometry: List[List[Tuple[float, float, int]]],
) -> np.ndarray:
    """Mutual inductance between passive loops and active coils, derived from geometry.

    Parameters
    ----------
    passive_loop_geometry : sequence
        Geometry of each passive loop [-].
    coil_geometry : sequence
        Geometry and turns of each active coil [-].

    Returns
    -------
    np.ndarray
        Mutual inductance with one row per loop and one column per coil [H].

    Convention
    ----------
    Derived from the geometry given, which is what makes it a usable fallback when
    a packaged coupling matrix does not apply. The packaged matrix represents a
    historical coil geometry; a machine whose coils have since moved needs this
    instead, and the two must not be mixed within one calculation.

    Applicability
    -------------
    Machine-independent.  The geometry is the argument; the packaged matrix it
    substitutes for is VEST's.

    Limitations
    -----------
    Filament idealization, so it is less accurate than a matrix computed from the
    real conductor sections. Use it when the packaged matrix does not match the
    geometry, not in preference to one that does.

    Provenance
    ----------
    .. [1] The Green's functions in :mod:`vaft.formula.green`; the packaged
       alternative is the machine layer's ``em_coupling`` asset.
    """
    coupling = np.zeros((len(passive_loop_geometry), len(coil_geometry)))

    for i_loop, (_, r_loop, z_loop, coefficient) in enumerate(
        passive_loop_geometry
    ):
        for i_coil, elements in enumerate(coil_geometry):
            if not elements:
                raise ValueError(f"active coil {i_coil} has no geometry elements")
            total_turns = sum(element[2] for element in elements)
            mean_response = sum(
                green_r(r_loop, z_loop, r_coil, z_coil)
                for r_coil, z_coil, _ in elements
            ) / len(elements)
            coupling[i_loop, i_coil] = coefficient * total_turns * mean_response

    return coupling

def compute_response_vector(
    coil_data: List[Dict[str, Any]],
    passive_loop_data: List[Dict[str, Any]],
    plasma_points: List[List[float]],
    observation_points: List[List[float]]
    ) -> Tuple[ndarray, ndarray, ndarray]:
    """The same response as the matrix form, with the arguments in the legacy order.

    Parameters
    ----------
    coil_data : sequence
        Active coil geometry and turns [-].
    passive_loop_data : sequence
        Passive loop geometry [-].
    plasma_points : array_like
        Plasma filament positions [m].
    observation_points : array_like
        Points at which to evaluate [m].

    Returns
    -------
    tuple of np.ndarray
        Exactly what :func:`compute_response_matrix` returns [-].

    Convention
    ----------
    A thin wrapper that reorders its arguments and delegates. It exists because
    the legacy call order put the observation points last; new code should call
    the matrix form directly.

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

    Provenance
    ----------
    .. [1] :func:`compute_response_matrix`, which does the work.
    """
    return compute_response_matrix(
        observation_points=observation_points,
        coil_data=coil_data,
        passive_loop_data=passive_loop_data,
        plasma_points=plasma_points
    )

def compute_vacuum_fields_1d(
    coil_plus_loop_currents: np.ndarray,  # shape (n_times, nb_coil+nb_loop)
    coil_plus_loop_psi_resp: np.ndarray,  # shape (n_points, nb_coil+nb_loop)
    coil_plus_loop_br_resp: np.ndarray,   # shape (n_points, nb_coil+nb_loop)
    coil_plus_loop_bz_resp: np.ndarray,   # shape (n_points, nb_coil+nb_loop)
    ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Contract current histories with response matrices into field histories.

    Parameters
    ----------
    coil_plus_loop_currents : array_like
        Current history, one row per time and one column per source [A].
    coil_plus_loop_psi_resp : array_like
        Flux response per unit current [Wb/A].
    coil_plus_loop_br_resp : array_like
        Radial field response per unit current [T/A].
    coil_plus_loop_bz_resp : array_like
        Vertical field response per unit current [T/A].

    Returns
    -------
    tuple of np.ndarray
        Flux in weber and the two field components in tesla, each with one row per
        time and one column per observation point [-].

    Convention
    ----------
    The linear step the response matrices exist for: each output is the current
    history contracted against its own response. **Vacuum fields only** -- the
    plasma's own contribution is not here unless it entered as a filament current
    in the input.

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

    Limitations
    -----------
    Column order in the currents must match the response matrices' own, and
    nothing checks it; a mismatched ordering gives plausible numbers that are
    wrong.

    Provenance
    ----------
    .. [1] The response matrices built by :func:`compute_response_matrix` or
       :func:`calc_grid`.
    """
    n_times = coil_plus_loop_currents.shape[0]
    n_points = coil_plus_loop_psi_resp.shape[0]

    psi_out = np.zeros((n_times, n_points))
    br_out = np.zeros((n_times, n_points))
    bz_out = np.zeros((n_times, n_points))

    for i_time in range(n_times):
        ix = coil_plus_loop_currents[i_time]  # shape (nb_coil+nb_loop,)
        psi_out[i_time] = coil_plus_loop_psi_resp @ ix
        br_out[i_time]  = coil_plus_loop_br_resp  @ ix
        bz_out[i_time]  = coil_plus_loop_bz_resp  @ ix

    return psi_out, br_out, bz_out