Field Line Tracing

Curve containers and analysis for field lines and particle trajectories. trace_field_line integrates with fixed-step classical RK4; trace_field_line_adaptive uses Dormand-Prince 5(4) with error-norm step control. Both feed the same geometric analysis (curvature, arc length, mirror points).

traces

Curve containers and analysis for field lines and particle trajectories.

FieldLine dataclass

Ordered sequence of 3D points tracing a vector field at a fixed time.

A field line is parameterized by arc length, not time. It represents the spatial structure of a vector field (typically \(\mathbf{B}\)) at a single instant.

Parameters:

Name	Type	Description	Default
points	FloatArray	Ordered positions along the line, shape (N, 3) with N >= 2.	required
field_name	str	Name of the traced vector field (e.g. "B").	required
seed_point	Vector3	Starting point for the trace.	required
normalization	Normalization	Unit conversion for this data.	required
time	float | None	Simulation time at which the field was sampled.	None
step	int | None	Timestep index.	None
direction	str	Trace direction: "both", "forward", or "backward".	'both'
scalars	dict[str, FloatArray]	Named scalar quantities sampled along the line, each shape (N,).	dict()
metadata	dict[str, Any]	Arbitrary metadata (reader info, integration parameters, etc.).	dict()

Examples:

>>> import numpy as np
>>> from pypic.units import Normalization
>>> pts = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]])
>>> fl = FieldLine(
...     points=pts, field_name="B", seed_point=(0.0, 0.0, 0.0),
...     normalization=Normalization.identity(),
... )
>>> fl.n_points
3
>>> fl.start_point
(0.0, 0.0, 0.0)
>>> fl.end_point
(2.0, 0.0, 0.0)
>>> len(fl)
3

Source code in src/pypic/traces/_fieldline.py

@dataclass(frozen=True, slots=True)
class FieldLine:
    r"""Ordered sequence of 3D points tracing a vector field at a fixed time.

    A field line is parameterized by arc length, not time. It represents
    the spatial structure of a vector field (typically $\mathbf{B}$) at a
    single instant.

    Parameters
    ----------
    points : FloatArray
        Ordered positions along the line, shape ``(N, 3)`` with ``N >= 2``.
    field_name : str
        Name of the traced vector field (e.g. ``"B"``).
    seed_point : Vector3
        Starting point for the trace.
    normalization : Normalization
        Unit conversion for this data.
    time : float | None
        Simulation time at which the field was sampled.
    step : int | None
        Timestep index.
    direction : str
        Trace direction: ``"both"``, ``"forward"``, or ``"backward"``.
    scalars : dict[str, FloatArray]
        Named scalar quantities sampled along the line, each shape ``(N,)``.
    metadata : dict[str, Any]
        Arbitrary metadata (reader info, integration parameters, etc.).

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.units import Normalization
    >>> pts = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]])
    >>> fl = FieldLine(
    ...     points=pts, field_name="B", seed_point=(0.0, 0.0, 0.0),
    ...     normalization=Normalization.identity(),
    ... )
    >>> fl.n_points
    3
    >>> fl.start_point
    (0.0, 0.0, 0.0)
    >>> fl.end_point
    (2.0, 0.0, 0.0)
    >>> len(fl)
    3
    """

    points: FloatArray
    field_name: str
    seed_point: Vector3
    normalization: Normalization
    time: float | None = None
    step: int | None = None
    direction: TraceDirection = "both"
    scalars: dict[str, FloatArray] = field(default_factory=dict)
    metadata: dict[str, Any] = field(default_factory=dict)

    def __post_init__(self) -> None:
        if self.points.ndim != 2 or self.points.shape[1] != 3:
            msg = f"points must have shape (N, 3), got {self.points.shape}"
            raise ValueError(msg)
        if self.points.shape[0] < 2:
            msg = f"points must have at least 2 rows, got {self.points.shape[0]}"
            raise ValueError(msg)
        if self.direction not in _VALID_DIRECTIONS:
            msg = (
                f"direction must be one of {sorted(_VALID_DIRECTIONS)}, "
                f"got {self.direction!r}"
            )
            raise ValueError(msg)
        n = self.points.shape[0]
        for name, arr in self.scalars.items():
            if arr.shape != (n,):
                msg = f"scalar {name!r} has shape {arr.shape}, expected ({n},)"
                raise ValueError(msg)
        object.__setattr__(self, "scalars", MappingProxyType(dict(self.scalars)))
        object.__setattr__(self, "metadata", MappingProxyType(dict(self.metadata)))

    @property
    def n_points(self) -> int:
        """Number of points along the field line."""
        return int(self.points.shape[0])

    @property
    def start_point(self) -> Vector3:
        """First point on the line as a 3-tuple."""
        p = self.points[0]
        return (float(p[0]), float(p[1]), float(p[2]))

    @property
    def end_point(self) -> Vector3:
        """Last point on the line as a 3-tuple."""
        p = self.points[-1]
        return (float(p[0]), float(p[1]), float(p[2]))

    def __len__(self) -> int:
        return int(self.points.shape[0])

    def __repr__(self) -> str:
        scalar_names = sorted(self.scalars) if self.scalars else []
        parts = [
            f"FieldLine({self.field_name!r}, n={self.n_points}",
            f"direction={self.direction!r}",
        ]
        if self.time is not None:
            parts.append(f"t={self.time}")
        if scalar_names:
            parts.append(f"scalars={scalar_names}")
        return ", ".join(parts) + ")"

    def with_scalars(self, **new_scalars: FloatArray) -> FieldLine:
        """Return a new FieldLine with additional or replaced scalars.

        Parameters
        ----------
        **new_scalars : FloatArray
            Scalar arrays to merge, each shape ``(N,)``.

        Returns
        -------
        FieldLine
        """
        n = self.n_points
        for name, arr in new_scalars.items():
            if arr.shape != (n,):
                msg = f"Scalar {name!r} has shape {arr.shape}, expected ({n},)"
                raise ValueError(msg)
        merged = dict(self.scalars)
        merged.update(new_scalars)
        return copy.replace(self, scalars=merged)

n_points property

Number of points along the field line.

start_point property

First point on the line as a 3-tuple.

end_point property

Last point on the line as a 3-tuple.

with_scalars(**new_scalars)

Return a new FieldLine with additional or replaced scalars.

Parameters:

Name	Type	Description	Default
**new_scalars	FloatArray	Scalar arrays to merge, each shape (N,).	{}

Returns:

Type	Description
FieldLine

Source code in src/pypic/traces/_fieldline.py

def with_scalars(self, **new_scalars: FloatArray) -> FieldLine:
    """Return a new FieldLine with additional or replaced scalars.

    Parameters
    ----------
    **new_scalars : FloatArray
        Scalar arrays to merge, each shape ``(N,)``.

    Returns
    -------
    FieldLine
    """
    n = self.n_points
    for name, arr in new_scalars.items():
        if arr.shape != (n,):
            msg = f"Scalar {name!r} has shape {arr.shape}, expected ({n},)"
            raise ValueError(msg)
    merged = dict(self.scalars)
    merged.update(new_scalars)
    return copy.replace(self, scalars=merged)

ParticleTrace dataclass

Ordered sequence of spacetime points tracing a particle's worldline.

Unlike FieldLine, a particle trace is parameterized by time (monotonic) and carries velocity at each point. Represents a physical particle's trajectory through the simulation domain.

Parameters:

Name	Type	Description	Default
points	FloatArray	Ordered positions, shape (N, 3) with N >= 2.	required
time	FloatArray	Time at each point, shape (N,), monotonically increasing.	required
velocity	FloatArray	Velocity at each point, shape (N, 3).	required
species_name	str	Species name (e.g. "electrons").	required
normalization	Normalization	Unit conversion for this data.	required
species	SpeciesInfo | None	Species charge/mass info for derived calculations.	None
particle_id	int | None	Tracking ID, if available.	None
scalars	dict[str, FloatArray]	Named scalar quantities along the trajectory, each shape (N,).	dict()
metadata	dict[str, Any]	Arbitrary metadata.	dict()

Examples:

>>> import numpy as np
>>> from pypic.units import Normalization
>>> pts = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]])
>>> t = np.array([0.0, 0.5, 1.0])
>>> vel = np.array([[2.0, 0.0, 0.0], [2.0, 0.0, 0.0], [2.0, 0.0, 0.0]])
>>> tr = ParticleTrace(
...     points=pts, time=t, velocity=vel, species_name="electrons",
...     normalization=Normalization.identity(),
... )
>>> tr.n_points
3
>>> tr.duration
1.0
>>> len(tr)
3

Source code in src/pypic/traces/_particletrace.py

@dataclass(frozen=True, slots=True)
class ParticleTrace:
    r"""Ordered sequence of spacetime points tracing a particle's worldline.

    Unlike `FieldLine`, a particle trace is parameterized by time (monotonic)
    and carries velocity at each point. Represents a physical particle's
    trajectory through the simulation domain.

    Parameters
    ----------
    points : FloatArray
        Ordered positions, shape ``(N, 3)`` with ``N >= 2``.
    time : FloatArray
        Time at each point, shape ``(N,)``, monotonically increasing.
    velocity : FloatArray
        Velocity at each point, shape ``(N, 3)``.
    species_name : str
        Species name (e.g. ``"electrons"``).
    normalization : Normalization
        Unit conversion for this data.
    species : SpeciesInfo | None
        Species charge/mass info for derived calculations.
    particle_id : int | None
        Tracking ID, if available.
    scalars : dict[str, FloatArray]
        Named scalar quantities along the trajectory, each shape ``(N,)``.
    metadata : dict[str, Any]
        Arbitrary metadata.

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.units import Normalization
    >>> pts = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]])
    >>> t = np.array([0.0, 0.5, 1.0])
    >>> vel = np.array([[2.0, 0.0, 0.0], [2.0, 0.0, 0.0], [2.0, 0.0, 0.0]])
    >>> tr = ParticleTrace(
    ...     points=pts, time=t, velocity=vel, species_name="electrons",
    ...     normalization=Normalization.identity(),
    ... )
    >>> tr.n_points
    3
    >>> tr.duration
    1.0
    >>> len(tr)
    3
    """

    points: FloatArray
    time: FloatArray
    velocity: FloatArray
    species_name: str
    normalization: Normalization
    species: SpeciesInfo | None = None
    particle_id: int | None = None
    scalars: dict[str, FloatArray] = field(default_factory=dict)
    metadata: dict[str, Any] = field(default_factory=dict)

    def __post_init__(self) -> None:
        if self.points.ndim != 2 or self.points.shape[1] != 3:
            msg = f"points must have shape (N, 3), got {self.points.shape}"
            raise ValueError(msg)
        n = self.points.shape[0]
        if n < 2:
            msg = f"points must have at least 2 rows, got {n}"
            raise ValueError(msg)
        if self.time.shape != (n,):
            msg = f"time must have shape ({n},), got {self.time.shape}"
            raise ValueError(msg)
        if self.velocity.ndim != 2 or self.velocity.shape != (n, 3):
            msg = f"velocity must have shape ({n}, 3), got {self.velocity.shape}"
            raise ValueError(msg)
        if not np.all(np.diff(self.time) > 0):
            msg = "time must be strictly monotonically increasing"
            raise ValueError(msg)
        for name, arr in self.scalars.items():
            if arr.shape != (n,):
                msg = f"scalar {name!r} has shape {arr.shape}, expected ({n},)"
                raise ValueError(msg)
        object.__setattr__(self, "scalars", MappingProxyType(dict(self.scalars)))
        object.__setattr__(self, "metadata", MappingProxyType(dict(self.metadata)))

    @property
    def n_points(self) -> int:
        """Number of points along the trajectory."""
        return int(self.points.shape[0])

    @property
    def start_point(self) -> Vector3:
        """First position as a 3-tuple."""
        p = self.points[0]
        return (float(p[0]), float(p[1]), float(p[2]))

    @property
    def end_point(self) -> Vector3:
        """Last position as a 3-tuple."""
        p = self.points[-1]
        return (float(p[0]), float(p[1]), float(p[2]))

    @property
    def start_time(self) -> float:
        """Time at first point."""
        return float(self.time[0])

    @property
    def end_time(self) -> float:
        """Time at last point."""
        return float(self.time[-1])

    @property
    def duration(self) -> float:
        """Total elapsed time."""
        return float(self.time[-1] - self.time[0])

    def __len__(self) -> int:
        return int(self.points.shape[0])

    def __repr__(self) -> str:
        scalar_names = sorted(self.scalars) if self.scalars else []
        parts = [
            f"ParticleTrace({self.species_name!r}, n={self.n_points}",
            f"t=[{self.start_time:.4g}, {self.end_time:.4g}]",
        ]
        if self.particle_id is not None:
            parts.append(f"id={self.particle_id}")
        if scalar_names:
            parts.append(f"scalars={scalar_names}")
        return ", ".join(parts) + ")"

    def with_scalars(self, **new_scalars: FloatArray) -> ParticleTrace:
        """Return a new ParticleTrace with additional or replaced scalars.

        Parameters
        ----------
        **new_scalars : FloatArray
            Scalar arrays to merge, each shape ``(N,)``.

        Returns
        -------
        ParticleTrace
        """
        n = self.n_points
        for name, arr in new_scalars.items():
            if arr.shape != (n,):
                msg = f"Scalar {name!r} has shape {arr.shape}, expected ({n},)"
                raise ValueError(msg)
        merged = dict(self.scalars)
        merged.update(new_scalars)
        return copy.replace(self, scalars=merged)

n_points property

Number of points along the trajectory.

start_point property

First position as a 3-tuple.

end_point property

Last position as a 3-tuple.

start_time property

Time at first point.

end_time property

Time at last point.

duration property

Total elapsed time.

with_scalars(**new_scalars)

Return a new ParticleTrace with additional or replaced scalars.

Parameters:

