def connection_length_weight(length, *, scale: float, form: str) -> np.ndarray:
    """Parallel-loss weight from a field line's connection length.

    Parameters
    ----------
    length : array_like
        Connection length per line; negative entries are treated as zero,
        and non-finite ones are given zero weight [m].
    scale : float
        The length the weight is measured against; must be positive [m].
    form : str
        ``"inverse"`` for ``1 / (1 + L / L_0)``, ``"exponential"`` for
        ``exp(-L / L_0)``, ``"none"`` for a flat one [n/a].

    Returns
    -------
    ndarray
        Weights in ``[0, 1]`` [-].

    Raises
    ------
    ValueError
        ``scale`` is not positive, or ``form`` is not one of
        :data:`CONNECTION_LENGTH_FORMS`.

    Convention
    ----------
    A longer line is weighted less: the length stands in for how much of the
    flux tube's content was lost on the way. The two forms are empirical and
    disagree by construction -- the inverse one keeps a long line at a small
    but finite weight, the exponential one takes it to zero -- so which is
    used is part of the model, never a default.

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

    Provenance
    ----------
    .. [legacy] the donor's ``connection_length_weight``, whose
       ``"exp_loss"`` is this ``"exponential"``.
    """
    if form not in CONNECTION_LENGTH_FORMS:
        raise ValueError(
            f"form must be one of {list(CONNECTION_LENGTH_FORMS)}, not {form!r}"
        )
    values = np.asarray(length, dtype=float)
    if form == "none":
        return np.ones(values.shape)
    if not scale > 0.0:
        raise ValueError(f"scale must be positive, not {scale!r}")
    finite = np.isfinite(values)
    positive = np.maximum(np.where(finite, values, 0.0), 0.0)
    if form == "inverse":
        weight = 1.0 / (1.0 + positive / float(scale))
    else:
        weight = np.exp(-positive / float(scale))
    return np.where(finite, weight, 0.0)

def footprint_incident_total(proxy: FootprintProxy) -> float:
    """The proxy integrated over one target.

    Parameters
    ----------
    proxy : FootprintProxy
        As :func:`footprint_heat_load_proxy` returns it [-].

    Returns
    -------
    float
        ``sum(density * area)``, in the proxy's own relative scale [-].

    Processing steps
    ----------------
    Non-finite samples are skipped rather than propagated: a line the tracer
    could not follow contributes nothing, and one NaN would otherwise erase
    the target.

    Convention
    ----------
    This is a Riemann sum over the sample areas, so it converges as the grid
    is refined -- which is the property a fraction between targets needs, and
    the one the legacy's extra ``1/area`` factor destroyed. With it the sum
    reduces to ``sum(w)``, which grows roughly with the number of samples the
    caller chose to trace rather than settling on the target's own value.

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

    Provenance
    ----------
    .. [legacy] the donor's ``summarize_heat_flux_proxy``, whose
       ``integral`` column is this same sum -- taken over a density that
       carried an extra ``1/area`` factor, which is what reduced it to the
       bare sum of the weights.
    """
    values = proxy.density * proxy.area
    return float(np.nansum(np.where(np.isfinite(values), values, 0.0)))

def footprint_proxy_model(name: str) -> FootprintProxyModel:
    """One named footprint proxy model.

    Parameters
    ----------
    name : str
        A key of :data:`FOOTPRINT_PROXY_MODELS` [n/a].

    Returns
    -------
    FootprintProxyModel
        The preset, which is frozen and safe to share [n/a].

    Raises
    ------
    KeyError
        No model carries that name; the message lists the ones that do.

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

    Provenance
    ----------
    .. [legacy] the donor's ``FootprintProxyConfig``, whose
       field defaults these presets carry as named values.
    """
    try:
        return FOOTPRINT_PROXY_MODELS[name]
    except KeyError:
        raise KeyError(
            f"unknown footprint proxy model {name!r}; "
            f"registered: {sorted(FOOTPRINT_PROXY_MODELS)}"
        ) from None

def incidence_factor(angle_deg) -> np.ndarray:
    """The projection of a field line onto the surface it strikes.

    Parameters
    ----------
    angle_deg : array_like
        Angle between the field line and the target surface, as FLARE's
        ``alphaS`` reports it [deg].

    Returns
    -------
    ndarray
        ``|sin alpha|``, in ``[0, 1]`` [-].

    Convention
    ----------
    The angle is measured from the *surface*, not from its normal, so a
    grazing line has a small angle and a small factor. That is FLARE's own
    definition: ``alphaS = asin(|b.n| / |b|)``, which makes ``|sin alpha|``
    the projection exactly rather than an approximation of it. A line that
    never struck the boundary is written as zero by the tracer and so
    receives zero weight, which is the right answer for a different reason --
    it deposits nothing because it does not arrive.

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

    Provenance
    ----------
    .. [1] ``FLARE/src/fortran/tasks/fieldline_connection.f90:158``.
    .. [legacy] the donor's ``incidence_factor``, which returned
       ones when the column was absent; see
       :func:`footprint_heat_load_proxy`.
    """
    alpha = np.deg2rad(np.abs(np.asarray(angle_deg, dtype=float)))
    return np.clip(np.abs(np.sin(alpha)), 0.0, 1.0)

