def boundary_curve(boundary: Boundary, sweep: str, values, swap_axes: bool = False, **fixed) -> BoundaryCurve:
    r"""A boundary sampled along one input, as $(x, y) = (x_\mathrm{sweep}, b)$ arrays.

    $$y_k = b\left(x_\mathrm{sweep} = v_k,\ \text{other inputs fixed}\right)$$

    Parameters
    ----------
    boundary : Boundary
        The boundary to sample [-].
    sweep : str
        Name of the input placed on the x axis [-].
    values : array_like
        Values of the swept input, in its declared unit [varies].
    swap_axes : bool
        Put the boundary's target on x and the swept input on y, as in a
        Hugill diagram, whose x axis is the density-like quantity [-].

    Returns
    -------
    BoundaryCurve
        ``x`` is the swept input and ``y`` the boundary value. The curve
        carries both quantities (identity and unit), the permitted side along
        y, and the fixed inputs [-].

    Raises
    ------
    KeyError
        ``sweep`` is not an input of the boundary.
    TypeError
        A fixed input is missing, unknown, an array, or duplicates the swept
        one; or ``values`` is not 1-D.

    Convention
    ----------
    The permitted side is the boundary's own ``allowed_side``, read along the
    y axis: ``"below"`` means operating points under the curve are permitted.
    With ``swap_axes`` the target is read along x, so ``"below"`` becomes
    ``"left"`` and ``"above"`` becomes ``"right"``.
    """
    x_quantity = boundary.input(sweep)
    if sweep in fixed:
        raise TypeError(f"{sweep!r} is swept; do not also fix it")
    non_scalar = [name for name, value in fixed.items() if np.ndim(value) != 0]
    if non_scalar:
        raise TypeError(f"fixed inputs must be scalars; {non_scalar} are arrays")
    x = np.atleast_1d(np.asarray(values, dtype=float))
    if x.ndim != 1:
        raise TypeError("values must be a 1-D sequence")
    y = np.broadcast_to(np.asarray(boundary_value(boundary, **{sweep: x}, **fixed), dtype=float), x.shape)
    fixed = {name: float(value) for name, value in fixed.items()}
    if swap_axes:
        side = {"below": "left", "above": "right"}[boundary.allowed_side]
        return BoundaryCurve(key=boundary.key, x=y, y=x, x_quantity=boundary.target, y_quantity=x_quantity,
                             allowed_side=side, fixed=fixed)
    return BoundaryCurve(key=boundary.key, x=x, y=y, x_quantity=x_quantity, y_quantity=boundary.target,
                         allowed_side=boundary.allowed_side, fixed=fixed)

def boundary_value(boundary: Boundary, **inputs):
    r"""Value of a boundary's target quantity at the given inputs.

    $$b = C\prod_i x_i^{\alpha_i}\quad\text{(power law)},\qquad b = C\quad\text{(threshold)},\qquad b = f(x_1,\ldots)\quad\text{(function)}$$

    Parameters
    ----------
    boundary : Boundary
        The boundary to evaluate [-].

    Returns
    -------
    float or np.ndarray
        Boundary value, in the unit of ``boundary.target`` [varies].

    Raises
    ------
    TypeError
        An input is missing or not one of the boundary's inputs.

    Convention
    ----------
    Each keyword input is a value in the unit its ``BoundaryQuantity``
    declares. No unit conversion happens here; a caller holding amperes for a
    boundary declared in MA converts first.

    Notes
    -----
    Array inputs broadcast with NumPy rules.
    """
    _check_inputs(boundary, inputs)
    if boundary.form == "threshold":
        return float(boundary.coefficient)
    if boundary.form == "function":
        return _scalar_or_array(boundary.function(**inputs))
    value = boundary.coefficient
    for name, exponent in boundary.exponents.items():
        value = value * np.asarray(inputs[name], dtype=float) ** exponent
    return _scalar_or_array(value)