Name	Type	Description	Default
**new_scalars	FloatArray	Scalar arrays to merge, each shape (N,).	{}

Returns:

Type	Description
ParticleTrace

Source code in src/pypic/traces/_particletrace.py

def with_scalars(self, **new_scalars: FloatArray) -> ParticleTrace:
    """Return a new ParticleTrace with additional or replaced scalars.

    Parameters
    ----------
    **new_scalars : FloatArray
        Scalar arrays to merge, each shape ``(N,)``.

    Returns
    -------
    ParticleTrace
    """
    n = self.n_points
    for name, arr in new_scalars.items():
        if arr.shape != (n,):
            msg = f"Scalar {name!r} has shape {arr.shape}, expected ({n},)"
            raise ValueError(msg)
    merged = dict(self.scalars)
    merged.update(new_scalars)
    return copy.replace(self, scalars=merged)

PoincareSection dataclass

Result of poincare_section: per-seed punctures and provenance.

Holds both the 3D crossing positions (for re-projection onto a different surface) and the projected 2D coordinates (for plotting). The underlying FieldLine traces are retained so the same trajectories can be re-punctured against a different surface without re-integrating.

Parameters:

Name	Type	Description	Default
surface	PoincareSurface	The plane the punctures were extracted on.	required
seeds	FloatArray	Original seed positions, shape (M, 3).	required
direction	str	"forward", "backward", or "both" — propagated from the trace call.	required
punctures_3d	tuple[FloatArray, ...]	Per-seed crossings in 3D, each shape (n_k, 3).	required
punctures_2d	tuple[FloatArray, ...]	Per-seed crossings in plane coordinates, each (n_k, 2).	required
field_lines	tuple[FieldLine, ...]	Underlying trace per seed.	required
metadata	dict[str, Any]	Provenance: n_crossings_per_seed, n_steps_per_seed, termination_reasons, surface_name.	dict()

Examples:

>>> import numpy as np
>>> from pypic.units import Normalization
>>> from pypic.dataset import FieldDataset
>>> from pypic.grid import GridInfo
>>> grid = GridInfo(dimensions=(16, 16, 16), spacing=(0.2, 0.2, 0.2),
...                 origin=(-1.5, -1.5, -1.5))
>>> x = np.linspace(-1.5, 1.5, 16)
>>> y = np.linspace(-1.5, 1.5, 16)
>>> z = np.linspace(-1.5, 1.5, 16)
>>> X, Y, _ = np.meshgrid(x, y, z, indexing="ij")
>>> data = FieldDataset.from_arrays(
...     {"B_1": -Y, "B_2": X, "B_3": 0.1 * np.ones_like(X)},
...     grid, Normalization.identity(),
... )
>>> surf = PoincareSurface.from_axis("y", 0.0)
>>> section = poincare_section(
...     data, np.array([[0.6, 0.0, 0.0]]), surf, max_steps=2000,
...     direction="forward",
... )
>>> bool(section.metadata["n_crossings_per_seed"][0] > 1)
True

Source code in src/pypic/traces/_poincare.py

@dataclass(frozen=True)
class PoincareSection:
    r"""Result of `poincare_section`: per-seed punctures and provenance.

    Holds both the 3D crossing positions (for re-projection onto a
    different surface) and the projected 2D coordinates (for plotting).
    The underlying `FieldLine` traces are retained so the same
    trajectories can be re-punctured against a different surface without
    re-integrating.

    Parameters
    ----------
    surface : PoincareSurface
        The plane the punctures were extracted on.
    seeds : FloatArray
        Original seed positions, shape ``(M, 3)``.
    direction : str
        ``"forward"``, ``"backward"``, or ``"both"`` — propagated from
        the trace call.
    punctures_3d : tuple[FloatArray, ...]
        Per-seed crossings in 3D, each shape ``(n_k, 3)``.
    punctures_2d : tuple[FloatArray, ...]
        Per-seed crossings in plane coordinates, each ``(n_k, 2)``.
    field_lines : tuple[FieldLine, ...]
        Underlying trace per seed.
    metadata : dict[str, Any]
        Provenance: ``n_crossings_per_seed``, ``n_steps_per_seed``,
        ``termination_reasons``, ``surface_name``.

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.units import Normalization
    >>> from pypic.dataset import FieldDataset
    >>> from pypic.grid import GridInfo
    >>> grid = GridInfo(dimensions=(16, 16, 16), spacing=(0.2, 0.2, 0.2),
    ...                 origin=(-1.5, -1.5, -1.5))
    >>> x = np.linspace(-1.5, 1.5, 16)
    >>> y = np.linspace(-1.5, 1.5, 16)
    >>> z = np.linspace(-1.5, 1.5, 16)
    >>> X, Y, _ = np.meshgrid(x, y, z, indexing="ij")
    >>> data = FieldDataset.from_arrays(
    ...     {"B_1": -Y, "B_2": X, "B_3": 0.1 * np.ones_like(X)},
    ...     grid, Normalization.identity(),
    ... )
    >>> surf = PoincareSurface.from_axis("y", 0.0)
    >>> section = poincare_section(
    ...     data, np.array([[0.6, 0.0, 0.0]]), surf, max_steps=2000,
    ...     direction="forward",
    ... )
    >>> bool(section.metadata["n_crossings_per_seed"][0] > 1)
    True
    """

    surface: PoincareSurface
    seeds: FloatArray
    direction: TraceDirection
    punctures_3d: tuple[FloatArray, ...]
    punctures_2d: tuple[FloatArray, ...]
    field_lines: tuple[FieldLine, ...]
    # Read-only after __post_init__ rewraps to MappingProxyType. Accept
    # a plain dict on construction for callers' convenience.
    metadata: Mapping[str, Any] = field(default_factory=dict)

    def __post_init__(self) -> None:
        object.__setattr__(self, "metadata", MappingProxyType(dict(self.metadata)))

    @cached_property
    def n_seeds(self) -> int:
        """Number of seeds in this section."""
        return int(self.seeds.shape[0])

    @cached_property
    def all_punctures_2d(self) -> FloatArray:
        r"""Concatenated $(\sum_k n_k, 2)$ cloud across all seeds."""
        if not self.punctures_2d:
            return np.empty((0, 2), dtype=np.float64)
        return np.concatenate(self.punctures_2d, axis=0)

    @cached_property
    def all_punctures_3d(self) -> FloatArray:
        r"""Concatenated $(\sum_k n_k, 3)$ cloud across all seeds."""
        if not self.punctures_3d:
            return np.empty((0, 3), dtype=np.float64)
        return np.concatenate(self.punctures_3d, axis=0)

n_seeds cached property

Number of seeds in this section.

all_punctures_2d cached property

Concatenated \((\sum_k n_k, 2)\) cloud across all seeds.

all_punctures_3d cached property

Concatenated \((\sum_k n_k, 3)\) cloud across all seeds.

PoincareSurface dataclass

Transverse plane \(\Sigma\) for a Poincaré section.

\(\Sigma = \{\mathbf{x} : \hat{\mathbf{n}} \cdot (\mathbf{x} - \mathbf{p}) = 0\}\), parameterized by an outward normal \(\mathbf{n}\) and an in-plane reference point \(\mathbf{p}\). The 2D plane coordinates \((u, v)\) produced by project are measured relative to \(\mathbf{p}\) in an orthonormal basis basis_2d spanning \(\Sigma\).

Parameters:

Name	Type	Description	Default
normal	Vector3	Plane normal. Normalized internally; need not be unit length.	required
point	Vector3	A point on the plane (origin of the 2D \((u, v)\) frame).	required
name	str | None	Optional label propagated into output metadata and plots.	None

Examples:

>>> surf = PoincareSurface.from_axis("y", 0.0, name="meridional")
>>> surf.offset
0.0
>>> import numpy as np
>>> u, v = surf.basis_2d
>>> float(np.dot(u, v))
0.0

Source code in src/pypic/traces/_poincare.py

@dataclass(frozen=True)
class PoincareSurface:
    r"""Transverse plane $\Sigma$ for a Poincaré section.

    $\Sigma = \{\mathbf{x} : \hat{\mathbf{n}} \cdot
    (\mathbf{x} - \mathbf{p}) = 0\}$, parameterized by an outward
    normal $\mathbf{n}$ and an in-plane reference point $\mathbf{p}$.
    The 2D plane coordinates $(u, v)$ produced by `project` are
    measured relative to $\mathbf{p}$ in an orthonormal basis
    `basis_2d` spanning $\Sigma$.

    Parameters
    ----------
    normal : Vector3
        Plane normal. Normalized internally; need not be unit length.
    point : Vector3
        A point on the plane (origin of the 2D $(u, v)$ frame).
    name : str | None
        Optional label propagated into output metadata and plots.

    Examples
    --------
    >>> surf = PoincareSurface.from_axis("y", 0.0, name="meridional")
    >>> surf.offset
    0.0
    >>> import numpy as np
    >>> u, v = surf.basis_2d
    >>> float(np.dot(u, v))
    0.0
    """

    normal: Vector3
    point: Vector3
    name: str | None = None

    @classmethod
    def from_axis(
        cls,
        axis: Literal["x", "y", "z"],
        value: float,
        *,
        name: str | None = None,
    ) -> Self:
        r"""Build an axis-aligned surface $x_i = \text{value}$.

        Covers tokamak poloidal sections (``from_axis("y", 0.0)`` for a
        $\phi = 0$ cut, with the toroidal direction along $\hat y$) and
        the standard magnetotail meridional plane
        (``from_axis("y", 0.0)`` in GSM).

        Examples
        --------
        >>> surf = PoincareSurface.from_axis("z", 1.5)
        >>> surf.normal
        (0.0, 0.0, 1.0)
        >>> surf.point
        (0.0, 0.0, 1.5)
        """
        if axis not in _AXIS_INDEX:
            msg = f"axis must be one of {sorted(_AXIS_INDEX)}, got {axis!r}"
            raise ValueError(msg)
        i = _AXIS_INDEX[axis]
        normal = [0.0, 0.0, 0.0]
        normal[i] = 1.0
        point = [0.0, 0.0, 0.0]
        point[i] = float(value)
        return cls(
            normal=(normal[0], normal[1], normal[2]),
            point=(point[0], point[1], point[2]),
            name=name,
        )

    @cached_property
    def _normal_arr(self) -> FloatArray:
        n = np.asarray(self.normal, dtype=np.float64)
        norm = float(np.linalg.norm(n))
        if norm == 0.0:
            msg = "PoincareSurface.normal must be non-zero"
            raise ValueError(msg)
        return n / norm

    @cached_property
    def _point_arr(self) -> FloatArray:
        return np.asarray(self.point, dtype=np.float64)

    @cached_property
    def offset(self) -> float:
        r"""Signed scalar $d = \hat{\mathbf{n}} \cdot \mathbf{p}$.

        The form consumed by
        [`pypic.traces.plane_crossings`][pypic.traces.plane_crossings].
        """
        return float(self._normal_arr @ self._point_arr)

    @cached_property
    def basis_2d(self) -> tuple[FloatArray, FloatArray]:
        r"""Orthonormal $(\hat{\mathbf{u}}, \hat{\mathbf{v}})$ spanning $\Sigma$.

        Constructed via Gram--Schmidt against the world axis least
        aligned with $\hat{\mathbf{n}}$ for numerical stability — the
        standard oblique-section basis shared by FLARE
        [@Frerichs2024] and most field-mapping tools.
        """
        n = self._normal_arr
        # Pick the world axis least parallel to n
        axis = int(np.argmin(np.abs(n)))
        seed = np.zeros(3, dtype=np.float64)
        seed[axis] = 1.0
        u = np.cross(n, seed)
        u = u / np.linalg.norm(u)
        v = np.cross(n, u)
        # cross of two unit-orthogonal vectors is already unit-length;
        # renormalize defensively against accumulated FP error.
        v = v / np.linalg.norm(v)
        return u, v

    def project(self, points_3d: FloatArray) -> FloatArray:
        r"""Project 3D points onto plane coordinates $(u, v)$.

        For each input point $\mathbf{x}$:
        $u = \hat{\mathbf{u}} \cdot (\mathbf{x} - \mathbf{p})$,
        $v = \hat{\mathbf{v}} \cdot (\mathbf{x} - \mathbf{p})$.

        Parameters
        ----------
        points_3d : FloatArray
            Input points, shape ``(M, 3)``. Empty ``(0, 3)`` allowed.

        Returns
        -------
        FloatArray
            Plane coordinates, shape ``(M, 2)``.

        Examples
        --------
        Plane coordinates are oriented by the Gram--Schmidt basis; for the
        z-axis surface the basis is $\hat{\mathbf{u}} = (0, 1, 0)$,
        $\hat{\mathbf{v}} = (-1, 0, 0)$ (right-handed about $\hat{\mathbf{n}}$):

        >>> import numpy as np
        >>> surf = PoincareSurface.from_axis("z", 0.0)
        >>> pts = np.array([[1.0, 2.0, 0.0], [-1.0, 3.0, 0.0]])
        >>> surf.project(pts)
        array([[ 2., -1.],
               [ 3.,  1.]])
        """
        if points_3d.shape[0] == 0:
            return np.empty((0, 2), dtype=np.float64)
        u, v = self.basis_2d
        rel = points_3d - self._point_arr
        return np.column_stack([rel @ u, rel @ v])

offset cached property

Signed scalar \(d = \hat{\mathbf{n}} \cdot \mathbf{p}\).

The form consumed by pypic.traces.plane_crossings.

basis_2d cached property

Orthonormal \((\hat{\mathbf{u}}, \hat{\mathbf{v}})\) spanning \(\Sigma\).

Constructed via Gram--Schmidt against the world axis least aligned with \(\hat{\mathbf{n}}\) for numerical stability — the standard oblique-section basis shared by FLARE [@Frerichs2024] and most field-mapping tools.

from_axis(axis, value, *, name=None) classmethod

Build an axis-aligned surface \(x_i = \text{value}\).