def reduce_traced_directions(backward, forward, *, how: str) -> np.ndarray:
    """Reduce a quantity traced in both directions to one value per line.

    Parameters
    ----------
    backward : array_like
        The value the tracer reported going one way [any].
    forward : array_like
        The value it reported going the other [any].
    how : str
        ``"shortest"`` takes the smaller magnitude, ``"total"`` the sum,
        ``"deepest"`` the smaller signed value [n/a].

    Returns
    -------
    ndarray
        One value per line, in the inputs' own unit [any].

    Raises
    ------
    ValueError
        The arrays disagree in shape, or ``how`` is none of the three.

    Convention
    ----------
    ``"shortest"`` compares magnitudes, because a tracer may report a
    backward length as negative; ``"deepest"`` compares signed values,
    because a flux label that ran further inward is genuinely smaller.
    Non-finite entries propagate rather than being skipped: a direction the
    tracer could not follow leaves the reduction undefined, and that is the
    honest answer for that line.

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

    Provenance
    ----------
    .. [1] ``FLARE/src/fortran/tasks/fieldline_connection.f90:90-96``, whose
       own ``Lc = Lc_bwd + Lc_fwd``, ``Lcs = min(Lc_bwd, Lc_fwd)`` and
       ``minPsiN = min(minPsiN_bwd, minPsiN_fwd)`` are these three
       reductions. FLARE's arclengths are non-negative, so its plain ``min``
       and the magnitude comparison here agree on its own output.
    """
    first = np.asarray(backward, dtype=float)
    second = np.asarray(forward, dtype=float)
    if first.shape != second.shape:
        raise ValueError(
            f"the two directions must share a shape; got {first.shape} and "
            f"{second.shape}"
        )
    if how == "shortest":
        return np.where(np.abs(first) <= np.abs(second), first, second)
    if how == "total":
        return first + second
    if how == "deepest":
        return np.minimum(first, second)
    raise ValueError(
        f'how must be "shortest", "total" or "deepest", not {how!r}'
    )

def target_incident_fractions(totals: Mapping[str, float]) -> dict[str, float]:
    """Each target's share of the proxy incident on all of them.

    Parameters
    ----------
    totals : mapping of str to float
        Per-target integrals, from :func:`footprint_incident_total`; each
        must be finite and not negative [-].

    Returns
    -------
    dict of str to float
        The same keys, each divided by the sum over all of them, so the
        values add to one [-].

    Raises
    ------
    ValueError
        ``totals`` is empty, an entry is negative or not finite, or they sum
        to zero -- in which case there is no share to take, and returning
        equal shares would invent one.

    Convention
    ----------
    The denominator is the total over the targets *in this mapping*, so the
    fraction says what it says only if that set is every target the run
    covered. A fraction over a subset is a fraction of the subset.

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

    Provenance
    ----------
    .. [owner] The IMAS stand-in decided on vaft#1099 (2026-09-21): the proxy
       integrated over each target, divided by the total over all targets,
       written to ``divertors.divertor[].target[].power_incident_fraction``
       because it is dimensionless and so is the proxy.
    """
    if not totals:
        raise ValueError("no targets to take a share of")
    for key, value in totals.items():
        if not np.isfinite(value) or value < 0.0:
            raise ValueError(
                f"target {key!r} has an incident total of {value!r}; a share "
                "is defined only for finite, non-negative totals"
            )
    total = float(sum(float(value) for value in totals.values()))
    if total <= 0.0:
        raise ValueError(
            "every target's incident total is zero, so there is no share to "
            "take; equal shares would be an invention, not a result"
        )
    return {key: float(value) / total for key, value in totals.items()}