def cylindrical_kink_coordinates(a, R0, B0, kappa, I_p):
    r"""Menard's cylindrical safety factor, $q^* = \pi a^2 B_{T0}(1+\kappa^2)/(\mu_0 R_0 I_P)$.

    $$q^* = \epsilon\,(1+\kappa^2)\,\frac{\pi a B_{T0}}{\mu_0 I_P}
          = \frac{\pi a^2 B_{T0}\,(1+\kappa^2)}{\mu_0 R_0 I_P}$$

    Parameters
    ----------
    a : float or np.ndarray
        Minor radius [m].
    R0 : float or np.ndarray
        Major radius [m].
    B0 : float or np.ndarray
        Vacuum toroidal field at ``R0``; its sign is dropped [T].
    kappa : float or np.ndarray
        Elongation [-].
    I_p : float or np.ndarray
        Plasma current; its sign is dropped [MA].

    Returns
    -------
    float or np.ndarray
        Cylindrical safety factor $q^*$, NaN where an input is NaN [-].

    Raises
    ------
    ValueError
        A finite non-positive (or infinite) minor radius, major radius, field
        magnitude or elongation.

    Convention
    ----------
    A global, shape-weighted proxy of the edge safety factor, not $q_a$ or $q_{95}$:
    Menard et al. use it because $q(1)$ and $q(0.95)$ at the current limit vary by a
    factor two with aspect ratio and shape while $q^*$ does not. It is the
    $(1+\kappa^2)/2$ form of Freidberg Eq. (13.171), and differs from Freidberg's
    Eq. (13.160) (``kink_coordinates``) by the factor $(1+\kappa^2)/(2\kappa)$.

    References
    ----------
    .. [1] J. E. Menard et al., Phys. Plasmas 11 (2004) 639; preprint PPPL-3908, p. 9.
    .. [2] D. A. D'Ippolito, J. P. Freidberg, J. P. Goedbloed and J. Rem, Phys. Fluids 21 (1978) 1600.
    """
    g = _positive_geometry(a=a, R0=R0, B0=np.abs(np.asarray(B0, dtype=float)), kappa=kappa)
    current = np.abs(np.asarray(I_p, dtype=float)) * 1e6
    with np.errstate(divide="ignore"):
        q = np.pi * g["a"] ** 2 * g["B0"] * (1.0 + g["kappa"] ** 2) / (MU0 * g["R0"] * current)
    return _scalar_or_array(q)

def evaluate_window(window: BoundaryWindow, operating_value, **inputs) -> WindowEvaluation:
    r"""Position of an operating value relative to an access window.

    $$b_\mathrm{lower} < x < b_\mathrm{upper}$$

    Parameters
    ----------
    window : BoundaryWindow
        The access window [-].
    operating_value : float or np.ndarray
        Operating value of the window's target quantity, in its unit [varies].

    Returns
    -------
    WindowEvaluation
        One evaluation per edge and ``inside``, true where both are
        permitted [-].

    Raises
    ------
    TypeError
        An input is missing or unused by both edges.

    Convention
    ----------
    ``inputs`` holds the union of both edges' inputs. Each edge receives only
    the inputs it declares.
    """
    known = set(window.lower.input_names) | set(window.upper.input_names)
    extra = [name for name in inputs if name not in known]
    if extra:
        raise TypeError(f"window {window.key!r} does not use inputs {extra}")
    lower = evaluate_boundary(
        window.lower, operating_value, **{k: inputs[k] for k in window.lower.input_names if k in inputs}
    )
    upper = evaluate_boundary(
        window.upper, operating_value, **{k: inputs[k] for k in window.upper.input_names if k in inputs}
    )
    inside = np.logical_and(lower.allowed, upper.allowed)
    return WindowEvaluation(
        key=window.key, lower=lower, upper=upper, inside=bool(inside) if np.ndim(inside) == 0 else inside
    )

def get_boundary(key: str):
    """A registered boundary or window by key.

    Parameters
    ----------
    key : str
        Registry key, for example ``"greenwald"`` [-].

    Returns
    -------
    Boundary or BoundaryWindow
        The registered entry [-].

    Raises
    ------
    KeyError
        No entry has that key; the message lists the keys that exist.
    """
    try:
        return _REGISTRY[key]
    except KeyError:
        raise KeyError(f"no boundary {key!r}; registered: {sorted(_REGISTRY)}") from None