Covers tokamak poloidal sections (from_axis("y", 0.0) for a \(\phi = 0\) cut, with the toroidal direction along \(\hat y\)) and the standard magnetotail meridional plane (from_axis("y", 0.0) in GSM).

Examples:

>>> surf = PoincareSurface.from_axis("z", 1.5)
>>> surf.normal
(0.0, 0.0, 1.0)
>>> surf.point
(0.0, 0.0, 1.5)

Source code in src/pypic/traces/_poincare.py

@classmethod
def from_axis(
    cls,
    axis: Literal["x", "y", "z"],
    value: float,
    *,
    name: str | None = None,
) -> Self:
    r"""Build an axis-aligned surface $x_i = \text{value}$.

    Covers tokamak poloidal sections (``from_axis("y", 0.0)`` for a
    $\phi = 0$ cut, with the toroidal direction along $\hat y$) and
    the standard magnetotail meridional plane
    (``from_axis("y", 0.0)`` in GSM).

    Examples
    --------
    >>> surf = PoincareSurface.from_axis("z", 1.5)
    >>> surf.normal
    (0.0, 0.0, 1.0)
    >>> surf.point
    (0.0, 0.0, 1.5)
    """
    if axis not in _AXIS_INDEX:
        msg = f"axis must be one of {sorted(_AXIS_INDEX)}, got {axis!r}"
        raise ValueError(msg)
    i = _AXIS_INDEX[axis]
    normal = [0.0, 0.0, 0.0]
    normal[i] = 1.0
    point = [0.0, 0.0, 0.0]
    point[i] = float(value)
    return cls(
        normal=(normal[0], normal[1], normal[2]),
        point=(point[0], point[1], point[2]),
        name=name,
    )

project(points_3d)

Project 3D points onto plane coordinates \((u, v)\).

For each input point \(\mathbf{x}\): \(u = \hat{\mathbf{u}} \cdot (\mathbf{x} - \mathbf{p})\), \(v = \hat{\mathbf{v}} \cdot (\mathbf{x} - \mathbf{p})\).

Parameters:

Name	Type	Description	Default
points_3d	FloatArray	Input points, shape (M, 3). Empty (0, 3) allowed.	required

Returns:

Type	Description
FloatArray	Plane coordinates, shape (M, 2).

Examples:

Plane coordinates are oriented by the Gram--Schmidt basis; for the z-axis surface the basis is \(\hat{\mathbf{u}} = (0, 1, 0)\), \(\hat{\mathbf{v}} = (-1, 0, 0)\) (right-handed about \(\hat{\mathbf{n}}\)):

>>> import numpy as np
>>> surf = PoincareSurface.from_axis("z", 0.0)
>>> pts = np.array([[1.0, 2.0, 0.0], [-1.0, 3.0, 0.0]])
>>> surf.project(pts)
array([[ 2., -1.],
       [ 3.,  1.]])

Source code in src/pypic/traces/_poincare.py

def project(self, points_3d: FloatArray) -> FloatArray:
    r"""Project 3D points onto plane coordinates $(u, v)$.

    For each input point $\mathbf{x}$:
    $u = \hat{\mathbf{u}} \cdot (\mathbf{x} - \mathbf{p})$,
    $v = \hat{\mathbf{v}} \cdot (\mathbf{x} - \mathbf{p})$.

    Parameters
    ----------
    points_3d : FloatArray
        Input points, shape ``(M, 3)``. Empty ``(0, 3)`` allowed.

    Returns
    -------
    FloatArray
        Plane coordinates, shape ``(M, 2)``.

    Examples
    --------
    Plane coordinates are oriented by the Gram--Schmidt basis; for the
    z-axis surface the basis is $\hat{\mathbf{u}} = (0, 1, 0)$,
    $\hat{\mathbf{v}} = (-1, 0, 0)$ (right-handed about $\hat{\mathbf{n}}$):

    >>> import numpy as np
    >>> surf = PoincareSurface.from_axis("z", 0.0)
    >>> pts = np.array([[1.0, 2.0, 0.0], [-1.0, 3.0, 0.0]])
    >>> surf.project(pts)
    array([[ 2., -1.],
           [ 3.,  1.]])
    """
    if points_3d.shape[0] == 0:
        return np.empty((0, 2), dtype=np.float64)
    u, v = self.basis_2d
    rel = points_3d - self._point_arr
    return np.column_stack([rel @ u, rel @ v])

TerminationReason

Bases: StrEnum

Why a field line trace stopped.

Source code in src/pypic/traces/_tracing.py

class TerminationReason(StrEnum):
    """Why a field line trace stopped."""

    MAX_STEPS = "max_steps"  # reached step limit
    DOMAIN_EXIT = "domain_exit"  # left interpolation domain (NaN)
    NULL_POINT = "null_point"  # |B| below null_threshold
    CALLBACK = "callback"  # user terminate() returned True
    CLOSED_LOOP = "closed_loop"  # trace re-entered a loop_tol ball of a past point

VectorFieldInterpolator dataclass

Pre-built trilinear interpolator for a 3-component vector field.

Wraps a single RegularGridInterpolator over a stacked (..., 3) value array so each __call__ dispatches once instead of three times. Built once and reused across all RK4 stages and seed points.

Source code in src/pypic/traces/_tracing.py

@dataclass(frozen=True, slots=True)
class VectorFieldInterpolator:
    r"""Pre-built trilinear interpolator for a 3-component vector field.

    Wraps a single ``RegularGridInterpolator`` over a stacked
    ``(..., 3)`` value array so each ``__call__`` dispatches once
    instead of three times. Built once and reused across all RK4
    stages and seed points.
    """

    _interp: RegularGridInterpolator

    @classmethod
    def from_dataset(
        cls,
        data: FieldDataset,
        components: tuple[str, str, str] = ("B_1", "B_2", "B_3"),
    ) -> Self:
        """Build from a FieldDataset.

        Parameters
        ----------
        data : FieldDataset
            Gridded field data.
        components : tuple[str, str, str]
            Names of the three vector components.

        Returns
        -------
        VectorFieldInterpolator
        """
        coord_arrays = data.grid.coordinate_arrays()
        stacked = np.stack(
            [np.asarray(data[c], dtype=np.float64) for c in components],
            axis=-1,
        )
        interp = RegularGridInterpolator(
            coord_arrays,
            stacked,
            method="linear",
            bounds_error=False,
            fill_value=np.nan,
        )
        return cls(_interp=interp)

    def __call__(self, point: FloatArray) -> FloatArray:
        """Evaluate the vector field at a single point.

        Parameters
        ----------
        point : FloatArray
            Position, shape ``(3,)``.

        Returns
        -------
        FloatArray
            Field vector, shape ``(3,)``. NaN if outside domain.
        """
        return self._interp(point.reshape(1, 3))[0]  # type: ignore[no-any-return]

    def batch(self, points: FloatArray) -> FloatArray:
        """Evaluate the vector field at ``N`` points in one call.

        Used by the batched adaptive tracer so all seeds in a step share
        a single ``RegularGridInterpolator`` dispatch. Roughly N× faster
        than calling `__call__` N times because per-call Python /
        argument-marshaling overhead amortizes over the batch.

        Parameters
        ----------
        points : FloatArray
            Positions, shape ``(N, 3)``.

        Returns
        -------
        FloatArray
            Field vectors, shape ``(N, 3)``. Rows are NaN where the
            corresponding point is outside the interpolation domain.
        """
        return self._interp(points)

from_dataset(data, components=('B_1', 'B_2', 'B_3')) classmethod

Build from a FieldDataset.

Parameters:

Name	Type	Description	Default
data	FieldDataset	Gridded field data.	required
components	tuple[str, str, str]	Names of the three vector components.	('B_1', 'B_2', 'B_3')

Returns:

Type	Description
VectorFieldInterpolator

Source code in src/pypic/traces/_tracing.py

@classmethod
def from_dataset(
    cls,
    data: FieldDataset,
    components: tuple[str, str, str] = ("B_1", "B_2", "B_3"),
) -> Self:
    """Build from a FieldDataset.

    Parameters
    ----------
    data : FieldDataset
        Gridded field data.
    components : tuple[str, str, str]
        Names of the three vector components.

    Returns
    -------
    VectorFieldInterpolator
    """
    coord_arrays = data.grid.coordinate_arrays()
    stacked = np.stack(
        [np.asarray(data[c], dtype=np.float64) for c in components],
        axis=-1,
    )
    interp = RegularGridInterpolator(
        coord_arrays,
        stacked,
        method="linear",
        bounds_error=False,
        fill_value=np.nan,
    )
    return cls(_interp=interp)

__call__(point)

Evaluate the vector field at a single point.

Parameters:

Name	Type	Description	Default
point	FloatArray	Position, shape (3,).	required

Returns:

Type	Description
FloatArray	Field vector, shape (3,). NaN if outside domain.

Source code in src/pypic/traces/_tracing.py

def __call__(self, point: FloatArray) -> FloatArray:
    """Evaluate the vector field at a single point.

    Parameters
    ----------
    point : FloatArray
        Position, shape ``(3,)``.

    Returns
    -------
    FloatArray
        Field vector, shape ``(3,)``. NaN if outside domain.
    """
    return self._interp(point.reshape(1, 3))[0]  # type: ignore[no-any-return]

batch(points)

Evaluate the vector field at N points in one call.

Used by the batched adaptive tracer so all seeds in a step share a single RegularGridInterpolator dispatch. Roughly N× faster than calling __call__ N times because per-call Python / argument-marshaling overhead amortizes over the batch.

Parameters:

Name	Type	Description	Default
points	FloatArray	Positions, shape (N, 3).	required

Returns:

Type	Description
FloatArray	Field vectors, shape (N, 3). Rows are NaN where the corresponding point is outside the interpolation domain.

Source code in src/pypic/traces/_tracing.py

def batch(self, points: FloatArray) -> FloatArray:
    """Evaluate the vector field at ``N`` points in one call.

    Used by the batched adaptive tracer so all seeds in a step share
    a single ``RegularGridInterpolator`` dispatch. Roughly N× faster
    than calling `__call__` N times because per-call Python /
    argument-marshaling overhead amortizes over the batch.

    Parameters
    ----------
    points : FloatArray
        Positions, shape ``(N, 3)``.

    Returns
    -------
    FloatArray
        Field vectors, shape ``(N, 3)``. Rows are NaN where the
        corresponding point is outside the interpolation domain.
    """
    return self._interp(points)

arc_length_cumulative(points)

Cumulative arc length along a curve.

\[s_i = \sum_{k=1}^{i} \|\mathbf{r}_k - \mathbf{r}_{k-1}\|\]

Parameters:

Name	Type	Description	Default
points	FloatArray	Ordered positions, shape (N, 3).	required

Returns:

Type	Description
FloatArray	Cumulative arc length, shape (N,). First element is 0.

Examples:

>>> import numpy as np
>>> pts = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]])
>>> arc_length_cumulative(pts)
array([0., 1., 2.])

Source code in src/pypic/traces/_analysis.py

def arc_length_cumulative(points: FloatArray) -> FloatArray:
    r"""Cumulative arc length along a curve.

    $$s_i = \sum_{k=1}^{i} \|\mathbf{r}_k - \mathbf{r}_{k-1}\|$$

    Parameters
    ----------
    points : FloatArray
        Ordered positions, shape ``(N, 3)``.

    Returns
    -------
    FloatArray
        Cumulative arc length, shape ``(N,)``. First element is 0.

    Examples
    --------
    >>> import numpy as np
    >>> pts = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]])
    >>> arc_length_cumulative(pts)
    array([0., 1., 2.])
    """
    segments = np.linalg.norm(np.diff(points, axis=0), axis=1)
    return np.concatenate(([0.0], np.cumsum(segments)))

arc_length_total(points)

Total arc length of a curve.

Parameters:

Name	Type	Description	Default
points	FloatArray	Ordered positions, shape (N, 3).	required

Returns:

Type	Description
float	Sum of segment lengths.

Examples:

>>> import numpy as np
>>> pts = np.array([[0.0, 0.0, 0.0], [3.0, 4.0, 0.0]])
>>> arc_length_total(pts)
5.0

Source code in src/pypic/traces/_analysis.py

def arc_length_total(points: FloatArray) -> float:
    r"""Total arc length of a curve.

    Parameters
    ----------
    points : FloatArray
        Ordered positions, shape ``(N, 3)``.

    Returns
    -------
    float
        Sum of segment lengths.

    Examples
    --------
    >>> import numpy as np
    >>> pts = np.array([[0.0, 0.0, 0.0], [3.0, 4.0, 0.0]])
    >>> arc_length_total(pts)
    5.0
    """
    return float(np.sum(np.linalg.norm(np.diff(points, axis=0), axis=1)))

closest_approach(points, target)

Find the point on the curve nearest to target.

Parameters:

Name	Type	Description	Default
points	FloatArray	Ordered positions, shape (N, 3).	required
target	Vector3	Reference point.	required

Returns:

Type	Description
tuple[int, float]	(index, distance) of the closest point.

Examples:

>>> import numpy as np
>>> pts = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]])
>>> closest_approach(pts, (0.9, 0.0, 0.0))
(1, 0.09999999999999998)

Source code in src/pypic/traces/_analysis.py

def closest_approach(points: FloatArray, target: Vector3) -> tuple[int, float]:
    """Find the point on the curve nearest to *target*.

    Parameters
    ----------
    points : FloatArray
        Ordered positions, shape ``(N, 3)``.
    target : Vector3
        Reference point.

    Returns
    -------
    tuple[int, float]
        ``(index, distance)`` of the closest point.

    Examples
    --------
    >>> import numpy as np
    >>> pts = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]])
    >>> closest_approach(pts, (0.9, 0.0, 0.0))
    (1, 0.09999999999999998)
    """
    distances = np.linalg.norm(points - np.asarray(target), axis=1)
    idx = int(np.argmin(distances))
    return idx, float(distances[idx])

