Skip to content

Components

device_inductance.coils

Coil dataclass

Coil(
    name: str,
    resistance: float,
    self_inductance: float,
    filaments: list[CoilFilament],
)

An axisymmetric magnet, which may not have a rectangular cross-section

extent cached property

extent: tuple[float, float, float, float]

[m] rmin, rmax, zmin, zmax extent of filament centers

filaments instance-attribute

filaments: list[CoilFilament]

Discretized circular filaments describing the coil's winding pattern

grids cached property

grids: tuple[NDArray, NDArray]

Generate a grid that spans the winding pack, landing exactly on the windings if possible.

local_field_interpolators cached property

local_field_interpolators: tuple[
    MulticubicRegular, MulticubicRegular, MulticubicRegular
]

Build interpolators over the solved local fields, if available

local_fields cached property

local_fields: LocalFields

Solve the local self-field flux and flux density per amp by mapping the coil section to a continuous current density distribution and solving the continuous flux field via 4th-order finite difference, then extracting the B-field from the flux field via 4th-order finite difference.

Returns:

Type Description
LocalFields

Structure containing local field map components and interpolators

meshes cached property

meshes: tuple[NDArray, NDArray]

Generate a meshgrid that spans the winding pack, landing exactly on the windings if possible.

name instance-attribute

name: str

This name should match the name used in the device description ODS

ns cached property

ns: NDArray

[dimensionless] Filament number of turns

resistance instance-attribute

resistance: float

[ohm] total effective resistance of the coil; for superconducting coils, this will be small

rs cached property

rs: NDArray

[m] Filament r-coordinates

self_inductance instance-attribute

self_inductance: float

[H] total scalar self-inductance of this coil

zs cached property

zs: NDArray

[m] Filament z-coordinates

CoilFilament dataclass

CoilFilament(
    r: float, z: float, n: float, self_inductance: float
)

A discretized element of an axisymmetric magnet. Self-inductance is calculated based on conductor geometry.

n instance-attribute

n: float

[dimensionless] number of turns

r instance-attribute

r: float

[m] radial location

self_inductance instance-attribute

self_inductance: float

[H] scalar self-inductance of this filament

z instance-attribute

z: float

[m] z location

device_inductance.circuits

CoilSeriesCircuit dataclass

CoilSeriesCircuit(
    name: str, supply: int, coils: list[tuple[int, float]]
)

A circuit with a single power supply and one or more coils connected in series.

coils instance-attribute

coils: list[tuple[int, float]]

Map of the coils that are in series on this circuit to their series orientation which should be either -1.0 or 1.0

name instance-attribute

name: str

Human-friendly name of the circuit

supply instance-attribute

supply: int

Index of power supply for this circuit

device_inductance.structures

filament

The lowest-level discretization of conducting structure.

PassiveStructureFilament dataclass

PassiveStructureFilament(
    parent_name: str,
    r: float,
    z: float,
    area: float,
    resistance: float,
    self_inductance: float,
    polygon: Polygon,
)

A chunk of a cylindrically-symmetric passive conductor

area instance-attribute
area: float

[m^2] cross-section area

parent_name instance-attribute
parent_name: str

Name of structure that this is associated with

polygon instance-attribute
polygon: Polygon

[m] Outline of mesh element

r instance-attribute
r: float

[m] radial location of filament center

resistance instance-attribute
resistance: float

Ohm] loop resistance

self_inductance instance-attribute
self_inductance: float

[H] self-inductance assuming one full loop

z instance-attribute
z: float

[m] axial location of filament center

heuristics

MAX_EDGE_LENGTH_M module-attribute

MAX_EDGE_LENGTH_M = 0.1

Default maximum length of an edge in the lowest-level discretization

poly_angle

poly_angle(
    p: Polygon,
    centroid: tuple[float, float],
    max_edge_length_m: float = MAX_EDGE_LENGTH_M,
) -> float