def hugill_coordinates(n_e, R_geo, B_t, a, kappa_a, I_p):
    r"""Hugill-diagram coordinates of an operating state: $(\bar n_e R/B_T,\ 1/q_\mathrm{cyl})$.

    $$x = \frac{\bar n_e R_{geo}}{B_T},\qquad y = \frac{1}{q_{cyl}} = \frac{R_{geo}\,|I_p|}{5\,a^2\kappa_a B_T}\quad(I_p\ \mathrm{in\ MA})$$

    Parameters
    ----------
    n_e : float or np.ndarray
        Line-averaged electron density [1e19 m^-3].
    R_geo : float or np.ndarray
        Geometric major radius [m].
    B_t : float or np.ndarray
        Vacuum toroidal field at ``R_geo``, magnitude [T].
    a : float or np.ndarray
        Minor radius [m].
    kappa_a : float or np.ndarray
        Area elongation $S/(\pi a^2)$ [-].
    I_p : float or np.ndarray
        Plasma current; its sign is dropped [MA].

    Returns
    -------
    tuple of (float or np.ndarray)
        Murakami parameter $\bar n_e R/B_T$ [1e19 m^-2 T^-1] and $1/q_{cyl}$ [-].

    Raises
    ------
    ValueError
        A negative density, or a non-positive field, radius, minor radius or
        elongation.

    Convention
    ----------
    $q_{cyl}$ is ``equilibrium.q_cyl_from_B_R_epsilon_kappa_I`` with
    $\varepsilon = a/R_{geo}$. The ``"greenwald_hugill"`` boundary uses the
    same convention, so a state at $\bar n_e = n_G$ lies exactly on it. This
    is the cylindrical $q$, not $q_{95}$. At VEST aspect ratio the two differ
    substantially, and plotting $1/q_{95}$ against this boundary is only
    approximate.

    Numerical notes
    ---------------
    $1/q_{cyl}$ is formed directly rather than as ``1/q_cyl_from_...``, so a
    zero current maps to the origin $y = 0$ instead of raising. A whole time
    trace, including start-up and termination, can therefore be projected
    without cropping.
    """
    n = np.asarray(n_e, dtype=float)
    B_abs = np.abs(np.asarray(B_t, dtype=float))
    R = np.asarray(R_geo, dtype=float)
    a_arr = np.asarray(a, dtype=float)
    kappa = np.asarray(kappa_a, dtype=float)
    if np.any(n < 0):
        raise ValueError("n_e must be non-negative")
    for name, value in (("B_t", B_abs), ("R_geo", R), ("a", a_arr), ("kappa_a", kappa)):
        if np.any(~(value > 0)):
            raise ValueError(f"{name} must be positive and finite")
    x = n * R / B_abs
    y = np.abs(np.asarray(I_p, dtype=float)) * R / (5.0 * a_arr**2 * kappa * B_abs)
    return _scalar_or_array(x), _scalar_or_array(y)

def hugill_coordinates_st(n_e, R_geo, B_t, a, kappa, I_p):
    r"""Spherical-tokamak Hugill coordinates $(\bar n_e R/B_T,\ 1/q_\mathrm{cyl}^{ST})$ of Sykes et al. (2000).

    $$x = \frac{\bar n_e R}{B_T},\qquad
      y = \frac{1}{q^{ST}_{cyl}} = \frac{R\,|I_p|}{2.5\,a^2(1+\kappa^2)B_T}\quad(I_p\ \mathrm{in\ MA})$$

    Parameters
    ----------
    n_e : float or np.ndarray
        Line-averaged electron density [1e19 m^-3].
    R_geo : float or np.ndarray
        Major radius [m].
    B_t : float or np.ndarray
        Vacuum toroidal field at ``R_geo``; its sign is dropped [T].
    a : float or np.ndarray
        Minor radius [m].
    kappa : float or np.ndarray
        Elongation [-].
    I_p : float or np.ndarray
        Plasma current; its sign is dropped [MA].

    Returns
    -------
    tuple of (float or np.ndarray)
        Murakami parameter $\bar n_e R/B_T$ [1e19 m^-2 T^-1] and $1/q^{ST}_{cyl}$ [-].

    Raises
    ------
    ValueError
        A negative density, or a finite non-positive (or infinite) radius, minor
        radius, field magnitude or elongation.

    Convention
    ----------
    $q^{ST}_{cyl}$ is Menard's cylindrical $q^*$ (``cylindrical_kink_coordinates``),
    as Sykes et al. plot it for MAST; it differs from the conventional Hugill
    $q_{cyl}$ (``hugill_coordinates``) by $(1+\kappa^2)/(2\kappa_a)$, and the two are
    separate quantities. NaN inputs propagate.

    References
    ----------
    .. [1] A. Sykes et al., First results from MAST, IAEA FEC 2000, Sec. 6, Fig. 12;
           Nucl. Fusion 41 (2001) 1423.
    """
    g = _positive_geometry(R_geo=R_geo, B_t=np.abs(np.asarray(B_t, dtype=float)), a=a, kappa=kappa)
    n = np.asarray(n_e, dtype=float)
    if np.any(np.isfinite(n) & (n < 0)):
        raise ValueError("n_e must be non-negative where it is given")
    x = n * g["R_geo"] / g["B_t"]
    y = g["R_geo"] * np.abs(np.asarray(I_p, dtype=float)) / (2.5 * g["a"] ** 2 * (1.0 + g["kappa"] ** 2) * g["B_t"])
    return _scalar_or_array(x), _scalar_or_array(y)