curvature(points)

Curvature \(\kappa = \|d\hat{T}/ds\|\) along a curve.

Parameters:

Name	Type	Description	Default
points	FloatArray	Ordered positions, shape (N, 3).	required

Returns:

Type	Description
FloatArray	Curvature at each point, shape (N,).

Examples:

>>> import numpy as np
>>> pts = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]])
>>> np.testing.assert_allclose(curvature(pts), [0.0, 0.0, 0.0], atol=1e-15)

Source code in src/pypic/traces/_analysis.py

def curvature(points: FloatArray) -> FloatArray:
    r"""Curvature $\kappa = \|d\hat{T}/ds\|$ along a curve.

    Parameters
    ----------
    points : FloatArray
        Ordered positions, shape ``(N, 3)``.

    Returns
    -------
    FloatArray
        Curvature at each point, shape ``(N,)``.

    Examples
    --------
    >>> import numpy as np
    >>> pts = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]])
    >>> np.testing.assert_allclose(curvature(pts), [0.0, 0.0, 0.0], atol=1e-15)
    """
    dr = np.gradient(points, axis=0)
    ds = np.linalg.norm(dr, axis=1, keepdims=True)
    ds = np.maximum(ds, np.finfo(dr.dtype).tiny)
    t_hat = dr / ds
    dt = np.gradient(t_hat, axis=0)
    # dt/ds, where ds is the arc-length increment per index step
    return np.linalg.norm(dt / ds, axis=1)  # type: ignore[no-any-return]

displacement(points)

End-to-end displacement \(\|\mathbf{r}_N - \mathbf{r}_0\|\).

Parameters:

Name	Type	Description	Default
points	FloatArray	Ordered positions, shape (N, 3).	required

Returns:

Type	Description
float

Examples:

>>> import numpy as np
>>> pts = np.array([[0.0, 0.0, 0.0], [1.0, 1.0, 0.0], [3.0, 4.0, 0.0]])
>>> displacement(pts)
5.0

Source code in src/pypic/traces/_analysis.py

def displacement(points: FloatArray) -> float:
    r"""End-to-end displacement $\|\mathbf{r}_N - \mathbf{r}_0\|$.

    Parameters
    ----------
    points : FloatArray
        Ordered positions, shape ``(N, 3)``.

    Returns
    -------
    float

    Examples
    --------
    >>> import numpy as np
    >>> pts = np.array([[0.0, 0.0, 0.0], [1.0, 1.0, 0.0], [3.0, 4.0, 0.0]])
    >>> displacement(pts)
    5.0
    """
    return float(np.linalg.norm(points[-1] - points[0]))

drift_velocity(points, time, *, window=5)

Running-average guiding-center drift velocity.

Smooths the instantaneous velocity \(d\mathbf{r}/dt\) with a uniform window to approximate the guiding-center drift, filtering out the gyromotion.

Parameters:

Name	Type	Description	Default
points	FloatArray	Positions, shape (N, 3).	required
time	FloatArray	Time at each point, shape (N,).	required
window	int	Averaging window size (must be odd and >= 1).	5

Returns:

Type	Description
FloatArray	Smoothed velocity, shape (N, 3).

Examples:

>>> import numpy as np
>>> points = np.zeros((5, 3))
>>> points[:, 0] = np.arange(5.0)
>>> drift_velocity(points, np.arange(5.0), window=1)[0]
array([1., 0., 0.])

Source code in src/pypic/traces/_analysis.py

def drift_velocity(
    points: FloatArray,
    time: FloatArray,
    *,
    window: int = 5,
) -> FloatArray:
    r"""Running-average guiding-center drift velocity.

    Smooths the instantaneous velocity $d\mathbf{r}/dt$ with a uniform
    window to approximate the guiding-center drift, filtering out the
    gyromotion.

    Parameters
    ----------
    points : FloatArray
        Positions, shape ``(N, 3)``.
    time : FloatArray
        Time at each point, shape ``(N,)``.
    window : int
        Averaging window size (must be odd and >= 1).

    Returns
    -------
    FloatArray
        Smoothed velocity, shape ``(N, 3)``.

    Examples
    --------
    >>> import numpy as np
    >>> points = np.zeros((5, 3))
    >>> points[:, 0] = np.arange(5.0)
    >>> drift_velocity(points, np.arange(5.0), window=1)[0]
    array([1., 0., 0.])
    """
    if window < 1:
        msg = f"window must be >= 1, got {window}"
        raise ValueError(msg)
    dt = np.gradient(time)
    dt = np.maximum(dt, np.finfo(dt.dtype).tiny)
    v_inst = np.gradient(points, axis=0) / dt[:, np.newaxis]
    if window == 1:
        return v_inst  # type: ignore[no-any-return]
    kernel = np.ones(window) / window
    return np.column_stack(
        [np.convolve(v_inst[:, i], kernel, mode="same") for i in range(3)]
    )

equatorial_crossings(points)

Find positions where a curve crosses the \(z = 0\) plane.

Shortcut for plane_crossings(points, (0, 0, 1), 0.0).

Parameters:

Name	Type	Description	Default
points	FloatArray	Ordered positions, shape (N, 3).	required

Returns:

Type	Description
FloatArray	Crossing positions, shape (M, 3).

Examples:

>>> import numpy as np
>>> pts = np.array([[0.0, 0.0, -1.0], [0.0, 0.0, 1.0]])
>>> equatorial_crossings(pts)
array([[0., 0., 0.]])

Source code in src/pypic/traces/_analysis.py

def equatorial_crossings(points: FloatArray) -> FloatArray:
    r"""Find positions where a curve crosses the $z = 0$ plane.

    Shortcut for ``plane_crossings(points, (0, 0, 1), 0.0)``.

    Parameters
    ----------
    points : FloatArray
        Ordered positions, shape ``(N, 3)``.

    Returns
    -------
    FloatArray
        Crossing positions, shape ``(M, 3)``.

    Examples
    --------
    >>> import numpy as np
    >>> pts = np.array([[0.0, 0.0, -1.0], [0.0, 0.0, 1.0]])
    >>> equatorial_crossings(pts)
    array([[0., 0., 0.]])
    """
    return plane_crossings(points, (0.0, 0.0, 1.0), 0.0)

gyroradius_estimate(points, velocity, b_magnitude, charge, mass)

Estimate local gyroradius \(r_g = m v_\perp / (|q| B)\) along a path.

The perpendicular velocity is estimated by subtracting the field-aligned component using the local tangent direction as a proxy for the field direction.

Parameters:

Name	Type	Description	Default
points	FloatArray	Positions along the path, shape (N, 3).	required
velocity	FloatArray	Velocity at each point, shape (N, 3).	required
b_magnitude	FloatArray	Magnetic field magnitude at each point, shape (N,).	required
charge	float	Particle charge (absolute value used).	required
mass	float	Particle mass.	required

Returns:

Type	Description
FloatArray	Estimated gyroradius at each point, shape (N,). NaN where \(B = 0\).

Examples:

>>> import numpy as np
>>> points = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]])
>>> velocity = np.tile([0.0, 3.0, 0.0], (3, 1))
>>> b = np.full(3, 2.0)
>>> gyroradius_estimate(points, velocity, b, charge=1.0, mass=1.0)
array([1.5, 1.5, 1.5])

Source code in src/pypic/traces/_analysis.py

def gyroradius_estimate(
    points: FloatArray,
    velocity: FloatArray,
    b_magnitude: FloatArray,
    charge: float,
    mass: float,
) -> FloatArray:
    r"""Estimate local gyroradius $r_g = m v_\perp / (|q| B)$ along a path.

    The perpendicular velocity is estimated by subtracting the field-aligned
    component using the local tangent direction as a proxy for the field
    direction.

    Parameters
    ----------
    points : FloatArray
        Positions along the path, shape ``(N, 3)``.
    velocity : FloatArray
        Velocity at each point, shape ``(N, 3)``.
    b_magnitude : FloatArray
        Magnetic field magnitude at each point, shape ``(N,)``.
    charge : float
        Particle charge (absolute value used).
    mass : float
        Particle mass.

    Returns
    -------
    FloatArray
        Estimated gyroradius at each point, shape ``(N,)``.
        NaN where $B = 0$.

    Examples
    --------
    >>> import numpy as np
    >>> points = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]])
    >>> velocity = np.tile([0.0, 3.0, 0.0], (3, 1))
    >>> b = np.full(3, 2.0)
    >>> gyroradius_estimate(points, velocity, b, charge=1.0, mass=1.0)
    array([1.5, 1.5, 1.5])
    """
    t_hat = tangent_vectors(points)
    v_par = np.sum(velocity * t_hat, axis=1, keepdims=True) * t_hat
    v_perp_mag = np.linalg.norm(velocity - v_par, axis=1)
    abs_q = abs(charge)
    with np.errstate(divide="ignore", invalid="ignore"):
        rg = mass * v_perp_mag / (abs_q * b_magnitude)
    rg[b_magnitude == 0] = np.nan
    return rg  # type: ignore[no-any-return]

kinetic_energy(velocity, mass)

Non-relativistic kinetic energy \(\frac{1}{2} m v^2\) at each point.

Parameters:

Name	Type	Description	Default
velocity	FloatArray	Velocity vectors, shape (N, 3).	required
mass	float	Particle mass in code units.	required

Returns:

Type	Description
FloatArray	Kinetic energy at each point, shape (N,).

Examples:

>>> import numpy as np
>>> v = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0]])
>>> kinetic_energy(v, 2.0)
array([1., 4.])

Source code in src/pypic/traces/_analysis.py

def kinetic_energy(velocity: FloatArray, mass: float) -> FloatArray:
    r"""Non-relativistic kinetic energy $\frac{1}{2} m v^2$ at each point.

    Parameters
    ----------
    velocity : FloatArray
        Velocity vectors, shape ``(N, 3)``.
    mass : float
        Particle mass in code units.

    Returns
    -------
    FloatArray
        Kinetic energy at each point, shape ``(N,)``.

    Examples
    --------
    >>> import numpy as np
    >>> v = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0]])
    >>> kinetic_energy(v, 2.0)
    array([1., 4.])
    """
    return 0.5 * mass * np.sum(velocity**2, axis=1)  # type: ignore[no-any-return]

mirror_points(b_magnitude)

Find mirror points (local maxima in \(|\mathbf{B}|\)) along a path.

A mirror point is where a trapped particle reverses direction due to the magnetic mirror force, occurring at local maxima of the field magnitude along the particle's guiding-center path.

Parameters:

Name	Type	Description	Default
b_magnitude	FloatArray	Magnetic field magnitude along the path, shape (N,).	required

Returns:

Type	Description
NDArray[intp]	Integer indices of local maxima. Empty if none found.

Examples:

>>> import numpy as np
>>> b = np.array([1.0, 3.0, 1.0, 4.0, 2.0])
>>> mirror_points(b)
array([1, 3])

Source code in src/pypic/traces/_analysis.py

def mirror_points(b_magnitude: FloatArray) -> NDArray[np.intp]:
    r"""Find mirror points (local maxima in $|\mathbf{B}|$) along a path.

    A mirror point is where a trapped particle reverses direction due to
    the magnetic mirror force, occurring at local maxima of the field
    magnitude along the particle's guiding-center path.

    Parameters
    ----------
    b_magnitude : FloatArray
        Magnetic field magnitude along the path, shape ``(N,)``.

    Returns
    -------
    NDArray[np.intp]
        Integer indices of local maxima. Empty if none found.

    Examples
    --------
    >>> import numpy as np
    >>> b = np.array([1.0, 3.0, 1.0, 4.0, 2.0])
    >>> mirror_points(b)
    array([1, 3])
    """
    if len(b_magnitude) < 3:
        return np.array([], dtype=np.intp)
    left = b_magnitude[:-2]
    center = b_magnitude[1:-1]
    right = b_magnitude[2:]
    is_max = (center > left) & (center > right)
    return np.where(is_max)[0] + 1

plane_crossings(points, normal, offset=0.0)

Find positions where a curve crosses a plane.

The plane is defined by \(\mathbf{n} \cdot \mathbf{r} = d\) where \(\mathbf{n}\) is the normal and \(d\) is the offset. Crossing positions are linearly interpolated between consecutive points. Segments tangent to the plane (same signed distance at both endpoints) are not counted as crossings.

Parameters:

Name	Type	Description	Default
points	FloatArray	Ordered positions, shape (N, 3).	required
normal	Vector3	Plane normal vector (will be normalized internally).	required
offset	float	Signed distance from origin along the normal direction.	0.0

Returns:

Type	Description
FloatArray	Crossing positions, shape (M, 3). Empty (0, 3) if none.

Examples:

>>> import numpy as np
>>> pts = np.array([[-1.0, 0.0, 0.0], [1.0, 0.0, 0.0]])
>>> plane_crossings(pts, (1.0, 0.0, 0.0), 0.0)
array([[0., 0., 0.]])

Source code in src/pypic/traces/_analysis.py

def plane_crossings(
    points: FloatArray,
    normal: Vector3,
    offset: float = 0.0,
) -> FloatArray:
    r"""Find positions where a curve crosses a plane.

    The plane is defined by $\mathbf{n} \cdot \mathbf{r} = d$ where
    $\mathbf{n}$ is the normal and $d$ is the offset. Crossing positions
    are linearly interpolated between consecutive points. Segments
    tangent to the plane (same signed distance at both endpoints)
    are not counted as crossings.

    Parameters
    ----------
    points : FloatArray
        Ordered positions, shape ``(N, 3)``.
    normal : Vector3
        Plane normal vector (will be normalized internally).
    offset : float
        Signed distance from origin along the normal direction.

    Returns
    -------
    FloatArray
        Crossing positions, shape ``(M, 3)``. Empty ``(0, 3)`` if none.

    Examples
    --------
    >>> import numpy as np
    >>> pts = np.array([[-1.0, 0.0, 0.0], [1.0, 0.0, 0.0]])
    >>> plane_crossings(pts, (1.0, 0.0, 0.0), 0.0)
    array([[0., 0., 0.]])
    """
    n = np.asarray(normal, dtype=np.float64)
    n = n / np.linalg.norm(n)
    signed_dist = points @ n - offset
    crossings = []
    for i in range(len(signed_dist) - 1):
        d0, d1 = signed_dist[i], signed_dist[i + 1]
        if d0 * d1 < 0:
            t = d0 / (d0 - d1)
            crossings.append(points[i] + t * (points[i + 1] - points[i]))
    if not crossings:
        return np.empty((0, 3), dtype=points.dtype)
    return np.array(crossings)