Get peak winding number in [rad] of the polygon boundary w.r.t. the centroid. This is not necessarily the same as the maximum angle subtended by the polygon, but is a good heuristic for it, and runs in O(N) instead of O(N^2).

Source code in src/device_inductance/structures/heuristics.py
def poly_angle(
    p: Polygon,
    centroid: tuple[float, float],
    max_edge_length_m: float = MAX_EDGE_LENGTH_M,
) -> float:
    """
    Get peak winding number in [rad] of the polygon boundary w.r.t. the centroid.
    This is not necessarily the same as the maximum angle subtended by the polygon,
    but is a good heuristic for it, and runs in O(N) instead of O(N^2).
    """
    # Expand to 3D and segmentize to avoid missing big sections
    pr, pz = p.boundary.segmentize(max_edge_length_m).xy
    xyz = np.array([pr, np.zeros_like(pr), pz])
    centroid_arr = np.array([centroid[0], 0.0, centroid[1]])

    # Get aggregate winding number w.r.t. centroid along polygon boundary
    _wx, wy, _wz = winding_number(xyz, centroid_arr)

    return np.max(np.abs(wy))  # [rad]

unroundness

unroundness(p: Polygon) -> float

[dimensionless] ratio of perimeter to perimeter of a circle with the same area. Gives a heuristic estimate for how un-like a circle the shape is.

Source code in src/device_inductance/structures/heuristics.py
def unroundness(p: Polygon) -> float:
    """
    [dimensionless] ratio of perimeter to perimeter of a circle with the same area.
    Gives a heuristic estimate for how un-like a circle the shape is.
    """
    perimeter = p.boundary.length  # [m]
    circle_perimeter = 2.0 * np.pi * np.sqrt(p.area / np.pi)  # [m]
    return perimeter / circle_perimeter  # [dimensionless]

winding_number

winding_number(
    path: NDArray, centroid: NDArray
) -> tuple[NDArray, NDArray, NDArray]

Calculate winding number about x, y, and z axes, w.r.t. the filament centroid. This estimates the angle subtended by each segment of the coil path, which is 2*pi radians per full revolution, and may be either positive or negative depending on the orientation of the coil winding path.

The choice of sign convention is arbitrary, as these results exist to compare to each other, and for no other purpose.

Source code in src/device_inductance/structures/heuristics.py
def winding_number(path: NDArray, centroid: NDArray) -> tuple[NDArray, NDArray, NDArray]:
    """
    Calculate winding number about x, y, and z axes, w.r.t. the filament centroid.
    This estimates the angle subtended by each segment of the coil path, which is
    2*pi radians per full revolution, and may be either positive or negative depending
    on the orientation of the coil winding path.

    The choice of sign convention is arbitrary, as these results exist to compare to each other,
    and for no other purpose.
    """

    # path = xyz  # [m]
    n = path.shape[1]

    # centroid = np.mean(path, axis=1)  # [m]
    x_winding = np.zeros(n)  # [rad]
    y_winding = np.zeros(n)
    z_winding = np.zeros(n)

    for i in range(n - 1):
        # Translate the target points to center on the centroid
        first_point = path[:, i] - centroid  # [m]
        second_point = path[:, i + 1] - centroid

        # Extract angle subtended by the two points about each axis
        cross_1_2 = np.cross(first_point, second_point)
        dot_1_2 = np.dot(first_point, second_point)
        x_winding[i + 1] = x_winding[i] + np.atan2(cross_1_2[0], dot_1_2)
        y_winding[i + 1] = y_winding[i] + np.atan2(cross_1_2[1], dot_1_2)
        z_winding[i + 1] = z_winding[i] + np.atan2(cross_1_2[2], dot_1_2)

    return x_winding, y_winding, z_winding  # [rad]

input

CONTACT_DETECTION_DISTANCE module-attribute

CONTACT_DETECTION_DISTANCE = 0.0001