def iter_q95_coordinates(a, R0, B0, kappa, delta, I_p):
    r"""The ITER design-guideline $q_{95}$ estimate from global shape, $q_{95} \approx q_* f(\epsilon)$.

    $$q_{95} \approx \frac{5a^2B}{R\,I_p[\mathrm{MA}]}\,
      \frac{1+\kappa^2(1+2\delta^2-1.2\delta^3)}{2}\,\frac{1.17-0.65\epsilon}{(1-\epsilon^2)^2},
      \qquad \epsilon = a/R$$

    Parameters
    ----------
    a : float or np.ndarray
        Minor radius [m].
    R0 : float or np.ndarray
        Major radius [m].
    B0 : float or np.ndarray
        Vacuum toroidal field at ``R0``; its sign is dropped [T].
    kappa : float or np.ndarray
        Elongation of the 95 % flux surface, kappa_95 [-].
    delta : float or np.ndarray
        Triangularity of the 95 % flux surface, delta_95; may be zero or negative [-].
    I_p : float or np.ndarray
        Plasma current; its sign is dropped [MA].

    Returns
    -------
    float or np.ndarray
        Estimated $q_{95}$, NaN where an input is NaN [-].

    Raises
    ------
    ValueError
        A finite non-positive (or infinite) minor radius, major radius, field
        magnitude or elongation, or $a \ge R_0$.

    Convention
    ----------
    The formula is calibrated on the 95 % flux-surface shape ($\kappa_{95}$,
    $\delta_{95}$); fed the last-closed-surface $\kappa$ and $\delta$, which are larger,
    it over-estimates $q_{95}$ (for the ITER design point by about 27 %).
    A fit for conventional aspect ratio used for ITER design; at a spherical
    tokamak's $A \approx 1.3$ it is an extrapolation and over-estimates $q_{95}$
    (Akers et al. 2000 give a START-based correction, not registered here). It is an
    estimate from global parameters, not an equilibrium $q_{95}$.

    References
    ----------
    .. [1] D. E. Post et al., *ITER Physics*, ITER Documentation Series No. 21, IAEA (1991),
           Table 1-2.
    .. [2] N. A. Uckan and ITER Physics Group, *ITER Physics Design Guidelines: 1989*,
           IAEA/ITER/DS-10, IAEA (1990), p. 10.
    """
    g = _positive_geometry(a=a, R0=R0, B0=np.abs(np.asarray(B0, dtype=float)), kappa=kappa)
    if np.any(np.isfinite(g["a"] / g["R0"]) & (g["a"] / g["R0"] >= 1.0)):
        raise ValueError("a must be smaller than R0")
    d = np.asarray(delta, dtype=float)
    if np.any(np.isinf(d)):
        raise ValueError("delta must be finite where it is given")
    current = np.abs(np.asarray(I_p, dtype=float))
    with np.errstate(divide="ignore"):
        q = _iter_q95_per_ma(g["a"], g["R0"], g["B0"], g["kappa"], d) / current
    return _scalar_or_array(q)