resample_by_arc_length(points, n_out, *, scalars=None)

Resample a curve to uniform arc-length spacing via linear interpolation.

Parameters:

Name	Type	Description	Default
points	FloatArray	Ordered positions, shape (N, 3).	required
n_out	int	Number of output points (must be >= 2).	required
scalars	dict[str, FloatArray] | None	Optional scalar arrays to resample, each shape (N,).	None

Returns:

Type	Description
tuple[FloatArray, dict[str, FloatArray]]	(resampled_points, resampled_scalars) where resampled_points has shape (n_out, 3) and each scalar has shape (n_out,).

Examples:

>>> import numpy as np
>>> pts = np.array([[0.0, 0.0, 0.0], [2.0, 0.0, 0.0], [4.0, 0.0, 0.0]])
>>> new_pts, _ = resample_by_arc_length(pts, 5)
>>> new_pts[:, 0]
array([0., 1., 2., 3., 4.])

Source code in src/pypic/traces/_analysis.py

def resample_by_arc_length(
    points: FloatArray,
    n_out: int,
    *,
    scalars: dict[str, FloatArray] | None = None,
) -> tuple[FloatArray, dict[str, FloatArray]]:
    """Resample a curve to uniform arc-length spacing via linear interpolation.

    Parameters
    ----------
    points : FloatArray
        Ordered positions, shape ``(N, 3)``.
    n_out : int
        Number of output points (must be >= 2).
    scalars : dict[str, FloatArray] | None
        Optional scalar arrays to resample, each shape ``(N,)``.

    Returns
    -------
    tuple[FloatArray, dict[str, FloatArray]]
        ``(resampled_points, resampled_scalars)`` where ``resampled_points``
        has shape ``(n_out, 3)`` and each scalar has shape ``(n_out,)``.

    Examples
    --------
    >>> import numpy as np
    >>> pts = np.array([[0.0, 0.0, 0.0], [2.0, 0.0, 0.0], [4.0, 0.0, 0.0]])
    >>> new_pts, _ = resample_by_arc_length(pts, 5)
    >>> new_pts[:, 0]
    array([0., 1., 2., 3., 4.])
    """
    if n_out < 2:
        msg = f"n_out must be >= 2, got {n_out}"
        raise ValueError(msg)
    s = arc_length_cumulative(points)
    s_new = np.linspace(s[0], s[-1], n_out)
    resampled = np.column_stack([np.interp(s_new, s, points[:, i]) for i in range(3)])
    resampled_scalars: dict[str, FloatArray] = {}
    if scalars:
        for name, arr in scalars.items():
            resampled_scalars[name] = np.interp(s_new, s, arr)
    return resampled, resampled_scalars

speed(velocity)

Speed (magnitude of velocity) at each point.

Parameters:

Name	Type	Description	Default
velocity	FloatArray	Velocity vectors, shape (N, 3).	required

Returns:

Type	Description
FloatArray	Speed at each point, shape (N,).

Examples:

>>> import numpy as np
>>> v = np.array([[3.0, 4.0, 0.0], [0.0, 0.0, 5.0]])
>>> speed(v)
array([5., 5.])

Source code in src/pypic/traces/_analysis.py

def speed(velocity: FloatArray) -> FloatArray:
    r"""Speed (magnitude of velocity) at each point.

    Parameters
    ----------
    velocity : FloatArray
        Velocity vectors, shape ``(N, 3)``.

    Returns
    -------
    FloatArray
        Speed at each point, shape ``(N,)``.

    Examples
    --------
    >>> import numpy as np
    >>> v = np.array([[3.0, 4.0, 0.0], [0.0, 0.0, 5.0]])
    >>> speed(v)
    array([5., 5.])
    """
    return np.linalg.norm(velocity, axis=1)  # type: ignore[no-any-return]

tangent_vectors(points)

Compute unit tangent vectors along a curve via central differences.

Uses np.gradient for second-order central differences at interior points and one-sided differences at endpoints.

Parameters:

Name	Type	Description	Default
points	FloatArray	Ordered positions, shape (N, 3).	required

Returns:

Type	Description
FloatArray	Unit tangent vectors, shape (N, 3).

Examples:

>>> import numpy as np
>>> pts = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]])
>>> tangent_vectors(pts)
array([[1., 0., 0.],
       [1., 0., 0.],
       [1., 0., 0.]])

Source code in src/pypic/traces/_analysis.py

def tangent_vectors(points: FloatArray) -> FloatArray:
    r"""Compute unit tangent vectors along a curve via central differences.

    Uses ``np.gradient`` for second-order central differences at interior
    points and one-sided differences at endpoints.

    Parameters
    ----------
    points : FloatArray
        Ordered positions, shape ``(N, 3)``.

    Returns
    -------
    FloatArray
        Unit tangent vectors, shape ``(N, 3)``.

    Examples
    --------
    >>> import numpy as np
    >>> pts = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]])
    >>> tangent_vectors(pts)
    array([[1., 0., 0.],
           [1., 0., 0.],
           [1., 0., 0.]])
    """
    dr = np.gradient(points, axis=0)
    norms = np.linalg.norm(dr, axis=1, keepdims=True)
    norms = np.maximum(norms, np.finfo(dr.dtype).tiny)
    return dr / norms  # type: ignore[no-any-return]

poincare_section(data, seeds, surface, *, direction='forward', max_steps=50000, field_components=('B_1', 'B_2', 'B_3'), atol=1e-06, rtol=0.0001, step_size_init=0.5, min_step=1e-08, max_step=2.0, null_threshold=1e-12, interpolator=None)

Build a Poincaré section by tracing seeds and puncturing on \(\Sigma\).

Each seed is integrated forward (and/or backward) with the adaptive Dormand-Prince tracer; every accepted-step segment that crosses \(\Sigma\) contributes one puncture, located via linear interpolation on the signed-distance field \(\mathbf{n}\cdot\mathbf{r} - d\).

Internally forces loop_tol=None on the tracer — the auto closed-loop detector would otherwise terminate the very orbits we want to sample on the second puncture.

Parameters:

Name	Type	Description	Default
data	FieldDataset	Gridded vector field data.	required
seeds	FloatArray or Sequence[Vector3]	Seed positions, shape (M, 3) (or any iterable that np.asarray(..., dtype=float64).reshape(-1, 3) accepts).	required
surface	PoincareSurface	Transverse plane to puncture on.	required
direction	``"forward"`` | ``"backward"`` | ``"both"``	Trace direction. Forward only is typical for tokamak/stellarator cuts; "both" doubles puncture density for fixed max_steps.	'forward'
max_steps	int	Per-direction step budget for the adaptive tracer. Memory cost scales as M * max_steps * 3 * 8 bytes; 50_000 × 100 seeds is ~120 MB.	50000
field_components	tuple[str, str, str]	Field component names. Default ("B_1", "B_2", "B_3").	('B_1', 'B_2', 'B_3')
atol	float	Adaptive error tolerances.	1e-06
rtol	float	Adaptive error tolerances.	1e-06
step_size_init	float	Step-size controller knobs.	0.5
min_step	float	Step-size controller knobs.	0.5
max_step	float	Step-size controller knobs.	0.5
null_threshold	float	Field-magnitude threshold below which a point is a null.	1e-12
interpolator	VectorFieldInterpolator | None	Pre-built interpolator; constructed internally if None.	None

Returns:

Type	Description
PoincareSection

Raises:

Type	Description
ValueError	Propagated from trace_field_lines_adaptive (bad seed, invalid direction, bad tolerances).

Notes

For fusion poloidal sections, set direction="forward" and pick max_steps to span ~\(N\) poloidal transits per seed; N = 100 is typical for a quick island survey. For magnetotail X-line geometry, direction="both" captures the separatrix from both inflow regions.

Examples:

See PoincareSection for a closed-circle example.

Source code in src/pypic/traces/_poincare.py

def poincare_section(
    data: FieldDataset,
    seeds: FloatArray | Sequence[Vector3],
    surface: PoincareSurface,
    *,
    direction: TraceDirection = "forward",
    max_steps: int = 50_000,
    field_components: tuple[str, str, str] = ("B_1", "B_2", "B_3"),
    atol: float = 1e-6,
    rtol: float = 1e-4,
    step_size_init: float = 0.5,
    min_step: float = 1e-8,
    max_step: float = 2.0,
    null_threshold: float = 1e-12,
    interpolator: VectorFieldInterpolator | None = None,
) -> PoincareSection:
    r"""Build a Poincaré section by tracing seeds and puncturing on $\Sigma$.

    Each seed is integrated forward (and/or backward) with the adaptive
    Dormand-Prince tracer; every accepted-step segment that crosses
    $\Sigma$ contributes one puncture, located via linear interpolation
    on the signed-distance field $\mathbf{n}\cdot\mathbf{r} - d$.

    Internally forces ``loop_tol=None`` on the tracer — the auto
    closed-loop detector would otherwise terminate the very orbits we
    want to sample on the second puncture.

    Parameters
    ----------
    data : FieldDataset
        Gridded vector field data.
    seeds : FloatArray or Sequence[Vector3]
        Seed positions, shape ``(M, 3)`` (or any iterable that
        ``np.asarray(..., dtype=float64).reshape(-1, 3)`` accepts).
    surface : PoincareSurface
        Transverse plane to puncture on.
    direction : ``"forward"`` | ``"backward"`` | ``"both"``
        Trace direction. Forward only is typical for tokamak/stellarator
        cuts; ``"both"`` doubles puncture density for fixed ``max_steps``.
    max_steps : int
        Per-direction step budget for the adaptive tracer. Memory cost
        scales as ``M * max_steps * 3 * 8`` bytes; 50_000 × 100 seeds
        is ~120 MB.
    field_components : tuple[str, str, str]
        Field component names. Default ``("B_1", "B_2", "B_3")``.
    atol, rtol : float
        Adaptive error tolerances.
    step_size_init, min_step, max_step : float
        Step-size controller knobs.
    null_threshold : float
        Field-magnitude threshold below which a point is a null.
    interpolator : VectorFieldInterpolator | None
        Pre-built interpolator; constructed internally if ``None``.

    Returns
    -------
    PoincareSection

    Raises
    ------
    ValueError
        Propagated from `trace_field_lines_adaptive` (bad seed,
        invalid direction, bad tolerances).

    Notes
    -----
    For fusion poloidal sections, set ``direction="forward"`` and pick
    ``max_steps`` to span ~$N$ poloidal transits per seed; ``N = 100``
    is typical for a quick island survey. For magnetotail X-line
    geometry, ``direction="both"`` captures the separatrix from both
    inflow regions.

    Examples
    --------
    See `PoincareSection` for a closed-circle example.
    """
    seeds_arr = np.asarray(seeds, dtype=np.float64).reshape(-1, 3)

    field_lines = trace_field_lines_adaptive(
        data,
        seeds_arr,
        atol=atol,
        rtol=rtol,
        step_size_init=step_size_init,
        min_step=min_step,
        max_step=max_step,
        max_steps=max_steps,
        direction=direction,
        field_components=field_components,
        null_threshold=null_threshold,
        loop_tol=None,
        interpolator=interpolator,
    )

    offset = surface.offset

    punctures_3d_list: list[FloatArray] = []
    punctures_2d_list: list[FloatArray] = []
    for fl in field_lines:
        pts3 = plane_crossings(fl.points, surface.normal, offset)
        punctures_3d_list.append(pts3)
        punctures_2d_list.append(surface.project(pts3))

    return PoincareSection(
        surface=surface,
        seeds=seeds_arr,
        direction=direction,
        punctures_3d=tuple(punctures_3d_list),
        punctures_2d=tuple(punctures_2d_list),
        field_lines=tuple(field_lines),
        metadata={
            "n_crossings_per_seed": np.array(
                [p.shape[0] for p in punctures_3d_list], dtype=np.intp
            ),
            "n_steps_per_seed": np.array(
                [fl.metadata.get("n_steps", fl.n_points - 1) for fl in field_lines],
                dtype=np.intp,
            ),
            "termination_reasons": tuple(
                fl.metadata.get("reason", "") for fl in field_lines
            ),
            "surface_name": surface.name,
        },
    )

attach_scalars(field_line, data, fields, *, method='nearest')

Return a new FieldLine with sampled scalars merged in.

Parameters:

Name	Type	Description	Default
field_line	FieldLine	Input field line.	required
data	FieldDataset	Gridded data to sample from.	required
fields	list[str]	Field names to sample and attach.	required
method	str	Interpolation method.	'nearest'

Returns:

Type	Description
FieldLine	New instance with additional scalars.

Source code in src/pypic/traces/_sampling.py

def attach_scalars(
    field_line: FieldLine,
    data: FieldDataset,
    fields: list[str],
    *,
    method: str = "nearest",
) -> FieldLine:
    """Return a new FieldLine with sampled scalars merged in.

    Parameters
    ----------
    field_line : FieldLine
        Input field line.
    data : FieldDataset
        Gridded data to sample from.
    fields : list[str]
        Field names to sample and attach.
    method : str
        Interpolation method.

    Returns
    -------
    FieldLine
        New instance with additional scalars.
    """
    sampled = sample_fields(data, field_line.points, fields, method=method)
    merged = dict(field_line.scalars)
    merged.update(sampled)
    return copy.replace(field_line, scalars=merged)