[m] distance under which two passive structure objects are considered to be in electrical contact

PassiveStructureInput dataclass

PassiveStructureInput(
    parent_name: str, polygon: Polygon, resistivity: float
)

First level of structure discretization within device_inductance. Some discretization may have already been performed upstream.

parent_name instance-attribute
parent_name: str

The name of the source of the info from the ODS input

polygon instance-attribute
polygon: Polygon

[m] R-Z cross-section

resistivity instance-attribute
resistivity: float

[ohm-m] (effective) electrical resistivity of material.

loop

Top-level discretization of conducting structure.

PassiveStructureLoop dataclass

PassiveStructureLoop(
    parent_name: str,
    polygon: Polygon,
    filaments: list[PassiveStructureFilament],
)

A logical chunk of structure material that may be the result of slicing a larger object into smaller pieces. Composed of one or more filaments.

ns cached property
ns: NDArray

[dimensionless] (Fractional) number of turns of each filament, weighted according to their cross-sectional area as a fraction of this loop's total.

parent_name instance-attribute
parent_name: str

Name of the input structure this was chunked from

polygon instance-attribute
polygon: Polygon

Shape of the enclosing polygon before sub-discretization

resistance cached property
resistance: float

[ohm] Total loop resistance; effective parallel resistance over all filaments

rs cached property
rs: NDArray

[m] Filament radial coordinates

self_inductance cached property
self_inductance: float

[H] Self-inductance of this loop.

zs cached property
zs: NDArray

[m] Filament axial coordinates

from_input classmethod
from_input(
    inp: PassiveStructureInput,
    slicer: RadialSlicer,
    angle_thresh_deg: float = 20.0,
    unroundness_thresh: float = 2.0,
) -> list[PassiveStructureLoop]

Make one or more PassiveStructureLoop from an input element. The input element may be subdivided into multiple loops if it subtends a large angle relative to the centroid and has a large perimeter-to-area ratio in the section plane.

Source code in src/device_inductance/structures/loop.py
@classmethod
def from_input(
    cls,
    inp: PassiveStructureInput,
    slicer: RadialSlicer,
    angle_thresh_deg: float = 20.0,
    unroundness_thresh: float = 2.0,
) -> list[PassiveStructureLoop]:
    """
    Make one or more PassiveStructureLoop from an input element.
    The input element may be subdivided into multiple loops
    if it subtends a large angle relative to the centroid
    and has a large perimeter-to-area ratio in the section plane.
    """
    # If something spans a large region around the centroid
    # AND it's not a conceptually solid block of material,
    # subdivide it so that we get adequate detail about current in different regions.
    angle_thresh_met = poly_angle(inp.polygon, slicer.centroid) > np.deg2rad(angle_thresh_deg)
    unroundness_thresh_met = unroundness(inp.polygon) > unroundness_thresh

    if angle_thresh_met and unroundness_thresh_met:
        # Slice into angular chunks
        chunks: list[Polygon] = slicer.slice(inp.polygon)
    else:
        # If the original takes up a small angular region or it's a solid block, use it as-is
        chunks: list[Polygon] = [inp.polygon]

    # Each chunk represents a fraction of one contiguous loop;
    # if we were to treat each chunk as a whole loop, the inductance of the system
    # would diverge with increasing discretization.
    # Each sub-loop's fraction of loop is weighted based on its section area.
    return [cls.from_poly(inp.parent_name, p, inp.resistivity) for p in chunks]
from_poly classmethod
from_poly(
    parent_name: str,
    polygon: Polygon,
    resistivity: float,
    max_edge_length_m: float = MAX_EDGE_LENGTH_M,
) -> PassiveStructureLoop

Discretize a structure that is reasonably well-associated with a single centroid and may be the result of earlier sub-division of a larger structure, but may not be fully discretized yet.