def kink_coordinates(a, R0, B0, kappa, I_p):
    r"""Freidberg's kink safety factor of an elongated tokamak, $q_* = 2\pi a^2\kappa B_0/(\mu_0 R_0 I)$.

    $$q_* = \frac{2\pi a^2 \kappa B_0}{\mu_0 R_0 I_p} = \frac{5\,a^2\kappa B_0}{R_0\,I_p[\mathrm{MA}]}$$

    Parameters
    ----------
    a : float or np.ndarray
        Minor radius [m].
    R0 : float or np.ndarray
        Major radius [m].
    B0 : float or np.ndarray
        Vacuum toroidal field at ``R0``; its sign is dropped [T].
    kappa : float or np.ndarray
        Elongation [-].
    I_p : float or np.ndarray
        Plasma current; its sign is dropped [MA].

    Returns
    -------
    float or np.ndarray
        Kink safety factor $q_*$ of Eq. (13.160), NaN where an input is NaN [-].

    Raises
    ------
    ValueError
        A finite non-positive (or infinite) minor radius, major radius, field
        magnitude or elongation.

    Convention
    ----------
    Eq. (13.160) is the $G(\kappa) \approx 1$ form of Eq. (13.158), the definition with
    which the kink limit ``"freidberg_2008_kink_qstar"`` (Eq. 13.162) is stated. It is
    not ``vaft.formula.equilibrium.kink_safety_factor`` (#1524). A zero current maps to
    an infinite $q_*$.

    References
    ----------
    .. [1] J. P. Freidberg, *Plasma Physics and Fusion Energy*, Cambridge University
           Press (2008), Eqs. (13.158)-(13.160), p. 405.
    """
    a_arr, R, B_abs, k = (np.asarray(v, dtype=float) for v in (a, R0, np.abs(np.asarray(B0, dtype=float)), kappa))
    for name, value in (("a", a_arr), ("R0", R), ("B0", B_abs), ("kappa", k)):
        # NaN propagates (a population table has gaps); a finite non-positive value or an infinity is an error
        if np.any(np.isinf(value)) or np.any(np.isfinite(value) & ~(value > 0)):
            raise ValueError(f"{name} must be positive and finite where it is given")
    current = np.abs(np.asarray(I_p, dtype=float)) * 1e6
    with np.errstate(divide="ignore"):
        q = 2.0 * np.pi * a_arr**2 * k * B_abs / (MU0 * R * current)
    return _scalar_or_array(q)

def list_boundaries(family: Optional[str] = None) -> tuple:
    """Keys of the registered boundaries, optionally restricted to one family.

    Parameters
    ----------
    family : str, optional
        Family name, for example ``"density_limit"`` [-].

    Returns
    -------
    tuple of str
        Sorted registry keys [-].
    """
    return tuple(sorted(
        key for key, entry in _REGISTRY.items()
        if family is None or getattr(entry, "family", getattr(getattr(entry, "lower", None), "family", None)) == family
    ))

def same_quantity(a: BoundaryQuantity, b: BoundaryQuantity) -> bool:
    r"""Whether two quantities are the same physical quantity in the same unit.

    Parameters
    ----------
    a, b : BoundaryQuantity
        The quantities to compare [-].

    Returns
    -------
    bool
        ``True`` when both ``name`` (the identity) and ``unit`` match. The
        symbol and the free-text definition are not compared [-].

    Convention
    ----------
    Identity is the ``name``: ``edge_safety_factor`` ($q_\psi$),
    ``edge_safety_factor_95`` ($q_{95}$) and ``inverse_cylindrical_q`` are
    different quantities even where their values are close. A boundary may be
    drawn on an axis only when this returns ``True`` for that axis (#1425).
    """
    return a.name == b.name and a.unit == b.unit