attach_scalars_to_trace(trace, data, fields, *, method='nearest')

Return a new ParticleTrace with sampled scalars merged in.

Parameters:

Name	Type	Description	Default
trace	ParticleTrace	Input particle trace.	required
data	FieldDataset	Gridded data to sample from.	required
fields	list[str]	Field names to sample and attach.	required
method	str	Interpolation method.	'nearest'

Returns:

Type	Description
ParticleTrace	New instance with additional scalars.

Source code in src/pypic/traces/_sampling.py

def attach_scalars_to_trace(
    trace: ParticleTrace,
    data: FieldDataset,
    fields: list[str],
    *,
    method: str = "nearest",
) -> ParticleTrace:
    """Return a new ParticleTrace with sampled scalars merged in.

    Parameters
    ----------
    trace : ParticleTrace
        Input particle trace.
    data : FieldDataset
        Gridded data to sample from.
    fields : list[str]
        Field names to sample and attach.
    method : str
        Interpolation method.

    Returns
    -------
    ParticleTrace
        New instance with additional scalars.
    """
    sampled = sample_fields(data, trace.points, fields, method=method)
    merged = dict(trace.scalars)
    merged.update(sampled)
    return copy.replace(trace, scalars=merged)

sample_field(data, points, field, *, method='nearest')

Sample a scalar field from a FieldDataset at arbitrary 3D positions.

Parameters:

Name	Type	Description	Default
data	FieldDataset	Gridded field data.	required
points	FloatArray	Query positions, shape (N, 3).	required
field	str	Field name (canonical, alias, or computable via compute()).	required
method	str	Interpolation method: "nearest" (pure NumPy) or "linear" (uses scipy.interpolate.RegularGridInterpolator).	'nearest'

Returns:

Type	Description
FloatArray	Sampled values, shape (N,). NaN for points outside the domain.

Examples:

>>> import numpy as np
>>> from pypic import FieldDataset, GridInfo, Normalization
>>> grid = GridInfo(dimensions=(4, 4, 4), spacing=(1.0, 1.0, 1.0))
>>> ds = FieldDataset.from_arrays(
...     {"B_1": np.arange(64.0).reshape(4, 4, 4)},
...     grid, Normalization.identity(),
... )
>>> points = np.array([[0.5, 0.5, 0.5], [1.5, 0.5, 0.5], [99.0, 0.0, 0.0]])
>>> sample_field(ds, points, "B_1")  # last point is outside the domain
array([ 0., 16., nan])

Source code in src/pypic/traces/_sampling.py

def sample_field(
    data: FieldDataset,
    points: FloatArray,
    field: str,
    *,
    method: str = "nearest",
) -> FloatArray:
    """Sample a scalar field from a FieldDataset at arbitrary 3D positions.

    Parameters
    ----------
    data : FieldDataset
        Gridded field data.
    points : FloatArray
        Query positions, shape ``(N, 3)``.
    field : str
        Field name (canonical, alias, or computable via ``compute()``).
    method : str
        Interpolation method: ``"nearest"`` (pure NumPy) or ``"linear"``
        (uses ``scipy.interpolate.RegularGridInterpolator``).

    Returns
    -------
    FloatArray
        Sampled values, shape ``(N,)``. NaN for points outside the domain.

    Examples
    --------
    >>> import numpy as np
    >>> from pypic import FieldDataset, GridInfo, Normalization
    >>> grid = GridInfo(dimensions=(4, 4, 4), spacing=(1.0, 1.0, 1.0))
    >>> ds = FieldDataset.from_arrays(
    ...     {"B_1": np.arange(64.0).reshape(4, 4, 4)},
    ...     grid, Normalization.identity(),
    ... )
    >>> points = np.array([[0.5, 0.5, 0.5], [1.5, 0.5, 0.5], [99.0, 0.0, 0.0]])
    >>> sample_field(ds, points, "B_1")  # last point is outside the domain
    array([ 0., 16., nan])
    """
    values: np.ndarray[Any, Any] = data[field]
    coord_arrays = data.grid.coordinate_arrays()

    if method == "nearest":
        return _sample_nearest(coord_arrays, values, points)
    if method == "linear":
        return _sample_linear(coord_arrays, values, points)
    msg = f"Unknown interpolation method {method!r}. Use 'nearest' or 'linear'."
    raise ValueError(msg)

sample_fields(data, points, fields, *, method='nearest')

Sample multiple fields at the same positions.

The nearest-index lookup (method="nearest") or interpolator construction (method="linear") is performed once and reused across all requested fields — the common case for trace workflows that sample B_1/B_2/B_3 or several diagnostics at the same points.

Parameters:

Name	Type	Description	Default
data	FieldDataset	Gridded field data.	required
points	FloatArray	Query positions, shape (N, 3).	required
fields	list[str]	Field names to sample.	required
method	str	Interpolation method: "nearest" or "linear".	'nearest'

Returns:

Type	Description
dict[str, FloatArray]	Field name → sampled values.

Examples:

>>> import numpy as np
>>> from pypic import FieldDataset, GridInfo, Normalization
>>> grid = GridInfo(dimensions=(4, 4, 4), spacing=(1.0, 1.0, 1.0))
>>> ds = FieldDataset.from_arrays(
...     {"B_1": np.arange(64.0).reshape(4, 4, 4),
...      "B_2": np.zeros((4, 4, 4))},
...     grid, Normalization.identity(),
... )
>>> points = np.array([[0.5, 0.5, 0.5], [1.5, 0.5, 0.5]])
>>> sampled = sample_fields(ds, points, ["B_1", "B_2"])
>>> sampled["B_1"], sampled["B_2"]
(array([ 0., 16.]), array([0., 0.]))

Source code in src/pypic/traces/_sampling.py

def sample_fields(
    data: FieldDataset,
    points: FloatArray,
    fields: list[str],
    *,
    method: str = "nearest",
) -> dict[str, FloatArray]:
    """Sample multiple fields at the same positions.

    The nearest-index lookup (``method="nearest"``) or interpolator
    construction (``method="linear"``) is performed once and reused
    across all requested fields — the common case for trace workflows
    that sample ``B_1/B_2/B_3`` or several diagnostics at the same points.

    Parameters
    ----------
    data : FieldDataset
        Gridded field data.
    points : FloatArray
        Query positions, shape ``(N, 3)``.
    fields : list[str]
        Field names to sample.
    method : str
        Interpolation method: ``"nearest"`` or ``"linear"``.

    Returns
    -------
    dict[str, FloatArray]
        Field name → sampled values.

    Examples
    --------
    >>> import numpy as np
    >>> from pypic import FieldDataset, GridInfo, Normalization
    >>> grid = GridInfo(dimensions=(4, 4, 4), spacing=(1.0, 1.0, 1.0))
    >>> ds = FieldDataset.from_arrays(
    ...     {"B_1": np.arange(64.0).reshape(4, 4, 4),
    ...      "B_2": np.zeros((4, 4, 4))},
    ...     grid, Normalization.identity(),
    ... )
    >>> points = np.array([[0.5, 0.5, 0.5], [1.5, 0.5, 0.5]])
    >>> sampled = sample_fields(ds, points, ["B_1", "B_2"])
    >>> sampled["B_1"], sampled["B_2"]
    (array([ 0., 16.]), array([0., 0.]))
    """
    if not fields:
        return {}

    coord_arrays = data.grid.coordinate_arrays()
    arrays = {name: np.asarray(data[name]) for name in fields}

    if method == "nearest":
        indices, mask = _nearest_indices(coord_arrays, points)
        idx_tuple = tuple(indices[mask, d] for d in range(len(coord_arrays)))
        out: dict[str, FloatArray] = {}
        for name, values in arrays.items():
            result = np.full(points.shape[0], np.nan, dtype=np.float64)
            if np.any(mask):
                result[mask] = values[idx_tuple]
            out[name] = result
        return out

    if method == "linear":
        from scipy.interpolate import RegularGridInterpolator

        stacked = np.stack(list(arrays.values()), axis=-1)
        interp = RegularGridInterpolator(
            coord_arrays,
            stacked,
            method="linear",
            bounds_error=False,
            fill_value=np.nan,
        )
        sampled = interp(points[:, : len(coord_arrays)])
        return {name: sampled[..., i] for i, name in enumerate(arrays)}

    msg = f"Unknown interpolation method {method!r}. Use 'nearest' or 'linear'."
    raise ValueError(msg)

estimate_tracing_error(field_line, data, *, field_components=('B_1', 'B_2', 'B_3'), interpolator=None)

Estimate tracing error via Richardson extrapolation.

Re-traces at half the original step size and compares endpoints. For an order-\(p\) scheme, halving the step reduces truncation error by \(2^p\); a ratio of ~16 confirms RK4's 4th-order convergence.

Only works with fixed-step field lines (from trace_field_line). Adaptive traces store max_local_error in metadata instead.

Parameters:

Name	Type	Description	Default
field_line	FieldLine	Previously traced field line (needs step_size in metadata).	required
data	FieldDataset	Same data used for the original trace.	required
field_components	tuple[str, str, str]	Vector field component names.	('B_1', 'B_2', 'B_3')
interpolator	VectorFieldInterpolator | None	Pre-built interpolator. Built internally if None.	None

Returns:

Type	Description
float	L2 distance between original and refined endpoints.

Examples:

>>> import numpy as np
>>> from pypic.dataset import FieldDataset
>>> from pypic.grid import GridInfo
>>> from pypic.units import Normalization
>>> grid = GridInfo(dimensions=(8, 8, 8), spacing=(1.0, 1.0, 1.0))
>>> data = FieldDataset.from_arrays(
...     {
...         "B_1": np.ones((8, 8, 8)),
...         "B_2": np.zeros((8, 8, 8)),
...         "B_3": np.zeros((8, 8, 8)),
...     },
...     grid,
...     Normalization.identity(),
... )
>>> fl = trace_field_line(
...     data, (4.0, 4.0, 4.0), step_size=1.0, max_steps=3, direction="forward"
... )
>>> err = estimate_tracing_error(fl, data)
>>> err < 1e-9  # uniform field, RK4 is exact
True

Source code in src/pypic/traces/_tracing.py

def estimate_tracing_error(
    field_line: FieldLine,
    data: FieldDataset,
    *,
    field_components: tuple[str, str, str] = ("B_1", "B_2", "B_3"),
    interpolator: VectorFieldInterpolator | None = None,
) -> float:
    r"""Estimate tracing error via Richardson extrapolation.

    Re-traces at half the original step size and compares endpoints.
    For an order-$p$ scheme, halving the step reduces truncation error
    by $2^p$; a ratio of ~16 confirms RK4's 4th-order convergence.

    Only works with fixed-step field lines (from ``trace_field_line``).
    Adaptive traces store ``max_local_error`` in metadata instead.

    Parameters
    ----------
    field_line : FieldLine
        Previously traced field line (needs ``step_size`` in metadata).
    data : FieldDataset
        Same data used for the original trace.
    field_components : tuple[str, str, str]
        Vector field component names.
    interpolator : VectorFieldInterpolator | None
        Pre-built interpolator. Built internally if ``None``.

    Returns
    -------
    float
        L2 distance between original and refined endpoints.

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.dataset import FieldDataset
    >>> from pypic.grid import GridInfo
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(dimensions=(8, 8, 8), spacing=(1.0, 1.0, 1.0))
    >>> data = FieldDataset.from_arrays(
    ...     {
    ...         "B_1": np.ones((8, 8, 8)),
    ...         "B_2": np.zeros((8, 8, 8)),
    ...         "B_3": np.zeros((8, 8, 8)),
    ...     },
    ...     grid,
    ...     Normalization.identity(),
    ... )
    >>> fl = trace_field_line(
    ...     data, (4.0, 4.0, 4.0), step_size=1.0, max_steps=3, direction="forward"
    ... )
    >>> err = estimate_tracing_error(fl, data)
    >>> err < 1e-9  # uniform field, RK4 is exact
    True
    """
    for key in ("step_size", "n_steps"):
        if key not in field_line.metadata:
            msg = (
                f"FieldLine metadata missing required key {key!r}. "
                f"Only fixed-step traces (trace_field_line) support error estimation."
            )
            raise ValueError(msg)

    original_step = field_line.metadata["step_size"]
    n_steps = field_line.metadata["n_steps"]

    refined = trace_field_line(
        data,
        field_line.seed_point,
        step_size=original_step / 2.0,
        max_steps=n_steps * 2,
        direction=field_line.direction,
        field_components=field_components,
        interpolator=interpolator,
    )
    orig_end = np.array(field_line.end_point)
    ref_end = np.array(refined.end_point)
    return float(np.linalg.norm(orig_end - ref_end))

trace_field_line(data, seed, *, step_size=0.5, max_steps=10000, direction='both', field_components=('B_1', 'B_2', 'B_3'), null_threshold=1e-12, terminate=None, interpolator=None)

Trace a field line using classical (fixed-step) RK4.

Integrates \(d\mathbf{r}/ds = \hat{\mathbf{B}}(\mathbf{r})\) where \(\hat{\mathbf{B}} = \mathbf{B}/|\mathbf{B}|\) and \(s\) is arc length.

Parameters:

Name	Type	Description	Default
data	FieldDataset	Gridded vector field data.	required
seed	Vector3	Starting point (x, y, z) for the trace.	required
step_size	float	Arc-length step size in code units.	0.5
max_steps	int	Maximum integration steps per direction.	10000
direction	str	"forward", "backward", or "both".	'both'
field_components	tuple[str, str, str]	Names of the three vector field components.	('B_1', 'B_2', 'B_3')
null_threshold	float	Field magnitude below which the point is a null.	1e-12
terminate	Callable[[FloatArray], bool] | None	Optional callback; stops if it returns True.	None
interpolator	VectorFieldInterpolator | None	Pre-built interpolator. Built internally if None.	None