Source code in src/device_inductance/structures/loop.py
@classmethod
def from_poly(
    cls,
    parent_name: str,
    polygon: Polygon,
    resistivity: float,
    max_edge_length_m: float = MAX_EDGE_LENGTH_M,
) -> PassiveStructureLoop:
    """
    Discretize a structure that is reasonably well-associated with a single centroid
    and may be the result of earlier sub-division of a larger structure, but may not be
    fully discretized yet.
    """
    # Subdivide by meshing
    sub_polygons: list[Polygon] = mesh._mesh_region(
        np.array(polygon.boundary.segmentize(max_edge_length_m).xy).T
    )

    # Make a filament from each mesh cell
    filaments = [_mesh_elem_to_fil(p, resistivity, parent_name) for p in sub_polygons]

    # Call the collection of filaments a loop
    loop = PassiveStructureLoop(parent_name, polygon, filaments)

    return loop
mutual_inductance
mutual_inductance(other: PassiveStructureLoop) -> float

[H] Mutual inductance between two loops. If self passed as other, the precalculated self-inductance is returned.

Source code in src/device_inductance/structures/loop.py
def mutual_inductance(self, other: PassiveStructureLoop) -> float:
    """
    [H] Mutual inductance between two loops.
    If `self` passed as `other`, the precalculated self-inductance is returned.
    """

    if id(other) == id(self):
        return self.self_inductance

    # Each loop's `ns` includes accounting of both filament number of turns and overall number of turns
    rzn1 = np.array([self.rs, self.zs, self.ns])  # [m], [m], [dimensionless]
    rzn2 = np.array([other.rs, other.zs, other.ns])

    m = cfsem.mutual_inductance_of_cylindrical_coils(rzn1, rzn2)

    return m  # [H]

slicer

Tool for cutting a polyogon into pie-slice chunks

RadialSlicer

RadialSlicer(
    n_slices: int,
    extent: tuple[float, float, float, float],
    centroid: tuple[float, float],
)

Tool for cutting a polyogon into pie-slice chunks

Source code in src/device_inductance/structures/slicer.py
def __init__(
    self,
    n_slices: int,
    extent: tuple[float, float, float, float],
    centroid: tuple[float, float],
):
    # Unpack
    self.centroid = centroid
    rmid, zmid = centroid
    centroid_arr = np.array(centroid)

    # Choose a pie radius
    #   Actualize all the points at the corners of the extent
    extent_points = [np.array(x) for x in product([extent[0], extent[1]], [extent[2], extent[3]])]
    #   Find the largest radius from the centroid to any of the points
    extent_radii = np.linalg.norm(np.array([x - centroid_arr for x in extent_points]), axis=1)
    pie_radius = 1.05 * np.max(extent_radii)

    pie_angles = np.linspace(0.0, 2.0 * np.pi, n_slices + 1)
    pie_r = rmid + pie_radius * np.cos(pie_angles)
    pie_z = zmid + pie_radius * np.sin(pie_angles)
    pie_rz = [x for x in zip(pie_r, pie_z, strict=True)]
    mids = [(rmid, zmid)] * n_slices
    #  Sections start at the midpoint, go to two points on the circle, then end back at the midpoint
    pie_sections = zip(mids, pie_rz[:-1], pie_rz[1:], mids, strict=True)
    pie_slices = [Polygon(s) for s in pie_sections]

    self.polygons = pie_slices
centroid instance-attribute
centroid: tuple[float, float] = centroid

Centroid used to generate the slices

polygons instance-attribute
polygons: list[Polygon] = pie_slices

Contiguous but non-overlapping slices spanning the provided extent

slice
slice(p: Polygon) -> list[Polygon]

Cut polygon p into as many as n_slices new polygons or as few as 1 polygon, if no cuts are needed.