def toroidal_surface_cell_areas(r, z, phi_deg) -> np.ndarray:
    """Area of every quadrilateral cell of a surface mesh in ``(R, z, phi)``.

    Parameters
    ----------
    r : array_like
        Node major radius, shape ``(n_v, n_u)`` [m].
    z : array_like
        Node height, same shape [m].
    phi_deg : array_like
        Node toroidal angle, same shape [deg].

    Returns
    -------
    ndarray
        Cell areas, shape ``(n_v - 1, n_u - 1)`` [m^2].

    Raises
    ------
    ValueError
        The three arrays disagree in shape, or the mesh is not two
        dimensional with at least two nodes on each axis.

    Processing steps
    ----------------
    The nodes are placed in Cartesian space, and each cell is split along the
    diagonal from ``[j, i+1]`` to ``[j+1, i]`` into two triangles whose areas
    are summed. No small-angle expansion is taken, so the result is exact for
    a flat quadrilateral. A cell that is not flat has two diagonals that give
    slightly different answers, and this is the one FLARE splits along --
    which is why the two agree to parts in a trillion rather than merely
    closely.

    Convention
    ----------
    Node order is ``[i_v, i_u]``, the order FLARE's own mesh files lay their
    values out in -- see :attr:`vaft.data.flare_products.FlareMesh.node_shape`.
    Cell ``[j, i]`` is bounded by nodes ``[j, i]``, ``[j, i+1]``,
    ``[j+1, i+1]`` and ``[j+1, i]``.

    Applicability
    -------------
    Machine-independent. Any surface described by node positions in
    cylindrical coordinates; nothing here is specific to a divertor target.

    Provenance
    ----------
    .. [1] ``moose/src/fortran/geometry/hypermesh3d.f90:563``, the same
       two-triangle construction FLARE uses for the ``area`` column of its
       strike-density product.
    """
    x = _cartesian(r, z, phi_deg)
    corner_00, corner_01 = x[:-1, :-1], x[:-1, 1:]
    corner_11, corner_10 = x[1:, 1:], x[1:, :-1]
    first = np.linalg.norm(np.cross(corner_01 - corner_00, corner_10 - corner_00), axis=-1)
    second = np.linalg.norm(np.cross(corner_01 - corner_11, corner_10 - corner_11), axis=-1)
    return (first + second) / 2.0

def toroidal_surface_node_areas(r, z, phi_deg) -> np.ndarray:
    """Area each node of a surface mesh represents, as a dual cell.

    Parameters
    ----------
    r : array_like
        Node major radius, shape ``(n_v, n_u)`` [m].
    z : array_like
        Node height, same shape [m].
    phi_deg : array_like
        Node toroidal angle, same shape [deg].

    Returns
    -------
    ndarray
        Node areas, shape ``(n_v, n_u)``, summing to the total surface area
        [m^2].

    Raises
    ------
    ValueError
        As :func:`toroidal_surface_cell_areas`.

    Processing steps
    ----------------
    Each cell's area is split equally between its four corners. A corner node
    therefore receives a quarter of one cell, an edge node a quarter of two,
    and an interior node a quarter of four -- so the node areas sum to the
    surface area exactly, which is what makes an integral over them the same
    integral whether the values sit on nodes or on cells.

    Convention
    ----------
    Node order is ``[i_v, i_u]``, as :func:`toroidal_surface_cell_areas`.

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

    Provenance
    ----------
    .. [legacy] The donor's ``footprint_cell_areas`` built the
       same quantity from midpoint axis widths and ``R dphi ds``, which
       assumes the second axis is parameterized by its own arc length. This
       form reads the arc from the node positions instead, so a target curve
       whose parameter is not arc length is measured rather than assumed.
    """
    cells = toroidal_surface_cell_areas(r, z, phi_deg)
    quarter = cells / 4.0
    nodes = np.zeros(np.asarray(r, dtype=float).shape)
    nodes[:-1, :-1] += quarter
    nodes[:-1, 1:] += quarter
    nodes[1:, :-1] += quarter
    nodes[1:, 1:] += quarter
    return nodes

def upstream_flux_weight(psi_norm, *, separatrix_flux: float, decay_width: float,
                         clip_inside: bool = True) -> np.ndarray:
    """Exponential upstream heat-channel weight in normalized flux.

    Parameters
    ----------
    psi_norm : array_like
        The flux label each field line is weighted by -- for a connection
        product, the deepest normalized flux it reached [-].
    separatrix_flux : float
        Where the channel starts, normally 1 [-].
    decay_width : float
        e-folding width of the channel, in the same normalized flux; must be
        positive [-].
    clip_inside : bool, optional
        Whether a line that reached inside ``separatrix_flux`` is given the
        full weight rather than one growing exponentially. ``True`` clips
        [n/a].

    Returns
    -------
    ndarray
        Weights in ``(0, 1]`` when clipped, unbounded above otherwise [-].

    Raises
    ------
    ValueError
        ``decay_width`` is not positive.

    Convention
    ----------
    ``psi_norm`` increases outward, so a line that stayed further out is
    weighted less. With ``clip_inside`` false the weight grows without bound
    inside the separatrix, which is the model being read outside the region
    it describes; it is offered because a caller comparing against the
    unclipped legacy form needs it, not because it is ever the better choice.

    Applicability
    -------------
    Machine-independent. The channel is an upstream scrape-off-layer width,
    so the weight is meaningful only for lines that reach the plasma edge.

    Provenance
    ----------
    .. [legacy] the donor's ``upstream_psi_weight``.
    """
    if not decay_width > 0.0:
        raise ValueError(f"decay_width must be positive, not {decay_width!r}")
    delta = np.asarray(psi_norm, dtype=float) - float(separatrix_flux)
    if clip_inside:
        delta = np.maximum(delta, 0.0)
    return np.exp(-delta / float(decay_width))