Returns:

Type	Description
FieldLine

Raises:

Type	Description
ValueError	If seed is outside domain or at a null point.

Examples:

>>> import numpy as np
>>> from pypic.dataset import FieldDataset
>>> from pypic.grid import GridInfo
>>> from pypic.units import Normalization
>>> grid = GridInfo(dimensions=(8, 8, 8), spacing=(1.0, 1.0, 1.0))
>>> data = FieldDataset.from_arrays(
...     {
...         "B_1": np.ones((8, 8, 8)),
...         "B_2": np.zeros((8, 8, 8)),
...         "B_3": np.zeros((8, 8, 8)),
...     },
...     grid,
...     Normalization.identity(),
... )
>>> fl = trace_field_line(
...     data, (4.0, 4.0, 4.0), step_size=0.5, max_steps=4, direction="forward"
... )
>>> fl.n_points
5
>>> bool(fl.points[-1, 0] > fl.points[0, 0])  # advances along +x
True

Source code in src/pypic/traces/_tracing.py

def trace_field_line(
    data: FieldDataset,
    seed: Vector3,
    *,
    step_size: float = 0.5,
    max_steps: int = 10_000,
    direction: TraceDirection = "both",
    field_components: tuple[str, str, str] = ("B_1", "B_2", "B_3"),
    null_threshold: float = 1e-12,
    terminate: Callable[[FloatArray], bool] | None = None,
    interpolator: VectorFieldInterpolator | None = None,
) -> FieldLine:
    r"""Trace a field line using classical (fixed-step) RK4.

    Integrates $d\mathbf{r}/ds = \hat{\mathbf{B}}(\mathbf{r})$ where
    $\hat{\mathbf{B}} = \mathbf{B}/|\mathbf{B}|$ and $s$ is arc length.

    Parameters
    ----------
    data : FieldDataset
        Gridded vector field data.
    seed : Vector3
        Starting point ``(x, y, z)`` for the trace.
    step_size : float
        Arc-length step size in code units.
    max_steps : int
        Maximum integration steps per direction.
    direction : str
        ``"forward"``, ``"backward"``, or ``"both"``.
    field_components : tuple[str, str, str]
        Names of the three vector field components.
    null_threshold : float
        Field magnitude below which the point is a null.
    terminate : Callable[[FloatArray], bool] | None
        Optional callback; stops if it returns ``True``.
    interpolator : VectorFieldInterpolator | None
        Pre-built interpolator. Built internally if ``None``.

    Returns
    -------
    FieldLine

    Raises
    ------
    ValueError
        If seed is outside domain or at a null point.

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.dataset import FieldDataset
    >>> from pypic.grid import GridInfo
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(dimensions=(8, 8, 8), spacing=(1.0, 1.0, 1.0))
    >>> data = FieldDataset.from_arrays(
    ...     {
    ...         "B_1": np.ones((8, 8, 8)),
    ...         "B_2": np.zeros((8, 8, 8)),
    ...         "B_3": np.zeros((8, 8, 8)),
    ...     },
    ...     grid,
    ...     Normalization.identity(),
    ... )
    >>> fl = trace_field_line(
    ...     data, (4.0, 4.0, 4.0), step_size=0.5, max_steps=4, direction="forward"
    ... )
    >>> fl.n_points
    5
    >>> bool(fl.points[-1, 0] > fl.points[0, 0])  # advances along +x
    True
    """
    if direction not in _VALID_DIRECTIONS:
        msg = f"direction must be one of {sorted(_VALID_DIRECTIONS)}, got {direction!r}"
        raise ValueError(msg)

    if interpolator is None:
        interpolator = VectorFieldInterpolator.from_dataset(data, field_components)

    seed_arr = np.asarray(seed, dtype=np.float64)
    _validate_seed(seed_arr, interpolator, null_threshold)

    field_name = _field_name_from_components(field_components)
    meta: dict[str, Any] = {"step_size": step_size, "method": "rk4"}
    args = (step_size, max_steps, null_threshold, terminate)

    match direction:
        case "forward":
            fwd, fwd_r = _trace_single_direction(interpolator, seed_arr, 1.0, *args)
            return _assemble_field_line(
                fwd,
                fwd_r,
                _EMPTY_POINTS,
                TerminationReason.MAX_STEPS,
                seed,
                direction,
                field_name,
                data.normalization,
                meta,
            )
        case "backward":
            bwd, bwd_r = _trace_single_direction(interpolator, seed_arr, -1.0, *args)
            return _assemble_field_line(
                _EMPTY_POINTS,
                TerminationReason.MAX_STEPS,
                bwd,
                bwd_r,
                seed,
                direction,
                field_name,
                data.normalization,
                meta,
            )
        case "both":
            fwd, fwd_r = _trace_single_direction(interpolator, seed_arr, 1.0, *args)
            bwd, bwd_r = _trace_single_direction(interpolator, seed_arr, -1.0, *args)
            return _assemble_field_line(
                fwd,
                fwd_r,
                bwd,
                bwd_r,
                seed,
                direction,
                field_name,
                data.normalization,
                meta,
            )
        case _ as unreachable:
            assert_never(unreachable)

trace_field_line_adaptive(data, seed, *, atol=1e-06, rtol=0.001, step_size_init=0.5, min_step=1e-08, max_step=2.0, max_steps=10000, direction='both', field_components=('B_1', 'B_2', 'B_3'), null_threshold=1e-12, terminate=None, loop_tol='auto', loop_min_arclen=None, interpolator=None)

Trace a field line using adaptive Dormand-Prince 5(4).

Uses embedded error estimation to adapt step size, taking larger steps in smooth regions and smaller steps near strong curvature.

Parameters:

Name	Type	Description	Default
data	FieldDataset	Gridded vector field data.	required
seed	Vector3	Starting point (x, y, z) for the trace.	required
atol	float	Absolute error tolerance.	1e-06
rtol	float	Relative error tolerance.	0.001
step_size_init	float	Initial arc-length step size.	0.5
min_step	float	Minimum allowed step size.	1e-08
max_step	float	Maximum allowed step size.	2.0
max_steps	int	Maximum accepted steps per direction.	10000
direction	str	"forward", "backward", or "both".	'both'
field_components	tuple[str, str, str]	Names of the three vector field components.	('B_1', 'B_2', 'B_3')
null_threshold	float	Field magnitude below which the point is a null.	1e-12
terminate	Callable[[FloatArray], bool] | None	Optional callback; stops if it returns True.	None
loop_tol	float, None, or ``"auto"``	Proximity threshold for closed-loop detection, in code units. The trace terminates with TerminationReason.CLOSED_LOOP as soon as it re-enters a loop_tol-radius ball around any previously visited point separated by more than loop_min_arclen of arc length. The canonical use case is mirror-mode magnetic holes and O-type islands [@Ahmadi2024], where an open RK tracer would otherwise burn through max_steps on a single closed orbit. A sliding-window proximity check is the streaming termination criterion used by streamline-visualization tooling; the Poincaré-map / invariant-manifold approach used by fusion boundary codes [@Frerichs2024] is for systematic island characterization, not per-trace short-circuiting. Default "auto" derives loop_tol = 0.5 * min(data.grid.spacing) — half a cell along the tightest grid axis — so closed orbits terminate cleanly even when the user didn't anticipate them. Pass None to force-disable detection (recovers pre-v0.X behavior), or pass a float to override the auto value.	'auto'
loop_min_arclen	float or None	Minimum arc-length distance between the current point and a candidate past point before a proximity hit counts as a closed loop. Prevents self-trigger on the immediately preceding samples when loop_tol is comparable to the step size. When None and loop_tol is set, defaults to 10 * step_size_init. Has no effect when loop_tol is None.	None
interpolator	VectorFieldInterpolator | None	Pre-built interpolator. Built internally if None.	None

Returns:

Type	Description
FieldLine

Raises:

Type	Description
ValueError	If seed is outside domain or at a null point, if loop_min_arclen is set without loop_tol, or if either loop threshold is non-positive.

Examples:

>>> import numpy as np
>>> from pypic.dataset import FieldDataset
>>> from pypic.grid import GridInfo
>>> from pypic.units import Normalization
>>> grid = GridInfo(dimensions=(8, 8, 8), spacing=(1.0, 1.0, 1.0))
>>> data = FieldDataset.from_arrays(
...     {
...         "B_1": np.ones((8, 8, 8)),
...         "B_2": np.zeros((8, 8, 8)),
...         "B_3": np.zeros((8, 8, 8)),
...     },
...     grid,
...     Normalization.identity(),
... )
>>> fl = trace_field_line_adaptive(
...     data, (4.0, 4.0, 4.0), max_steps=4, direction="forward"
... )
>>> fl.metadata["method"]
'rk45_dopri'
>>> "max_local_error" in fl.metadata
True

Source code in src/pypic/traces/_tracing.py

def trace_field_line_adaptive(
    data: FieldDataset,
    seed: Vector3,
    *,
    atol: float = 1e-6,
    rtol: float = 1e-3,
    step_size_init: float = 0.5,
    min_step: float = 1e-8,
    max_step: float = 2.0,
    max_steps: int = 10_000,
    direction: TraceDirection = "both",
    field_components: tuple[str, str, str] = ("B_1", "B_2", "B_3"),
    null_threshold: float = 1e-12,
    terminate: Callable[[FloatArray], bool] | None = None,
    loop_tol: float | None | Literal["auto"] = "auto",
    loop_min_arclen: float | None = None,
    interpolator: VectorFieldInterpolator | None = None,
) -> FieldLine:
    r"""Trace a field line using adaptive Dormand-Prince 5(4).

    Uses embedded error estimation to adapt step size, taking larger
    steps in smooth regions and smaller steps near strong curvature.

    Parameters
    ----------
    data : FieldDataset
        Gridded vector field data.
    seed : Vector3
        Starting point ``(x, y, z)`` for the trace.
    atol : float
        Absolute error tolerance.
    rtol : float
        Relative error tolerance.
    step_size_init : float
        Initial arc-length step size.
    min_step : float
        Minimum allowed step size.
    max_step : float
        Maximum allowed step size.
    max_steps : int
        Maximum accepted steps per direction.
    direction : str
        ``"forward"``, ``"backward"``, or ``"both"``.
    field_components : tuple[str, str, str]
        Names of the three vector field components.
    null_threshold : float
        Field magnitude below which the point is a null.
    terminate : Callable[[FloatArray], bool] | None
        Optional callback; stops if it returns ``True``.
    loop_tol : float, None, or ``"auto"``
        Proximity threshold for closed-loop detection, in code units.
        The trace terminates with
        `TerminationReason.CLOSED_LOOP` as soon as it re-enters a
        ``loop_tol``-radius ball around any previously visited point
        separated by more than ``loop_min_arclen`` of arc length. The
        canonical use case is mirror-mode magnetic holes and O-type
        islands [@Ahmadi2024], where an open RK tracer would otherwise
        burn through ``max_steps`` on a single closed orbit. A
        sliding-window proximity check is the streaming termination
        criterion used by streamline-visualization tooling; the
        Poincaré-map / invariant-manifold approach used by fusion
        boundary codes [@Frerichs2024] is for systematic island
        characterization, not per-trace short-circuiting. Default
        ``"auto"`` derives ``loop_tol = 0.5 * min(data.grid.spacing)``
        — half a cell along the tightest grid axis — so closed orbits
        terminate cleanly even when the user didn't anticipate them.
        Pass ``None`` to force-disable detection (recovers pre-v0.X
        behavior), or pass a float to override the auto value.
    loop_min_arclen : float or None
        Minimum arc-length distance between the current point and a
        candidate past point before a proximity hit counts as a closed
        loop. Prevents self-trigger on the immediately preceding samples
        when ``loop_tol`` is comparable to the step size. When ``None``
        and ``loop_tol`` is set, defaults to ``10 * step_size_init``.
        Has no effect when ``loop_tol`` is ``None``.
    interpolator : VectorFieldInterpolator | None
        Pre-built interpolator. Built internally if ``None``.

    Returns
    -------
    FieldLine

    Raises
    ------
    ValueError
        If seed is outside domain or at a null point, if
        ``loop_min_arclen`` is set without ``loop_tol``, or if either
        loop threshold is non-positive.

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.dataset import FieldDataset
    >>> from pypic.grid import GridInfo
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(dimensions=(8, 8, 8), spacing=(1.0, 1.0, 1.0))
    >>> data = FieldDataset.from_arrays(
    ...     {
    ...         "B_1": np.ones((8, 8, 8)),
    ...         "B_2": np.zeros((8, 8, 8)),
    ...         "B_3": np.zeros((8, 8, 8)),
    ...     },
    ...     grid,
    ...     Normalization.identity(),
    ... )
    >>> fl = trace_field_line_adaptive(
    ...     data, (4.0, 4.0, 4.0), max_steps=4, direction="forward"
    ... )
    >>> fl.metadata["method"]
    'rk45_dopri'
    >>> "max_local_error" in fl.metadata
    True
    """
    if direction not in _VALID_DIRECTIONS:
        msg = f"direction must be one of {sorted(_VALID_DIRECTIONS)}, got {direction!r}"
        raise ValueError(msg)

    if loop_tol == "auto":
        loop_tol = _auto_loop_tol(data)
    loop_tol, loop_min_arclen = _resolve_loop_kwargs(
        loop_tol, loop_min_arclen, step_size_init
    )

    if interpolator is None:
        interpolator = VectorFieldInterpolator.from_dataset(data, field_components)

    seed_arr = np.asarray(seed, dtype=np.float64)
    _validate_seed(seed_arr, interpolator, null_threshold)

    field_name = _field_name_from_components(field_components)

    def _adapt(sign: float) -> tuple[FloatArray, TerminationReason, float]:
        return _trace_single_direction_adaptive(
            interpolator,
            seed_arr,
            sign,
            atol,
            rtol,
            step_size_init,
            min_step,
            max_step,
            max_steps,
            null_threshold,
            terminate,
            loop_tol,
            loop_min_arclen,
        )

    def _meta(err: float) -> dict[str, Any]:
        return {
            "method": "rk45_dopri",
            "atol": atol,
            "rtol": rtol,
            "max_local_error": err,
        }

    match direction:
        case "forward":
            fwd, fwd_r, fwd_err = _adapt(1.0)
            return _assemble_field_line(
                fwd,
                fwd_r,
                _EMPTY_POINTS,
                TerminationReason.MAX_STEPS,
                seed,
                direction,
                field_name,
                data.normalization,
                _meta(fwd_err),
            )
        case "backward":
            bwd, bwd_r, bwd_err = _adapt(-1.0)
            return _assemble_field_line(
                _EMPTY_POINTS,
                TerminationReason.MAX_STEPS,
                bwd,
                bwd_r,
                seed,
                direction,
                field_name,
                data.normalization,
                _meta(bwd_err),
            )
        case "both":
            fwd, fwd_r, fwd_err = _adapt(1.0)
            bwd, bwd_r, bwd_err = _adapt(-1.0)
            return _assemble_field_line(
                fwd,
                fwd_r,
                bwd,
                bwd_r,
                seed,
                direction,
                field_name,
                data.normalization,
                _meta(max(fwd_err, bwd_err)),
            )
        case _ as unreachable:
            assert_never(unreachable)