Source code in src/device_inductance/structures/slicer.py
def slice(self, p: Polygon) -> list[Polygon]:
    """Cut polygon `p` into as many as `n_slices` new
    polygons or as few as 1 polygon, if no cuts are needed."""

    # Do the intersections - this can balloon into a jumble of outputs
    intersections = [p.intersection(s) for s in self.polygons]

    # Remove edge and vertex intersections
    new_polygons = [g for g in intersections if isinstance(g, Polygon)]

    for g in intersections:
        # It's possible to have one polygon intersect more than once,
        # at an edge, at a point, etc. and we have to handle all cases.
        if isinstance(g, GeometryCollection | MultiPolygon):
            new_polygons.extend([x for x in g.geoms if isinstance(x, Polygon)])

    # Remove any empty polygons; these are common
    new_polygons = [x for x in new_polygons if not x.is_empty]

    # Make sure the area adds up
    area_orig = p.area
    area_new = sum([x.area for x in new_polygons])
    assert abs(1.0 - abs(area_new / area_orig)) < 1e-3

    return new_polygons

device_inductance.sensors

FullFluxLoop dataclass

FullFluxLoop(name: str, r: float, z: float)

Flux sensor wrapped around the entire toroid, encompassing all the poloidal projected area from this point inward. Idealized here as attached to an ideal integrator s.t. it senses flux integrated in time instead of rate of change. ::

        , - ~ ~ ~ - ,
    , '               ' ,
  ,                       ,         phi
 ,         A_pol           ,         ^
,                           ,        |
,              -----------> * p1     x
,                    R      ,        Z
 ,                         ,
  ,                       ,
    ,                  , '
      ' - , _ _ _ ,  '
             ||
            _||_
            -  +

name instance-attribute

name: str

Sensor name

r instance-attribute

r: float

[m] r-location

z instance-attribute

z: float

[m] z-location

response

response(
    grids: tuple[NDArray, NDArray], psi: NDArray
) -> float

Ideal integrated response to a given flux field

Parameters:

Name Type Description Default
grids tuple[NDArray, NDArray]

[m] 1D arrays describing rectilinear mesh

required
psi NDArray

[Wb], 2D array of rate of change of poloidal flux

required

Raises:

Type Description
ValueError

If the sensor is not inside the mesh

Returns:

Type Description
float

[Wb] ideal integrated response

Source code in src/device_inductance/sensors.py
def response(self, grids: tuple[NDArray, NDArray], psi: NDArray) -> float:
    """Ideal integrated response to a given flux field

    Args:
        grids: [m] 1D arrays describing rectilinear mesh
        psi: [Wb], 2D array of rate of change of poloidal flux

    Raises:
        ValueError: If the sensor is not inside the mesh

    Returns:
        [Wb] ideal integrated response
    """
    # Build interpolators - this is fast but not totally free
    psi_interpolator = MulticubicRectilinear.new(
        [x for x in grids], psi.flatten(), linearize_extrapolation=True
    )  # [V], from [Wb/s] or [V-s/s]

    # Unpack/repack geometry
    point_to_interp = [np.atleast_1d(self.r), np.atleast_1d(self.z)]

    # Check if we are extrapolating the fields
    out_of_bounds_flags = psi_interpolator.check_bounds(point_to_interp, atol=1e-3)
    if np.any(out_of_bounds_flags):
        raise ValueError(
            f"Full flux loop {self.name} outside of mesh bounds."
            + f" Mesh bounds flags: {out_of_bounds_flags}"
        )

    # Interpolate field
    psi_interped = psi_interpolator.eval(point_to_interp)[0]  # [Wb]

    # The response assumes an ideal integrator, so we get a flux value instead of a voltage
    response = psi_interped  # [Wb]

    return response  # [Wb]

PartialFluxLoop dataclass

PartialFluxLoop(
    name: str,
    r0: float,
    phi0: float,
    z0: float,
    r1: float,
    phi1: float,
    z1: float,
    n_discretization: int = 1000,
)

Flux sensor defined in terms of the opposite corners of a rectangle wrapped part-way around the toroid, enclosing some poloidal projected area but no toroidal projected area. Idealized here as attached to an ideal integrator s.t. it senses flux integrated in time instead of rate of change.