def start_q95_coordinates(a, R0, B0, kappa, delta, I_p, configuration="limiter"):
    r"""The START low-aspect-ratio $q_{95}$ estimate of Akers et al. (2000) from global shape.

    $$q_{95} \approx \frac{5a^2B}{R\,I_p[\mathrm{MA}]}\,
      \frac{1+\kappa^2(1+2\delta^2-1.2\delta^3)}{2}\;1.17\,C\sqrt{\frac{A}{A-1}},
      \qquad A = R/a$$

    Parameters
    ----------
    a : float or np.ndarray
        Minor radius [m].
    R0 : float or np.ndarray
        Major radius [m].
    B0 : float or np.ndarray
        Vacuum toroidal field at ``R0``; its sign is dropped [T].
    kappa : float or np.ndarray
        Elongation of the 95 % flux surface, kappa_95 [-].
    delta : float or np.ndarray
        Triangularity of the 95 % flux surface, delta_95 [-].
    I_p : float or np.ndarray
        Plasma current; its sign is dropped [MA].
    configuration : {"limiter", "double_null"}
        Sets C = 1.0 (natural limiter plasma) or 0.77 (double null) [-].

    Returns
    -------
    float or np.ndarray
        Estimated $q_{95}$, NaN where an input is NaN [-].

    Raises
    ------
    ValueError
        A finite non-positive (or infinite) minor radius, major radius, field
        magnitude or elongation, an infinite triangularity, $a \ge R_0$, or an
        unknown configuration.

    Convention
    ----------
    A fit to START equilibria (spherical tokamak), more conservative in aspect ratio
    than the ITER guideline formula (``iter_q95_coordinates``), whose shaping factor and
    normalisation it shares. Like it, it is written for the 95 % surface shape.

    References
    ----------
    .. [1] R. J. Akers et al., Nucl. Fusion 40 (2000) 1223, Sec. 2.1, p. 1227.
    """
    if configuration not in AKERS_2000_C:
        raise ValueError(f"configuration must be one of {sorted(AKERS_2000_C)}, not {configuration!r}")
    g = _positive_geometry(a=a, R0=R0, B0=np.abs(np.asarray(B0, dtype=float)), kappa=kappa)
    if np.any(np.isfinite(g["a"] / g["R0"]) & (g["a"] / g["R0"] >= 1.0)):
        raise ValueError("a must be smaller than R0")
    d = np.asarray(delta, dtype=float)
    if np.any(np.isinf(d)):
        raise ValueError("delta must be finite where it is given")
    current = np.abs(np.asarray(I_p, dtype=float))
    with np.errstate(divide="ignore"):
        q = _start_q95_per_ma(g["a"], g["R0"], g["B0"], g["kappa"], d, AKERS_2000_C[configuration]) / current
    return _scalar_or_array(q)

def threshold_curve(boundary: Boundary, axis: BoundaryQuantity, values, target_axis: str = "y") -> BoundaryCurve:
    r"""A threshold boundary drawn as a straight line across a 2-D projection.

    $$\text{target} = C \quad\text{for every value of the other axis}$$

    Parameters
    ----------
    boundary : Boundary
        A ``"threshold"`` boundary, for example ``"murakami_hugill"`` [-].
    axis : BoundaryQuantity
        The quantity on the other axis, which the threshold does not depend on [-].
    values : array_like
        Sample points along that axis, in its declared unit [varies].
    target_axis : str
        ``"y"`` draws a horizontal line (target on y); ``"x"`` a vertical one [-].

    Returns
    -------
    BoundaryCurve
        The line, carrying both quantities and the permitted side: ``"below"``
        / ``"above"`` for a horizontal line, ``"left"`` / ``"right"`` for a
        vertical one [-].

    Raises
    ------
    TypeError
        ``boundary`` is not a threshold, or ``values`` is not 1-D.
    ValueError
        ``target_axis`` is neither ``"x"`` nor ``"y"``, or ``axis`` is the
        boundary's own target.

    Convention
    ----------
    :func:`boundary_curve` sweeps one of a boundary's inputs, and a threshold
    has none, so this is its counterpart. The value is the registered
    coefficient; nothing is restated.
    """
    if boundary.form != "threshold":
        raise TypeError(f"boundary {boundary.key!r} is a {boundary.form!r}, not a threshold; use boundary_curve")
    if target_axis not in ("x", "y"):
        raise ValueError(f"target_axis must be 'x' or 'y', not {target_axis!r}")
    if same_quantity(axis, boundary.target):
        raise ValueError(f"the other axis cannot be the threshold's own target {boundary.target.name!r}")
    v = np.atleast_1d(np.asarray(values, dtype=float))
    if v.ndim != 1:
        raise TypeError("values must be a 1-D sequence")
    level = np.full(v.shape, float(boundary_value(boundary)))
    if target_axis == "x":
        side = {"below": "left", "above": "right"}[boundary.allowed_side]
        return BoundaryCurve(key=boundary.key, x=level, y=v, x_quantity=boundary.target, y_quantity=axis,
                             allowed_side=side)
    return BoundaryCurve(key=boundary.key, x=v, y=level, x_quantity=axis, y_quantity=boundary.target,
                         allowed_side=boundary.allowed_side)