trace_field_lines_adaptive(data, seeds, *, atol=1e-06, rtol=0.001, step_size_init=0.5, min_step=1e-08, max_step=2.0, max_steps=10000, direction='both', field_components=('B_1', 'B_2', 'B_3'), null_threshold=1e-12, terminate=None, loop_tol='auto', loop_min_arclen=None, interpolator=None)

Trace N field lines adaptively in parallel via the batched DP kernel.

Functionally equivalent to calling trace_field_line_adaptive N times in a Python loop, but each Butcher-stage RHS evaluation is amortized across all seeds in one VectorFieldInterpolator dispatch — ~10× faster for moderate N, larger speedups for N in the thousands. Memory cost is N * (max_steps + 1) * 24 bytes per direction; pick max_steps accordingly for large seed arrays.

Per-seed termination state (live mask, step count, reason, max local error) is tracked in vectorized NumPy; dead seeds report valid=False to the kernel on subsequent steps and don't contaminate the surviving seeds (Butcher contraction is over the stage axis, not the seed axis).

Parameters:

Name	Type	Description	Default
data	FieldDataset	Gridded vector field data.	required
seeds	FloatArray	Starting points, shape (N, 3).	required
atol	float	Adaptive integration controls; semantics identical to trace_field_line_adaptive. Applied per-seed.	1e-06
rtol	float	Adaptive integration controls; semantics identical to trace_field_line_adaptive. Applied per-seed.	1e-06
step_size_init	float	Adaptive integration controls; semantics identical to trace_field_line_adaptive. Applied per-seed.	1e-06
min_step	float	Adaptive integration controls; semantics identical to trace_field_line_adaptive. Applied per-seed.	1e-06
max_step	float	Adaptive integration controls; semantics identical to trace_field_line_adaptive. Applied per-seed.	1e-06
max_steps	float	Adaptive integration controls; semantics identical to trace_field_line_adaptive. Applied per-seed.	1e-06
direction	str	"forward", "backward", or "both".	'both'
field_components	tuple[str, str, str]	Names of the three vector field components.	('B_1', 'B_2', 'B_3')
null_threshold	float	Field magnitude below which a point is treated as a null.	1e-12
terminate	Callable[[FloatArray], bool] | None	Optional per-seed callback; stops a seed when it returns True on the newly accepted point.	None
loop_tol	float, None, or ``"auto"``	Proximity threshold for closed-loop detection, in code units. Applied per-seed: each trace terminates with TerminationReason.CLOSED_LOOP as soon as it re-enters a loop_tol-radius ball around one of its own previously visited points separated by more than loop_min_arclen of arc length. Default "auto" derives 0.5 * min(data.grid.spacing); pass None to disable or a float to override. See trace_field_line_adaptive for the underlying rationale (mirror-mode magnetic holes, O-type islands).	'auto'
loop_min_arclen	float or None	Minimum arc-length separation before a proximity hit counts. Defaults to 10 * step_size_init when loop_tol is set and this is left None.	None
interpolator	VectorFieldInterpolator | None	Pre-built interpolator. Built internally if None.	None

Returns:

Type	Description
list[FieldLine]	One FieldLine per input seed, in seed order.

Raises:

Type	Description
ValueError	If seeds has the wrong shape, if any seed is outside the interpolation domain, or if any seed is at a field null. The upfront validation matches the single-seed contract — invalid seeds fail loudly rather than producing a 1-point trace.

Examples:

>>> import numpy as np
>>> from pypic.dataset import FieldDataset
>>> from pypic.grid import GridInfo
>>> from pypic.units import Normalization
>>> grid = GridInfo(dimensions=(8, 8, 8), spacing=(1.0, 1.0, 1.0))
>>> data = FieldDataset.from_arrays(
...     {
...         "B_1": np.ones((8, 8, 8)),
...         "B_2": np.zeros((8, 8, 8)),
...         "B_3": np.zeros((8, 8, 8)),
...     },
...     grid,
...     Normalization.identity(),
... )
>>> seeds = np.array([[2.0, 2.0, 2.0], [4.0, 4.0, 4.0]])
>>> lines = trace_field_lines_adaptive(
...     data, seeds, max_steps=4, direction="forward"
... )
>>> len(lines)
2
>>> all(fl.metadata["method"] == "rk45_dopri" for fl in lines)
True

Source code in src/pypic/traces/_tracing.py

def trace_field_lines_adaptive(
    data: FieldDataset,
    seeds: FloatArray,
    *,
    atol: float = 1e-6,
    rtol: float = 1e-3,
    step_size_init: float = 0.5,
    min_step: float = 1e-8,
    max_step: float = 2.0,
    max_steps: int = 10_000,
    direction: TraceDirection = "both",
    field_components: tuple[str, str, str] = ("B_1", "B_2", "B_3"),
    null_threshold: float = 1e-12,
    terminate: Callable[[FloatArray], bool] | None = None,
    loop_tol: float | None | Literal["auto"] = "auto",
    loop_min_arclen: float | None = None,
    interpolator: VectorFieldInterpolator | None = None,
) -> list[FieldLine]:
    r"""Trace ``N`` field lines adaptively in parallel via the batched DP kernel.

    Functionally equivalent to calling `trace_field_line_adaptive`
    ``N`` times in a Python loop, but each Butcher-stage RHS evaluation
    is amortized across all seeds in one
    `VectorFieldInterpolator` dispatch — ~10× faster for
    moderate ``N``, larger speedups for ``N`` in the thousands. Memory
    cost is ``N * (max_steps + 1) * 24`` bytes per direction; pick
    ``max_steps`` accordingly for large seed arrays.

    Per-seed termination state (live mask, step count, reason, max
    local error) is tracked in vectorized NumPy; dead seeds report
    ``valid=False`` to the kernel on subsequent steps and don't
    contaminate the surviving seeds (Butcher contraction is over the
    stage axis, not the seed axis).

    Parameters
    ----------
    data : FieldDataset
        Gridded vector field data.
    seeds : FloatArray
        Starting points, shape ``(N, 3)``.
    atol, rtol, step_size_init, min_step, max_step, max_steps : float
        Adaptive integration controls; semantics identical to
        `trace_field_line_adaptive`. Applied per-seed.
    direction : str
        ``"forward"``, ``"backward"``, or ``"both"``.
    field_components : tuple[str, str, str]
        Names of the three vector field components.
    null_threshold : float
        Field magnitude below which a point is treated as a null.
    terminate : Callable[[FloatArray], bool] | None
        Optional per-seed callback; stops a seed when it returns
        ``True`` on the newly accepted point.
    loop_tol : float, None, or ``"auto"``
        Proximity threshold for closed-loop detection, in code units.
        Applied per-seed: each trace terminates with
        `TerminationReason.CLOSED_LOOP` as soon as it re-enters
        a ``loop_tol``-radius ball around one of its own previously
        visited points separated by more than ``loop_min_arclen`` of
        arc length. Default ``"auto"`` derives
        ``0.5 * min(data.grid.spacing)``; pass ``None`` to disable or
        a float to override. See `trace_field_line_adaptive` for
        the underlying rationale (mirror-mode magnetic holes, O-type
        islands).
    loop_min_arclen : float or None
        Minimum arc-length separation before a proximity hit counts.
        Defaults to ``10 * step_size_init`` when ``loop_tol`` is set
        and this is left ``None``.
    interpolator : VectorFieldInterpolator | None
        Pre-built interpolator. Built internally if ``None``.

    Returns
    -------
    list[FieldLine]
        One `FieldLine` per input seed, in seed order.

    Raises
    ------
    ValueError
        If ``seeds`` has the wrong shape, if any seed is outside the
        interpolation domain, or if any seed is at a field null. The
        upfront validation matches the single-seed contract — invalid
        seeds fail loudly rather than producing a 1-point trace.

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.dataset import FieldDataset
    >>> from pypic.grid import GridInfo
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(dimensions=(8, 8, 8), spacing=(1.0, 1.0, 1.0))
    >>> data = FieldDataset.from_arrays(
    ...     {
    ...         "B_1": np.ones((8, 8, 8)),
    ...         "B_2": np.zeros((8, 8, 8)),
    ...         "B_3": np.zeros((8, 8, 8)),
    ...     },
    ...     grid,
    ...     Normalization.identity(),
    ... )
    >>> seeds = np.array([[2.0, 2.0, 2.0], [4.0, 4.0, 4.0]])
    >>> lines = trace_field_lines_adaptive(
    ...     data, seeds, max_steps=4, direction="forward"
    ... )
    >>> len(lines)
    2
    >>> all(fl.metadata["method"] == "rk45_dopri" for fl in lines)
    True
    """
    if direction not in _VALID_DIRECTIONS:
        msg = f"direction must be one of {sorted(_VALID_DIRECTIONS)}, got {direction!r}"
        raise ValueError(msg)

    if loop_tol == "auto":
        loop_tol = _auto_loop_tol(data)
    loop_tol, loop_min_arclen = _resolve_loop_kwargs(
        loop_tol, loop_min_arclen, step_size_init
    )

    if interpolator is None:
        interpolator = VectorFieldInterpolator.from_dataset(data, field_components)

    seeds_arr = np.asarray(seeds, dtype=np.float64)
    if seeds_arr.ndim != 2 or seeds_arr.shape[1] != 3:
        msg = f"seeds must have shape (N, 3), got {seeds_arr.shape}"
        raise ValueError(msg)
    n_seeds = seeds_arr.shape[0]

    # Upfront per-seed validation: invalid seeds raise here rather than
    # producing 1-point traces that violate FieldLine's N >= 2 invariant.
    # O(N) Python overhead, negligible vs. the actual tracing.
    for i in range(n_seeds):
        _validate_seed(seeds_arr[i], interpolator, null_threshold)

    field_name = _field_name_from_components(field_components)
    args = (
        atol,
        rtol,
        step_size_init,
        min_step,
        max_step,
        max_steps,
        null_threshold,
        terminate,
        loop_tol,
        loop_min_arclen,
    )

    def _build(
        fwd_pts: FloatArray,
        fwd_reason: TerminationReason,
        bwd_pts: FloatArray,
        bwd_reason: TerminationReason,
        seed_i: int,
        max_err: float,
    ) -> FieldLine:
        meta: dict[str, Any] = {
            "method": "rk45_dopri",
            "atol": atol,
            "rtol": rtol,
            "max_local_error": max_err,
        }
        return _assemble_field_line(
            fwd_pts,
            fwd_reason,
            bwd_pts,
            bwd_reason,
            tuple(seeds_arr[seed_i].tolist()),
            direction,
            field_name,
            data.normalization,
            meta,
        )

    field_lines: list[FieldLine] = []

    match direction:
        case "forward":
            buf, n_steps, reasons, errs = _trace_batch_single_direction_adaptive(
                interpolator, seeds_arr, 1.0, *args
            )
            field_lines.extend(
                _build(
                    buf[i, : int(n_steps[i]) + 1],
                    _REASON_FROM_INT[int(reasons[i])],
                    _EMPTY_POINTS,
                    TerminationReason.MAX_STEPS,
                    i,
                    float(errs[i]),
                )
                for i in range(n_seeds)
            )
        case "backward":
            buf, n_steps, reasons, errs = _trace_batch_single_direction_adaptive(
                interpolator, seeds_arr, -1.0, *args
            )
            field_lines.extend(
                _build(
                    _EMPTY_POINTS,
                    TerminationReason.MAX_STEPS,
                    buf[i, : int(n_steps[i]) + 1],
                    _REASON_FROM_INT[int(reasons[i])],
                    i,
                    float(errs[i]),
                )
                for i in range(n_seeds)
            )
        case "both":
            fwd_buf, fwd_n, fwd_reasons, fwd_errs = (
                _trace_batch_single_direction_adaptive(
                    interpolator, seeds_arr, 1.0, *args
                )
            )
            bwd_buf, bwd_n, bwd_reasons, bwd_errs = (
                _trace_batch_single_direction_adaptive(
                    interpolator, seeds_arr, -1.0, *args
                )
            )
            field_lines.extend(
                _build(
                    fwd_buf[i, : int(fwd_n[i]) + 1],
                    _REASON_FROM_INT[int(fwd_reasons[i])],
                    bwd_buf[i, : int(bwd_n[i]) + 1],
                    _REASON_FROM_INT[int(bwd_reasons[i])],
                    i,
                    max(float(fwd_errs[i]), float(bwd_errs[i])),
                )
                for i in range(n_seeds)
            )
        case _ as unreachable:
            assert_never(unreachable)

    return field_lines