::

                           Z
   * --------- * p2        ^
   |           |           |
   |           |           x ---> phi
p1 * --------- *           R

The normal vector is taken as clockwise from the vector formed by the first and second points.

o (r1, z1)
|
|------> normal
|
|
o (r0, z0)

discretized_rz_section cached property

discretized_rz_section: tuple[NDArray, NDArray, NDArray]

Generate evenly-sampled points along the r-z section of the loop and the projected area associated with each point.

Because the projected area at each point is different, 1D trapezoid integration is inconvenient; instead, points are sampled entirely inside the span of the section, not including the endpoints, in order to support a Riemann sum over dot(B, dA) using a finer resolution.

loop_frac property

loop_frac: float

Portion of the full toroid that is swept by this loop

n_discretization class-attribute instance-attribute

n_discretization: int = 1000

Number of points to use for integrating the flux field. Meant to be a constant that is tuned manually to target an acceptable error level.

name instance-attribute

name: str

Sensor name

normal_vector property

normal_vector: NDArray

[dimensionless] (r, z) components of the normal vector of the surface enclosed by the loop

phi0 instance-attribute

phi0: float

[rad] angular location (about z) of first corner

phi1 instance-attribute

phi1: float

[rad] angular location (about z) of first corner

r0 instance-attribute

r0: float

[m] r-location of first corner

r1 instance-attribute

r1: float

[m] r-location of second corner

z0 instance-attribute

z0: float

[m] z-location of first corner

z1 instance-attribute

z1: float

[m] z-location of second corner

response

response(
    grids: tuple[NDArray, NDArray], br: NDArray, bz: NDArray
) -> float

Ideal integrated response to a given local B-field

Parameters:

Name Type Description Default
grids tuple[NDArray, NDArray]

[m] 1D arrays describing rectilinear mesh

required
br NDArray

[T]

required
bz NDArray

[T]

required

Raises:

Type Description
ValueError

If the sensor is not inside the mesh

Returns:

Type Description
float

[Wb] ideal integrated flux through the loop

Source code in src/device_inductance/sensors.py
def response(self, grids: tuple[NDArray, NDArray], br: NDArray, bz: NDArray) -> float:
    """Ideal integrated response to a given local B-field

    Args:
        grids: [m] 1D arrays describing rectilinear mesh
        br: [T]
        bz: [T]

    Raises:
        ValueError: If the sensor is not inside the mesh

    Returns:
        [Wb] ideal integrated flux through the loop
    """

    # Build interpolators - this is fast but not totally free
    br_interpolator = MulticubicRectilinear.new(
        [x for x in grids], br.flatten(), linearize_extrapolation=True
    )
    bz_interpolator = MulticubicRectilinear.new(
        [x for x in grids], bz.flatten(), linearize_extrapolation=True
    )

    # Unpack/repack geometry
    rs, zs, das = self.discretized_rz_section  # [m], [m], [m^2]
    points_to_interp = [rs, zs]  # [m]
    normal = self.normal_vector  # [dimensionless]

    # Check if we are extrapolating the fields
    out_of_bounds_flags = br_interpolator.check_bounds(points_to_interp, atol=1e-3)
    if np.any(out_of_bounds_flags):
        raise ValueError(
            f"Partial flux loop {self.name} outside of mesh bounds."
            + f" Mesh bounds flags: {out_of_bounds_flags}"
        )

    # Interpolate fields
    br_interped = br_interpolator.eval(points_to_interp)  # [T]
    bz_interped = bz_interpolator.eval(points_to_interp)

    # Aggregate response
    # This is the sum of B components aligned with direction of probe;
    # this is _not_ the normed magnitude, it's a sum, because the real probe integrates
    # flux components not normed field
    br_projected = br_interped * normal[0]  # [T]
    bz_projected = bz_interped * normal[1]
    b_projected = br_projected + bz_projected  # [T]

    # The response assumes an ideal integrator, so we get the integrated flux instead of a voltage
    response = float(np.sum(b_projected * das))  # [Wb]

    return response  # [Wb]

PoloidalFieldProbe dataclass

PoloidalFieldProbe(
    name: str,
    r: float,
    phi: float,
    z: float,
    poloidal_angle: float,
)

A poloidal B-field sensor at a point in space; aka Mirnov probe or pickup coil. Idealized here as attached to an ideal integrator s.t. it senses B integrated in time instead of rate of change.

In reality, the hardware senses dB/dt on some vector in the poloidal plane, but this is not so useful in a digital context without high-bandwidth integration.

name instance-attribute

name: str

Sensor name

phi instance-attribute

phi: float

[rad] angular location (about z)

poloidal_angle instance-attribute

poloidal_angle: float

[rad] poloidal orientation (about phi), clockwise from +R-axis

r instance-attribute

r: float

[m] r-location

unit_vector property

unit_vector: NDArray

Get the unit vector in the rz plane that the measurement is projected on. Uses cocos-11 coordinate system, with poloidal angle clockwise about into-the-board toroidal direction, starting from the outboard (r=1, z=0) direction.

z instance-attribute

z: float

[m] z-location

response

response(
    grids: tuple[NDArray, NDArray], br: NDArray, bz: NDArray
) -> float

Ideal integrated response to a given local B-field. This is the sum of local B-field components, not the normed field!

Parameters:

Name Type Description Default
grids tuple[NDArray, NDArray]

[m] 1D arrays describing rectilinear mesh

required
br NDArray

[T]

required
bz NDArray

[T]

required

Raises:

Type Description
ValueError

If the sensor is not inside the mesh

Returns:

Type Description
float

[T] ideal integrated B-field, sum of components projected on the axis of the probe

Source code in src/device_inductance/sensors.py
def response(self, grids: tuple[NDArray, NDArray], br: NDArray, bz: NDArray) -> float:
    """Ideal integrated response to a given local B-field.  This is the sum of local B-field
    components, not the normed field!

    Args:
        grids: [m] 1D arrays describing rectilinear mesh
        br: [T]
        bz: [T]

    Raises:
        ValueError: If the sensor is not inside the mesh

    Returns:
        [T] ideal integrated B-field, sum of components projected on the axis of the probe
    """

    # Build interpolators - this is fast but not totally free
    br_interpolator = MulticubicRectilinear.new(
        [x for x in grids], br.flatten(), linearize_extrapolation=True
    )
    bz_interpolator = MulticubicRectilinear.new(
        [x for x in grids], bz.flatten(), linearize_extrapolation=True
    )

    # Unpack/repack geometry
    unit_vector = self.unit_vector
    point_to_interp = [np.atleast_1d(self.r), np.atleast_1d(self.z)]

    # Check if we are extrapolating the fields
    out_of_bounds_flags = br_interpolator.check_bounds(point_to_interp, atol=1e-3)
    if np.any(out_of_bounds_flags):
        raise ValueError(
            f"Poloidal field probe {self.name} outside of mesh bounds."
            + f" Mesh bounds flags: {out_of_bounds_flags}"
        )

    # Interpolate fields
    br_interped = br_interpolator.eval(point_to_interp)[0]  # [T]
    bz_interped = bz_interpolator.eval(point_to_interp)[0]

    # Aggregate response
    # This is the sum of B components aligned with direction of probe;
    # this is _not_ the normed magnitude, it's a sum, because the real probe integrates
    # flux components not normed field
    br_projected = br_interped * unit_vector[0]  # [T]
    bz_projected = bz_interped * unit_vector[1]
    b_projected = br_projected + bz_projected  # [T]

    # The response is based on an ideal integrator, so we get the
    # exact B-field instead of a voltage
    response = b_projected  # [T]

    return response  # [T]