Core Containers

The data structures every other part of pypic is written against. FieldDataset is the central type: readers produce it, derived quantities and selections consume it, and the I/O layer serializes it.

Dependency direction is one-way — grid ← containers ← dataset — so the geometry description never has to know about the data it describes.

Type	Module	Role
FieldDataset	pypic.dataset	Field arrays plus grid, normalization, and species metadata. Wraps xarray.Dataset; computation happens on raw NumPy.
GridInfo	pypic.grid	Dimensions, spacing, extent, and geometry of the co-located grid.
SimulationConfig	pypic.containers	Validated simulation.toml contents in internal form.
TabularData	pypic.containers	Time-series and other column-oriented output (probes, energy histories).
ParticleData	pypic.containers	Per-macroparticle positions, velocities, and weights.
StaggerInfo	pypic.containers	Provenance record of the source code's mesh convention. Readers destagger on load; this documents what they destaggered from.

FieldDataset

dataset

FieldDataset: the xarray-backed container every reader returns.

compute and reductions sit above this module in the dependency order (grid ← containers ← dataset ← everything else), so the methods that delegate to them import inside the method body. Every other import is at module scope.

FieldDataset

Universal container for simulation field data.

Wraps an xr.Dataset with grid metadata, normalization info, species definitions, and geometry-aware field aliases (e.g. "Bx" → "B_1").

The full alias hierarchy — geometry/Cartesian aliases, species-name aliases (n_electrons→n_s0), and the e/i library-convenience forms (Pe↔P_s0) — is documented at docs/aliases.md. The e/i shortcuts are pypic-only ergonomics and are not part of the cross-tool schema contract in docs/schema.md.

Use from_arrays to construct from raw NumPy arrays.

Parameters:

Name	Type	Description	Default
dataset	Dataset	The underlying xarray dataset.	required
grid	GridInfo	Grid metadata.	required
normalization	Normalization	Unit normalization for this data.	required
species	Sequence[SpeciesInfo] | None	Species definitions, if applicable.	None
physics	PhysicsParams | None	Physics parameters (adiabatic index, speed of light, etc.).	None
metadata	Mapping[str, Any] | None	Arbitrary metadata (run name, code version, etc.).	None
aliases	Mapping[str, str] | None	Extra field-name aliases merged with geometry defaults. Aliases whose canonical target is absent from the dataset are silently dropped (they become inactive).	None
frame	str	Name of the reference frame these arrays are expressed in. Defaults to "simulation".	'simulation'
transforms	Mapping[str, FrameTransform] | None	Frame transforms reachable from frame, keyed by target-frame name. Consumed by transform_to and available_frames; chains resolve by breadth-first search.	None

Examples:

>>> import numpy as np
>>> from pypic.units import Normalization
>>> grid = GridInfo(
...     dimensions=(4, 3), spacing=(1.0, 1.0), origin=(0.0, 0.0),
...     geometry=CARTESIAN,
... )
>>> fields = {"B_1": np.ones((4, 3)), "rho_c": np.zeros((4, 3))}
>>> ds = FieldDataset.from_arrays(fields, grid, Normalization.identity())
>>> ds["B_1"].shape
(4, 3)
>>> ds.has_field("Bx")
True

Source code in src/pypic/dataset.py

class FieldDataset:
    r"""Universal container for simulation field data.

    Wraps an ``xr.Dataset`` with grid metadata, normalization info, species
    definitions, and geometry-aware field aliases (e.g. ``"Bx"`` → ``"B_1"``).

    The full alias hierarchy — geometry/Cartesian aliases, species-name
    aliases (``n_electrons``→``n_s0``), and the e/i library-convenience
    forms (``Pe``↔``P_s0``) — is documented at ``docs/aliases.md``. The
    e/i shortcuts are pypic-only ergonomics and are **not** part of the
    cross-tool schema contract in ``docs/schema.md``.

    Use `from_arrays` to construct from raw NumPy arrays.

    Parameters
    ----------
    dataset : xr.Dataset
        The underlying xarray dataset.
    grid : GridInfo
        Grid metadata.
    normalization : Normalization
        Unit normalization for this data.
    species : Sequence[SpeciesInfo] | None
        Species definitions, if applicable.
    physics : PhysicsParams | None
        Physics parameters (adiabatic index, speed of light, etc.).
    metadata : Mapping[str, Any] | None
        Arbitrary metadata (run name, code version, etc.).
    aliases : Mapping[str, str] | None
        Extra field-name aliases merged with geometry defaults.
        Aliases whose canonical target is absent from the dataset
        are silently dropped (they become inactive).
    frame : str
        Name of the reference frame these arrays are expressed in.
        Defaults to ``"simulation"``.
    transforms : Mapping[str, FrameTransform] | None
        Frame transforms reachable from *frame*, keyed by target-frame
        name.  Consumed by `transform_to` and `available_frames`;
        chains resolve by breadth-first search.

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(
    ...     dimensions=(4, 3), spacing=(1.0, 1.0), origin=(0.0, 0.0),
    ...     geometry=CARTESIAN,
    ... )
    >>> fields = {"B_1": np.ones((4, 3)), "rho_c": np.zeros((4, 3))}
    >>> ds = FieldDataset.from_arrays(fields, grid, Normalization.identity())
    >>> ds["B_1"].shape
    (4, 3)
    >>> ds.has_field("Bx")
    True
    """

    def __init__(
        self,
        dataset: xr.Dataset,
        grid: GridInfo,
        normalization: Normalization,
        *,
        species: Sequence[SpeciesInfo] | None = None,
        physics: PhysicsParams | None = None,
        metadata: Mapping[str, Any] | None = None,
        aliases: Mapping[str, str] | None = None,
        frame: str = "simulation",
        transforms: Mapping[str, FrameTransform] | None = None,
    ) -> None:
        _require_real_arrays(dataset)
        self._ds = dataset
        self._grid = grid
        self._normalization = normalization
        self._species = tuple(species) if species is not None else ()
        self._physics = physics if physics is not None else PhysicsParams()
        self._metadata = dict(metadata) if metadata is not None else {}
        self._frame = frame
        self._transforms = dict(transforms) if transforms is not None else {}

        merged = _default_aliases(grid.geometry)
        if aliases:
            merged.update(aliases)
        # Generate species-name aliases for every per-species canonical
        # actually in the dataset (n_electrons→n_s0, P_ions→P_s1, etc.).
        merged.update(
            species_name_aliases(
                tuple(sp.name for sp in self._species),
                [str(name) for name in self._ds.data_vars],
            )
        )
        # Bidirectional alias filter: ``alias→canonical`` when the canonical
        # is stored, and ``canonical→alias`` when the data sits under what is
        # now the alias (``Pe``, since v1.0 made ``P_s0`` canonical).  Without
        # the reverse pass a recipe asking for ``_sN`` would miss.
        data_vars_set = set(self._ds.data_vars)
        filtered: dict[str, str] = {}
        for alias_name, canonical_name in merged.items():
            if canonical_name in data_vars_set:
                filtered[alias_name] = canonical_name
            elif alias_name in data_vars_set:
                filtered.setdefault(canonical_name, alias_name)
        # Second pass — transitive chains: ``P_e → P_s0 → Pe`` collapses
        # to ``P_e → Pe`` when only ``Pe`` is stored.  Iterates to a
        # fixed point in O(|merged|) per round; chain depth is bounded
        # by the alias graph (≤2 in practice).
        for _round in range(3):
            changed = False
            for alias_name, canonical_name in merged.items():
                if alias_name in filtered or alias_name in data_vars_set:
                    continue
                if canonical_name in filtered:
                    filtered[alias_name] = filtered[canonical_name]
                    changed = True
            if not changed:
                break
        self._aliases = filtered

    @classmethod
    def from_arrays(
        cls,
        fields: Mapping[str, FloatArray],
        grid: GridInfo,
        normalization: Normalization | None = None,
        *,
        species: Sequence[SpeciesInfo] | None = None,
        physics: PhysicsParams | None = None,
        metadata: Mapping[str, Any] | None = None,
        aliases: Mapping[str, str] | None = None,
        frame: str = "simulation",
        transforms: Mapping[str, FrameTransform] | None = None,
        coords: Mapping[str, FloatArray] | None = None,
        strict_fields: bool = True,
    ) -> FieldDataset:
        r"""Build a FieldDataset from a dict of NumPy arrays.

        Parameters
        ----------
        fields : Mapping[str, FloatArray]
            Mapping of field names to arrays. Shapes must match
            ``grid.dimensions``.  Names must resolve through the field
            registry (canonical names, registered aliases, or
            species-templated patterns) when ``strict_fields=True``.
        grid : GridInfo
            Grid metadata.
        normalization : Normalization | None
            Unit normalization.  Defaults to
            ``Normalization.undeclared()`` — arrays are left alone, but
            `in_si` and `in_units` raise on dimensional quantities
            rather than return code units labelled SI.  Pass
            ``Normalization.identity()`` to assert the arrays already
            *are* SI.
        species : Sequence[SpeciesInfo] | None
            Species definitions, if applicable.
        physics : PhysicsParams | None
            Physics parameters.
        metadata : Mapping[str, Any] | None
            Arbitrary metadata.
        aliases : Mapping[str, str] | None
            Extra field-name aliases.
        frame : str
            Name of the reference frame these arrays are expressed in.
            Defaults to ``"simulation"``.
        transforms : Mapping[str, FrameTransform] | None
            Frame transforms reachable from *frame*, keyed by
            target-frame name.  Consumed by `transform_to`.
        coords : Mapping[str, FloatArray] | None
            Coordinate arrays keyed by axis name, replacing the uniform
            ones derived from *grid* on those axes.  For non-uniform
            meshes whose ``GridInfo`` spacing is only a mean; the true
            positions then live in the xarray coordinates.
        strict_fields : bool
            When ``True`` (default), every key in *fields* must resolve
            through ``pypic.fields.field_info`` — any unknown name
            raises ``KeyError`` with the full list of unresolved names.
            Set ``False`` to inject scratch fields (test fixtures, ad-hoc
            scalars) that aren't part of the canonical schema; the field
            is then stored without registry metadata.

        Returns
        -------
        FieldDataset

        Raises
        ------
        UnknownFieldError
            When ``strict_fields=True`` and one or more keys in *fields*
            do not resolve through the field registry.  Injection
            points fail loud on unknown names so reader bugs surface at
            construction time instead of later at ``compute()``.  Subclass of
            `KeyError`, so existing ``except KeyError`` callers
            keep working unchanged.
        ValueError
            When *coords* names an axis the grid does not have.

        Examples
        --------
        >>> import numpy as np
        >>> grid = GridInfo(dimensions=(2,), spacing=(1.0,))
        >>> ds = FieldDataset.from_arrays({"B_1": np.array([1.0, 2.0])}, grid)
        >>> ds["B_1"]
        array([1., 2.])
        """
        if normalization is None:
            normalization = Normalization.undeclared()
        dim_names = list(grid.surviving_axis_names)
        coord_arrays = grid.coordinate_arrays()
        axis_coords = dict(zip(dim_names, coord_arrays, strict=True))
        if coords:
            if unknown := sorted(set(coords) - set(dim_names)):
                msg = f"coords name axes the grid lacks: {unknown}; axes: {dim_names}"
                raise ValueError(msg)
            axis_coords.update(coords)

        data_vars: dict[str, xr.DataArray] = {}
        unresolved: list[str] = []
        for var_name, arr in fields.items():
            da = xr.DataArray(data=arr, dims=dim_names)
            try:
                info = _field_info(var_name, axis_names=grid.geometry.axis_names)
                da.attrs["long_name"] = info.long_name
                da.attrs["units"] = "normalized"
                da.attrs["quantity_type"] = info.quantity_type
                da.attrs["si_unit"] = info.si_unit
                if info.latex:
                    da.attrs["latex"] = info.latex
                ud = info.unit_dimension or quantity_dimension(info.quantity_type)
                da.attrs["unit_dimension"] = list(ud)
            except KeyError:
                if strict_fields:
                    unresolved.append(var_name)
            data_vars[var_name] = da
        if unresolved:
            msg = (
                f"Field names not in the registry: {unresolved!r}. "
                "Pass strict_fields=False to inject scratch fields, or "
                "use FieldDataset.with_field(name, data, quantity_type=...) "
                "to register an ad-hoc quantity_type."
            )
            raise UnknownFieldError(msg)
        dataset = xr.Dataset(data_vars, coords=axis_coords)
        return cls(
            dataset,
            grid,
            normalization,
            species=species,
            physics=physics,
            metadata=metadata,
            aliases=aliases,
            frame=frame,
            transforms=transforms,
        )

    @property
    def grid(self) -> GridInfo:
        """Grid metadata."""
        return self._grid

    @property
    def normalization(self) -> Normalization:
        """Unit normalization."""
        return self._normalization

    @property
    def species(self) -> tuple[SpeciesInfo, ...]:
        """Species definitions."""
        return self._species

    @property
    def physics(self) -> PhysicsParams:
        """Physics parameters."""
        return self._physics

    @property
    def metadata(self) -> MappingProxyType[str, Any]:
        """Arbitrary metadata (read-only view)."""
        return MappingProxyType(self._metadata)

    @property
    def time(self) -> float | None:
        """Snapshot time in code units, or ``None`` when unknown.

        ``metadata["time"]`` when the reader recorded one, otherwise
        ``metadata["step"] * grid.dt`` when both are known.

        Examples
        --------
        >>> import numpy as np
        >>> grid = GridInfo(dimensions=(2,), spacing=(1.0,), dt=0.5)
        >>> FieldDataset.from_arrays(
        ...     {"B_1": np.zeros(2)}, grid, metadata={"step": 4}
        ... ).time
        2.0
        >>> FieldDataset.from_arrays({"B_1": np.zeros(2)}, grid).time is None
        True
        """
        time = self._metadata.get("time")
        if time is not None:
            return float(time)
        step = self._metadata.get("step")
        if step is None or self._grid.dt is None:
            return None
        return float(step) * self._grid.dt

    @property
    def aliases(self) -> MappingProxyType[str, str]:
        """Active field-name aliases (read-only view)."""
        return MappingProxyType(self._aliases)

    @property
    def frame(self) -> str:
        """Current reference frame label."""
        return self._frame

    @property
    def transforms(self) -> MappingProxyType[str, FrameTransform]:
        """Registered frame transforms (read-only view)."""
        return MappingProxyType(self._transforms)

    @property
    def available_frames(self) -> list[str]:
        """Frame names reachable via registered transforms."""
        frames: set[str] = {self._frame}
        for t in self._transforms.values():
            frames.add(t.source_frame)
            frames.add(t.target_frame)
        return sorted(frames)

    @property
    def xr(self) -> xr.Dataset:
        """Raw xarray Dataset."""
        return self._ds

    def transform_to(self, target: str | FrameTransform) -> FieldDataset:
        r"""Transform this dataset to a different coordinate reference frame.

        Applies an affine transformation: translates the grid origin,
        scales coordinates, and rotates vector field components and
        pressure tensors. Scalar values pass through unchanged. When the
        rotation includes an axis swap or reflection, field arrays are
        copied to a contiguous buffer in the new axis order; pure
        translations and identity rotations only touch metadata.

        Parameters
        ----------
        target : str or FrameTransform
            Target frame name (looked up in the transform registry) or
            a `FrameTransform` instance applied directly.

        Returns
        -------
        FieldDataset
            New dataset in the target frame.

        Raises
        ------
        ValueError
            If no transform path exists from the current frame.
        KeyError
            If *target* is a string and no transforms are registered.

        Examples
        --------
        A -90° rotation about $z$ sends $\hat{x}$ to $-\hat{y}$, and
        swaps the two grid axes with it:

        >>> import numpy as np
        >>> from pypic.coordinates import FrameTransform
        >>> from pypic.units import Normalization
        >>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0))
        >>> gsm = FrameTransform(
        ...     "simulation", "GSM",
        ...     rotation=((0.0, 1.0, 0.0), (-1.0, 0.0, 0.0), (0.0, 0.0, 1.0)),
        ... )
        >>> ds = FieldDataset.from_arrays(
        ...     {"B_1": np.ones((4, 3, 2)),
        ...      "B_2": np.zeros((4, 3, 2)),
        ...      "B_3": np.zeros((4, 3, 2))},
        ...     grid, Normalization.identity(), transforms={"GSM": gsm},
        ... )
        >>> rotated = ds.transform_to("GSM")
        >>> rotated.frame, rotated.grid.dimensions
        ('GSM', (3, 4, 2))
        >>> float(rotated["B_2"][0, 0, 0])
        -1.0
        """
        if isinstance(target, FrameTransform):
            transform = target
            target_frame = transform.target_frame
        else:
            target_frame = target
            if target_frame == self._frame:
                return self
            if not self._transforms:
                msg = (
                    f"No transforms registered; cannot transform "
                    f"from {self._frame!r} to {target_frame!r}"
                )
                raise KeyError(msg)
            transform = resolve_transform(self._frame, target_frame, self._transforms)
        rotation = transform.rotation_matrix

        # Determine which axes survive and derive the axis permutation
        # and sign map from the rotation matrix up front. All field
        # arrays then flow through a single reorient → construct pass.
        if self._grid.surviving_axes is not None:
            indices = list(self._grid.surviving_axes)
        else:
            indices = list(range(len(self._grid.dimensions)))
        ndim = len(indices)

        # The grid transformation (transpose + flip) only works for
        # signed permutation matrices (axis swaps and reflections).
        # General rotations would require interpolation onto a new grid.
        rotation_sub = rotation[np.ix_(indices, indices)]
        abs_rotation = np.abs(rotation_sub)
        if not (
            np.allclose(np.sum(abs_rotation, axis=1), 1.0, atol=1e-6)
            and np.allclose(np.sum(abs_rotation, axis=0), 1.0, atol=1e-6)
        ):
            # GeometryUnsupportedError subclasses NotImplementedError, so
            # existing handlers still catch it, and the server maps it to
            # 400 rather than a 500 "internal" frame.
            raise GeometryUnsupportedError(
                "transform_to() only supports axis-swap/reflection "
                "rotation matrices (signed permutations). General "
                "rotations require grid interpolation (not yet implemented)."
            )

        # Axis permutation (target_i → source local index) and sign.
        axis_permutation: list[int] = []
        axis_signs: list[float] = []
        for target_i in indices:
            row = rotation[target_i, :]
            source_orig = int(np.argmax(np.abs(row[indices])))
            axis_permutation.append(source_orig)
            axis_signs.append(float(np.sign(row[indices[source_orig]])))

        needs_transpose = axis_permutation != list(range(ndim))
        flip_axes = [i for i, s in enumerate(axis_signs) if s < 0]

        def _reorient(arr: FloatArray) -> FloatArray:
            if needs_transpose:
                arr = np.transpose(arr, axis_permutation)
            for ax in flip_axes:
                arr = np.flip(arr, axis=ax)
            if needs_transpose or flip_axes:
                arr = np.ascontiguousarray(arr)
            return arr

        # Accumulate raw numpy arrays — no intermediate DataArray wrapping.
        new_arrays: dict[str, FloatArray] = {}
        rotated: set[str] = set()

        for n1, n2, n3 in find_vector_triplets(self.field_names()):
            r1, r2, r3 = rotate_vector_components(
                self[n1], self[n2], self[n3], rotation
            )
            for name, arr in ((n1, r1), (n2, r2), (n3, r3)):
                new_arrays[name] = _reorient(arr)
                rotated.add(name)

        for p11, p22, p33, p12, p13, p23 in find_pressure_tensor_groups(
            self.field_names()
        ):
            rp = rotate_pressure_tensor(
                self[p11],
                self[p22],
                self[p33],
                self[p12],
                self[p13],
                self[p23],
                rotation,
            )
            for name, arr in zip((p11, p22, p33, p12, p13, p23), rp, strict=True):
                new_arrays[name] = _reorient(arr)
                rotated.add(name)

        for raw_name in self._ds.data_vars:
            name = str(raw_name)
            if name not in rotated:
                new_arrays[name] = _reorient(self._ds[name].values)

        # Compute new origin, spacing, dimensions from permuted source
        transform_origin = np.array(transform.origin, dtype=np.float64)
        dx_scale = transform.scale
        old_coords = self._grid.coordinate_arrays()
        new_origin_list: list[float] = []
        new_spacing_list: list[float] = []
        new_dims_list: list[int] = []
        for src_i, sign in zip(axis_permutation, axis_signs, strict=True):
            src_coords = old_coords[src_i]
            coord_first = (
                dx_scale * sign * (src_coords[0] - transform_origin[indices[src_i]])
            )
            coord_last = (
                dx_scale * sign * (src_coords[-1] - transform_origin[indices[src_i]])
            )
            dx = dx_scale * self._grid.spacing[src_i]
            new_origin_list.append(float(min(coord_first, coord_last)) - 0.5 * dx)
            new_spacing_list.append(dx)
            new_dims_list.append(self._grid.dimensions[src_i])

        new_surviving = None
        if self._grid.surviving_axes is not None:
            new_surviving = tuple(
                self._grid.surviving_axes[axis_permutation[i]] for i in range(ndim)
            )

        new_geometry = self._grid.geometry
        if transform.target_axis_names is not None:
            new_geometry = copy.replace(
                new_geometry, axis_names=transform.target_axis_names
            )

        new_grid = copy.replace(
            self._grid,
            dimensions=tuple(new_dims_list),
            origin=tuple(new_origin_list),
            spacing=tuple(new_spacing_list),
            geometry=new_geometry,
            surviving_axes=new_surviving,
        )

        # Single Dataset construction: each field is wrapped exactly once.
        new_dim_names = list(new_grid.surviving_axis_names)
        coord_arrays = new_grid.coordinate_arrays()
        coords = {new_dim_names[i]: coord_arrays[i] for i in range(len(new_dim_names))}
        new_ds = xr.Dataset(
            data_vars={
                name: (new_dim_names, arr, dict(self._ds[name].attrs))
                for name, arr in new_arrays.items()
            },
            coords=coords,
        )

        return FieldDataset(
            new_ds,
            new_grid,
            self._normalization,
            species=self._species,
            physics=self._physics,
            metadata=self._metadata,
            frame=target_frame,
            transforms=self._transforms,
        )

    def _wrap_sliced(self, new_ds: Dataset) -> FieldDataset:
        """Wrap a sliced xr.Dataset in a new FieldDataset, preserving metadata."""
        new_grid = _build_grid_from_dataset(self._grid, new_ds)
        return FieldDataset(
            new_ds,
            new_grid,
            self._normalization,
            species=self._species,
            physics=self._physics,
            metadata=self._metadata,
            aliases={k: v for k, v in self._aliases.items() if v in new_ds.data_vars},
            frame=self._frame,
            transforms=self._transforms,
        )

    def resolve_key(self, key: str) -> str:
        """Resolve a field key through aliases to the canonical name.

        Returns the canonical (data-vars) name for *key*, which may itself
        be a canonical name or an alias from this dataset's alias table.
        Raises `KeyError` (with close-match suggestions) if *key*
        matches neither.

        Parameters
        ----------
        key : str
            Canonical field name or alias.

        Returns
        -------
        str
            Canonical name in ``self._ds.data_vars``.

        Examples
        --------
        >>> import numpy as np
        >>> from pypic.units import Normalization
        >>> grid = GridInfo(
        ...     dimensions=(2,), spacing=(1.0,), origin=(0.0,),
        ...     geometry=CARTESIAN,
        ... )
        >>> ds = FieldDataset.from_arrays(
        ...     {"B_1": np.ones(2)}, grid, Normalization.identity()
        ... )
        >>> ds.resolve_key("Bx")  # Cartesian alias
        'B_1'
        >>> ds.resolve_key("B_1")  # canonical
        'B_1'
        """
        if key in self._ds.data_vars:
            return key
        canonical = self._aliases.get(key)
        if canonical is not None:
            return canonical
        available = sorted(self._ds.data_vars, key=str)
        alias_keys = sorted(self._aliases)
        msg = f"Field {key!r} not found. Available: {available}. Aliases: {alias_keys}."
        import difflib

        candidates = [str(v) for v in self._ds.data_vars] + list(self._aliases)
        suggestions = difflib.get_close_matches(key, candidates, n=3, cutoff=0.5)
        if suggestions:
            msg += f" Did you mean: {suggestions}?"
        raise UnknownFieldError(msg)

    def __getitem__(self, key: str) -> FloatArray:
        """Return field data as a NumPy array (zero-copy when possible).

        Parameters
        ----------
        key : str
            Canonical field name or alias.

        Returns
        -------
        NDArray
        """
        resolved = self.resolve_key(key)
        return self._ds[resolved].values

    def has_field(self, key: str) -> bool:
        """Check whether a field exists (canonical or alias).

        Parameters
        ----------
        key : str
            Field name to check.

        Returns
        -------
        bool

        Examples
        --------
        >>> import numpy as np
        >>> from pypic.units import Normalization
        >>> grid = GridInfo(
        ...     dimensions=(2,), spacing=(1.0,), origin=(0.0,),
        ...     geometry=CARTESIAN,
        ... )
        >>> ds = FieldDataset.from_arrays(
        ...     {"B_1": np.array([1.0, 2.0])}, grid, Normalization.identity(),
        ... )
        >>> ds.has_field("B_1"), ds.has_field("Bx"), ds.has_field("rho")
        (True, True, False)
        """
        return key in self._ds.data_vars or key in self._aliases

    def field_names(self) -> list[str]:
        """Return canonical field names (no aliases).

        Returns
        -------
        list[str]

        Examples
        --------
        >>> import numpy as np
        >>> from pypic.units import Normalization
        >>> grid = GridInfo(
        ...     dimensions=(2,), spacing=(1.0,), origin=(0.0,),
        ...     geometry=CARTESIAN,
        ... )
        >>> ds = FieldDataset.from_arrays(
        ...     {"B_1": np.array([1.0, 2.0]), "rho_c": np.array([0.5, 0.5])},
        ...     grid, Normalization.identity(),
        ... )
        >>> sorted(ds.field_names())
        ['B_1', 'rho_c']
        """
        return list(self._ds.data_vars)  # type: ignore[arg-type]  # xarray types Hashable, always str

    def select_fields(self, names: Iterable[str]) -> FieldDataset:
        """Return a new FieldDataset containing only the specified fields.

        Parameters
        ----------
        names : Iterable[str]
            Field names to keep (canonical or alias).

        Returns
        -------
        FieldDataset

        Raises
        ------
        KeyError
            If any name is not found.

        Examples
        --------
        >>> import numpy as np
        >>> from pypic.units import Normalization
        >>> grid = GridInfo(
        ...     dimensions=(2,), spacing=(1.0,), origin=(0.0,),
        ...     geometry=CARTESIAN,
        ... )
        >>> ds = FieldDataset.from_arrays(
        ...     {"B_1": np.array([1.0, 2.0]), "B_2": np.array([3.0, 4.0]),
        ...      "rho_c": np.array([0.5, 0.5])},
        ...     grid, Normalization.identity(),
        ... )
        >>> sub = ds.select_fields(["rho_c", "B_1"])
        >>> sub.field_names()
        ['rho_c', 'B_1']
        """
        # Dict-as-ordered-set: preserves caller insertion order while
        # deduplicating on the resolved canonical name. Alphabetical
        # sorting hides both request order and dataset order from the
        # caller, so we keep the order the caller asked for.
        resolved: dict[str, None] = {}
        for name in names:
            resolved[self.resolve_key(name)] = None

        new_ds = self._ds[list(resolved)]
        return FieldDataset(
            new_ds,
            self._grid,
            self._normalization,
            species=self._species,
            physics=self._physics,
            metadata=self._metadata,
            aliases={k: v for k, v in self._aliases.items() if v in resolved},
            frame=self._frame,
            transforms=self._transforms,
        )

    def sel(
        self,
        indexers: dict[str, Any] | None = None,
        *,
        method: str | None = None,
        **kwargs: Any,  # noqa: ANN401 — xarray passthrough
    ) -> FieldDataset:
        """Label-based selection, returning a new FieldDataset.

        Accepts ``indexers`` as a dict (useful for Unicode dim names like
        ``"θ"`` that can't be keyword arguments) and/or ``**kwargs``.

        Parameters
        ----------
        indexers : dict[str, Any] | None
            Dimension-name → label mapping.
        method : str | None
            Passed to ``xr.Dataset.sel`` (e.g. ``"nearest"``).
        **kwargs
            Additional dimension selections.

        Returns
        -------
        FieldDataset

        Examples
        --------
        Cell centres sit at half-integer positions, so a label that
        falls between them needs ``method="nearest"``:

        >>> import numpy as np
        >>> from pypic.units import Normalization
        >>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0))
        >>> ds = FieldDataset.from_arrays(
        ...     {"B_1": np.arange(24.0).reshape(4, 3, 2)},
        ...     grid, Normalization.identity(),
        ... )
        >>> midplane = ds.sel({"z": 1.4}, method="nearest")
        >>> midplane.grid.dimensions, midplane.grid.surviving_axis_names
        ((4, 3), ('x', 'y'))
        """
        merged = dict(indexers) if indexers else {}
        merged.update(kwargs)
        return self._wrap_sliced(self._ds.sel(merged, method=method))

    def compute(self, name: str) -> FloatArray:
        """Compute a derived quantity by name. Returns code units.

        Parameters
        ----------
        name : str
            Derived quantity name (e.g. ``"|B|"``, ``"beta"``, ``"v_A"``).
            See ``pypic.compute.available_quantities()`` for the full list.

        Returns
        -------
        FloatArray
            Computed array in code units.

        Examples
        --------
        >>> import numpy as np
        >>> from pypic.units import Normalization
        >>> grid = GridInfo(
        ...     dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0),
        ...     geometry=CARTESIAN,
        ... )
        >>> ds = FieldDataset.from_arrays(
        ...     {"B_1": np.full((4,3,2), 3.0),
        ...      "B_2": np.full((4,3,2), 4.0),
        ...      "B_3": np.zeros((4,3,2))},
        ...     grid, Normalization.identity(),
        ... )
        >>> ds.compute("|B|")[0, 0, 0]
        np.float64(5.0)
        """
        from pypic.compute import compute_field  # layered above dataset

        return compute_field(name, self)

    def with_field(
        self,
        name: str,
        data: FloatArray,
        quantity_type: QuantityType | str | None = None,
        *,
        long_name: str = "",
        latex: str = "",
    ) -> FieldDataset:
        """Return a new FieldDataset with an additional custom field.

        The field's ``quantity_type`` is stored in xarray attrs, so
        ``in_si()``, ``field_info()``, and unit conversion work without
        global ``register_field()`` calls.

        When *quantity_type* is ``None``, metadata is looked up from the
        field registry automatically. For unregistered fields, provide
        *quantity_type* explicitly.

        Parameters
        ----------
        name : str
            Field name.
        data : FloatArray
            Array matching the grid dimensions.
        quantity_type : QuantityType | str | None
            Physical quantity type (e.g. ``QuantityType.VELOCITY``).
            When ``None``, auto-filled from the field registry.
        long_name : str
            Human-readable label for plot titles.
        latex : str
            LaTeX symbol for plot labels.

        Returns
        -------
        FieldDataset
            New dataset with the field added.

        Raises
        ------
        ValueError
            If *quantity_type* is not recognized, or is ``None`` and
            the field name is not in the registry.
        """
        expected = tuple(self._grid.dimensions)
        if data.shape != expected:
            msg = f"Array shape {data.shape} doesn't match grid dimensions {expected}"
            raise ValueError(msg)

        info_ud: tuple[int, int, int, int, int, int, int] | None = None
        if quantity_type is None:
            try:
                info = _field_info(name, axis_names=self._grid.geometry.axis_names)
            except KeyError:
                msg = f"Unknown field {name!r}; provide quantity_type explicitly"
                raise ValueError(msg) from None
            qt = info.quantity_type
            if not long_name:
                long_name = info.long_name
            if not latex:
                latex = info.latex
            info_ud = info.unit_dimension
        else:
            qt = str(quantity_type)

        if qt not in _QUANTITY_UNITS:
            valid = sorted(_QUANTITY_UNITS)
            msg = f"Unknown quantity_type {qt!r}. Valid: {valid}"
            raise ValueError(msg)

        si_unit = _QUANTITY_UNITS[qt]
        dim_names = list(self._grid.surviving_axis_names)
        da = xr.DataArray(data=data, dims=dim_names)
        da.attrs["quantity_type"] = qt
        da.attrs["si_unit"] = si_unit
        da.attrs["units"] = "normalized"
        da.attrs["long_name"] = long_name
        da.attrs["latex"] = latex
        da.attrs["unit_dimension"] = list(info_ud or quantity_dimension(qt))

        new_ds = self._ds.assign({name: da})
        return FieldDataset(
            new_ds,
            self._grid,
            self._normalization,
            species=self._species,
            physics=self._physics,
            metadata=self._metadata,
            aliases=dict(self._aliases),
            frame=self._frame,
            transforms=self._transforms,
        )

    def with_derived(self, *names: str) -> FieldDataset:
        """Return a new dataset with derived fields computed and stored.

        For each name, computes the field (if not already present),
        looks up metadata from the field registry, and attaches it
        with full metadata. Fields already in the dataset are skipped.

        When a vector-component recipe is encountered (e.g. ``"S_1"``
        from Poynting flux), all sibling components (``"S_2"``, ``"S_3"``)
        are computed from a single function call and stored together.

        Parameters
        ----------
        *names : str
            Derived quantity names (e.g. ``"|B|"``, ``"beta"``).

        Returns
        -------
        FieldDataset
            New dataset with the computed fields attached.

        Examples
        --------
        >>> import numpy as np
        >>> from pypic.units import Normalization
        >>> grid = GridInfo(
        ...     dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0),
        ...     geometry=CARTESIAN,
        ... )
        >>> ds = FieldDataset.from_arrays(
        ...     {"B_1": np.full((4,3,2), 3.0),
        ...      "B_2": np.full((4,3,2), 4.0),
        ...      "B_3": np.zeros((4,3,2))},
        ...     grid, Normalization.identity(),
        ... )
        >>> ds = ds.with_derived("|B|")
        >>> ds.has_field("|B|")
        True
        >>> ds["|B|"][0, 0, 0]
        np.float64(5.0)
        """
        from pypic.compute import compute_with_siblings  # layered above dataset

        result = self
        for name in names:
            if result.has_field(name):
                continue
            for field, data in compute_with_siblings(name, result).items():
                if not result.has_field(field):
                    result = result.with_field(field, data)
        return result

    def field_info(self, name: str) -> FieldInfo:
        """Return metadata for a field or derived quantity.

        Checks xarray DataArray attrs first (set by ``with_field()``
        or ``from_arrays()``), then falls back to the global registry.

        Parameters
        ----------
        name : str
            Field or derived quantity name (canonical or alias).

        Returns
        -------
        FieldInfo
        """
        if self.has_field(name):
            resolved = self.resolve_key(name)
            attrs = self._ds[resolved].attrs
            qt = attrs.get("quantity_type")
            if qt is not None:
                return FieldInfo(
                    quantity_type=qt,
                    long_name=attrs.get("long_name", ""),
                    si_unit=attrs.get("si_unit", ""),
                    latex=attrs.get("latex", ""),
                )
        return _field_info(name, axis_names=self._grid.geometry.axis_names)

    def in_si(self, name: str) -> FloatArray:
        r"""Return a field or derived quantity in SI units.

        Checks xarray DataArray attrs first (set by ``with_field()``
        or ``from_arrays()``), then falls back to the global registry.

        For fields that have been reduced via
        [`pypic.reductions.reduce`][pypic.reductions.reduce] with
        ``reduction="integrate"``,
        the ``attrs["reduction"]["length_axes"]`` provenance stamp
        records how many length-dimension factors the integration
        added; ``in_si`` multiplies the registry SI factor by
        ``normalization.length_ref ** length_axes`` so column
        densities, line-of-sight integrals, etc. come out in the
        right SI units (e.g. column density m\ :sup:`-2` instead of
        m\ :sup:`-3`).  Weighted ``integrate`` does not stamp
        ``length_axes`` because the length factor cancels in
        ``∫ f w dx / ∫ w dx``.

        Parameters
        ----------
        name : str
            Field or derived quantity name.

        Returns
        -------
        FloatArray
            Values in SI units.

        Raises
        ------
        ValueError
            If no SI conversion is registered for *name* — a field
            carrying no ``quantity_type`` attr whose name the global
            registry also cannot resolve.

        Examples
        --------
        >>> import numpy as np
        >>> from pypic.units import Normalization
        >>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0))
        >>> norm = Normalization.mhd_standard(6.371e6, 1.67e-17, 5.0e-9)
        >>> ds = FieldDataset.from_arrays(
        ...     {"B_1": np.full((4, 3, 2), 2.0)}, grid, norm,
        ... )
        >>> float(ds.in_si("B_1")[0, 0, 0])  # 2 code units at B_ref = 5 nT
        1e-08
        """
        from pypic.compute import compute_field, field_si_factor  # above dataset

        length_axes = 0
        factor: float | None = None
        if self.has_field(name):
            resolved = self.resolve_key(name)
            data = self._ds[resolved].values
            reduction_attr = self._ds[resolved].attrs.get("reduction") or {}
            length_axes = int(reduction_attr.get("length_axes", 0))
            qt = self._ds[resolved].attrs.get("quantity_type")
            if qt is not None:
                factor = self._normalization.si_factor(qt)
        else:
            data = compute_field(name, self)
        if factor is None:
            # Either a derived quantity, or a stored field with no
            # ``quantity_type`` attr — what
            # ``from_arrays(..., strict_fields=False)`` and
            # ``open_virtual`` produce. Both resolve through the registry.
            factor = field_si_factor(name, self._normalization)
        if length_axes:
            factor *= self._normalization.length_ref**length_axes
        return data if factor == 1.0 else data * factor

    def in_units(self, name: str, unit_str: str) -> FloatArray:
        """Return a field or derived quantity in display units.

        Parameters
        ----------
        name : str
            Field or derived quantity name.
        unit_str : str
            Target unit (e.g. ``"nT"``, ``"km/s"``, ``"cm^-3"``).

        Returns
        -------
        FloatArray
            Values in the requested units.

        Notes
        -----
        Temperature is stored in **energy units** (J), so use
        ``in_units(name, "eV")`` (or ``"keV"``, ``"MeV"``) for the
        plasma working unit, or ``in_units(name, "K")`` for the
        Boltzmann-factor-converted form. Magnetic field defaults to
        T in SI; ``in_units(name, "nT")`` is the space-physics
        idiom. See ``_DISPLAY_UNITS`` in ``pypic.units`` for the
        full vocabulary.

        Examples
        --------
        >>> import numpy as np
        >>> from pypic.units import Normalization
        >>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0))
        >>> norm = Normalization.mhd_standard(6.371e6, 1.67e-17, 5.0e-9)
        >>> ds = FieldDataset.from_arrays(
        ...     {"B_1": np.full((4, 3, 2), 2.0)}, grid, norm,
        ... )
        >>> float(ds.in_units("B_1", "nT")[0, 0, 0])
        10.0
        """
        from pypic.compute import display_unit_factor  # layered above dataset

        return self.in_si(name) / display_unit_factor(unit_str)

    def isel(
        self,
        indexers: dict[str, Any] | None = None,
        **kwargs: Any,  # noqa: ANN401 — xarray passthrough
    ) -> FieldDataset:
        """Integer-index selection, returning a new FieldDataset.

        Accepts ``indexers`` as a dict and/or ``**kwargs``.

        Parameters
        ----------
        indexers : dict[str, Any] | None
            Dimension-name → integer index or slice mapping.
        **kwargs
            Additional dimension selections.

        Returns
        -------
        FieldDataset

        Examples
        --------
        >>> import numpy as np
        >>> from pypic.units import Normalization
        >>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0))
        >>> ds = FieldDataset.from_arrays(
        ...     {"B_1": np.arange(24.0).reshape(4, 3, 2)},
        ...     grid, Normalization.identity(),
        ... )
        >>> face = ds.isel(x=0)
        >>> face.grid.dimensions, face.grid.surviving_axis_names
        ((3, 2), ('y', 'z'))
        """
        merged = dict(indexers) if indexers else {}
        merged.update(kwargs)
        return self._wrap_sliced(self._ds.isel(merged))

    def where(self, cond: BoolArray, other: float = np.nan) -> FieldDataset:
        r"""Mask fields where *cond* is ``False``.

        Returns a new `FieldDataset` with the same grid shape.
        Points where *cond* is ``False`` are set to *other* (default
        ``NaN``).  Useful for spatial masks (spherical cutouts, boundary
        regions) without reducing dimensions.

        Parameters
        ----------
        cond : BoolArray
            Boolean array with shape matching the grid dimensions.
            ``True`` keeps the value, ``False`` replaces with *other*.
        other : float
            Fill value for masked points (default ``NaN``).

        Returns
        -------
        FieldDataset

        Examples
        --------
        >>> import numpy as np
        >>> from pypic.units import Normalization
        >>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0))
        >>> ds = FieldDataset.from_arrays(
        ...     {"B_1": np.arange(24.0).reshape(4, 3, 2)},
        ...     grid, Normalization.identity(),
        ... )
        >>> x, _, _ = np.meshgrid(*grid.coordinate_arrays(), indexing="ij")
        >>> inner = ds.where(x < 2.0)
        >>> inner["B_1"].shape, int(np.isnan(inner["B_1"]).sum())
        ((4, 3, 2), 12)
        """
        axis_names = list(self._grid.surviving_axis_names)
        mask_da = xr.DataArray(cond, dims=axis_names)
        return self._wrap_sliced(self._ds.where(mask_da, other=other))

    def reduce(
        self,
        axis: str | tuple[str, ...],
        **kwargs: Any,  # noqa: ANN401
    ) -> FieldDataset:
        r"""Reduce this dataset along one or more axes.

        See [`pypic.reduce`][pypic.reduce]. Convenience method equivalent to
        ``pypic.reduce(self, axis, **kwargs)``.  Enables fluent
        chaining: ``ds.where(mask).reduce("z", reduction="integrate")``.

        Parameters
        ----------
        axis : str or tuple of str
            Surviving-axis name(s) to reduce away.
        **kwargs
            Forwarded to [`pypic.reductions.reduce`][pypic.reductions.reduce]:
            ``reduction``, ``selection``, ``fields``, ``nan_policy``.

        Returns
        -------
        FieldDataset
            With *axis* (or every name in the tuple) removed from the grid.

        Examples
        --------
        Integrating a unit field over two cells of unit spacing gives the
        trapezoidal path length between the two cell centres:

        >>> import numpy as np
        >>> from pypic.units import Normalization
        >>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0))
        >>> ds = FieldDataset.from_arrays(
        ...     {"B_1": np.ones((4, 3, 2))}, grid, Normalization.identity(),
        ... )
        >>> column = ds.reduce("z", reduction="integrate")
        >>> column.grid.dimensions, float(column["B_1"][0, 0])
        ((4, 3), 1.0)
        """
        from pypic.reductions import reduce as _reduce  # layered above dataset

        return _reduce(self, axis, **kwargs)

grid property

Grid metadata.

normalization property

Unit normalization.

species property

Species definitions.

physics property

Physics parameters.

metadata property

Arbitrary metadata (read-only view).

time property

Snapshot time in code units, or None when unknown.

metadata["time"] when the reader recorded one, otherwise metadata["step"] * grid.dt when both are known.

Examples:

>>> import numpy as np
>>> grid = GridInfo(dimensions=(2,), spacing=(1.0,), dt=0.5)
>>> FieldDataset.from_arrays(
...     {"B_1": np.zeros(2)}, grid, metadata={"step": 4}
... ).time
2.0
>>> FieldDataset.from_arrays({"B_1": np.zeros(2)}, grid).time is None
True

aliases property

Active field-name aliases (read-only view).

frame property

Current reference frame label.

transforms property

Registered frame transforms (read-only view).

available_frames property

Frame names reachable via registered transforms.

xr property

Raw xarray Dataset.

from_arrays(fields, grid, normalization=None, *, species=None, physics=None, metadata=None, aliases=None, frame='simulation', transforms=None, coords=None, strict_fields=True) classmethod

Build a FieldDataset from a dict of NumPy arrays.

Parameters:

Name	Type	Description	Default
fields	Mapping[str, FloatArray]	Mapping of field names to arrays. Shapes must match grid.dimensions. Names must resolve through the field registry (canonical names, registered aliases, or species-templated patterns) when strict_fields=True.	required
grid	GridInfo	Grid metadata.	required
normalization	Normalization | None	Unit normalization. Defaults to Normalization.undeclared() — arrays are left alone, but in_si and in_units raise on dimensional quantities rather than return code units labelled SI. Pass Normalization.identity() to assert the arrays already are SI.	None
species	Sequence[SpeciesInfo] | None	Species definitions, if applicable.	None
physics	PhysicsParams | None	Physics parameters.	None
metadata	Mapping[str, Any] | None	Arbitrary metadata.	None
aliases	Mapping[str, str] | None	Extra field-name aliases.	None
frame	str	Name of the reference frame these arrays are expressed in. Defaults to "simulation".	'simulation'
transforms	Mapping[str, FrameTransform] | None	Frame transforms reachable from frame, keyed by target-frame name. Consumed by transform_to.	None
coords	Mapping[str, FloatArray] | None	Coordinate arrays keyed by axis name, replacing the uniform ones derived from grid on those axes. For non-uniform meshes whose GridInfo spacing is only a mean; the true positions then live in the xarray coordinates.	None
strict_fields	bool	When True (default), every key in fields must resolve through pypic.fields.field_info — any unknown name raises KeyError with the full list of unresolved names. Set False to inject scratch fields (test fixtures, ad-hoc scalars) that aren't part of the canonical schema; the field is then stored without registry metadata.	True

Returns:

Type	Description
FieldDataset

Raises:

Type	Description
UnknownFieldError	When strict_fields=True and one or more keys in fields do not resolve through the field registry. Injection points fail loud on unknown names so reader bugs surface at construction time instead of later at compute(). Subclass of KeyError, so existing except KeyError callers keep working unchanged.
ValueError	When coords names an axis the grid does not have.

Examples:

>>> import numpy as np
>>> grid = GridInfo(dimensions=(2,), spacing=(1.0,))
>>> ds = FieldDataset.from_arrays({"B_1": np.array([1.0, 2.0])}, grid)
>>> ds["B_1"]
array([1., 2.])

Source code in src/pypic/dataset.py

@classmethod
def from_arrays(
    cls,
    fields: Mapping[str, FloatArray],
    grid: GridInfo,
    normalization: Normalization | None = None,
    *,
    species: Sequence[SpeciesInfo] | None = None,
    physics: PhysicsParams | None = None,
    metadata: Mapping[str, Any] | None = None,
    aliases: Mapping[str, str] | None = None,
    frame: str = "simulation",
    transforms: Mapping[str, FrameTransform] | None = None,
    coords: Mapping[str, FloatArray] | None = None,
    strict_fields: bool = True,
) -> FieldDataset:
    r"""Build a FieldDataset from a dict of NumPy arrays.

    Parameters
    ----------
    fields : Mapping[str, FloatArray]
        Mapping of field names to arrays. Shapes must match
        ``grid.dimensions``.  Names must resolve through the field
        registry (canonical names, registered aliases, or
        species-templated patterns) when ``strict_fields=True``.
    grid : GridInfo
        Grid metadata.
    normalization : Normalization | None
        Unit normalization.  Defaults to
        ``Normalization.undeclared()`` — arrays are left alone, but
        `in_si` and `in_units` raise on dimensional quantities
        rather than return code units labelled SI.  Pass
        ``Normalization.identity()`` to assert the arrays already
        *are* SI.
    species : Sequence[SpeciesInfo] | None
        Species definitions, if applicable.
    physics : PhysicsParams | None
        Physics parameters.
    metadata : Mapping[str, Any] | None
        Arbitrary metadata.
    aliases : Mapping[str, str] | None
        Extra field-name aliases.
    frame : str
        Name of the reference frame these arrays are expressed in.
        Defaults to ``"simulation"``.
    transforms : Mapping[str, FrameTransform] | None
        Frame transforms reachable from *frame*, keyed by
        target-frame name.  Consumed by `transform_to`.
    coords : Mapping[str, FloatArray] | None
        Coordinate arrays keyed by axis name, replacing the uniform
        ones derived from *grid* on those axes.  For non-uniform
        meshes whose ``GridInfo`` spacing is only a mean; the true
        positions then live in the xarray coordinates.
    strict_fields : bool
        When ``True`` (default), every key in *fields* must resolve
        through ``pypic.fields.field_info`` — any unknown name
        raises ``KeyError`` with the full list of unresolved names.
        Set ``False`` to inject scratch fields (test fixtures, ad-hoc
        scalars) that aren't part of the canonical schema; the field
        is then stored without registry metadata.

    Returns
    -------
    FieldDataset

    Raises
    ------
    UnknownFieldError
        When ``strict_fields=True`` and one or more keys in *fields*
        do not resolve through the field registry.  Injection
        points fail loud on unknown names so reader bugs surface at
        construction time instead of later at ``compute()``.  Subclass of
        `KeyError`, so existing ``except KeyError`` callers
        keep working unchanged.
    ValueError
        When *coords* names an axis the grid does not have.

    Examples
    --------
    >>> import numpy as np
    >>> grid = GridInfo(dimensions=(2,), spacing=(1.0,))
    >>> ds = FieldDataset.from_arrays({"B_1": np.array([1.0, 2.0])}, grid)
    >>> ds["B_1"]
    array([1., 2.])
    """
    if normalization is None:
        normalization = Normalization.undeclared()
    dim_names = list(grid.surviving_axis_names)
    coord_arrays = grid.coordinate_arrays()
    axis_coords = dict(zip(dim_names, coord_arrays, strict=True))
    if coords:
        if unknown := sorted(set(coords) - set(dim_names)):
            msg = f"coords name axes the grid lacks: {unknown}; axes: {dim_names}"
            raise ValueError(msg)
        axis_coords.update(coords)

    data_vars: dict[str, xr.DataArray] = {}
    unresolved: list[str] = []
    for var_name, arr in fields.items():
        da = xr.DataArray(data=arr, dims=dim_names)
        try:
            info = _field_info(var_name, axis_names=grid.geometry.axis_names)
            da.attrs["long_name"] = info.long_name
            da.attrs["units"] = "normalized"
            da.attrs["quantity_type"] = info.quantity_type
            da.attrs["si_unit"] = info.si_unit
            if info.latex:
                da.attrs["latex"] = info.latex
            ud = info.unit_dimension or quantity_dimension(info.quantity_type)
            da.attrs["unit_dimension"] = list(ud)
        except KeyError:
            if strict_fields:
                unresolved.append(var_name)
        data_vars[var_name] = da
    if unresolved:
        msg = (
            f"Field names not in the registry: {unresolved!r}. "
            "Pass strict_fields=False to inject scratch fields, or "
            "use FieldDataset.with_field(name, data, quantity_type=...) "
            "to register an ad-hoc quantity_type."
        )
        raise UnknownFieldError(msg)
    dataset = xr.Dataset(data_vars, coords=axis_coords)
    return cls(
        dataset,
        grid,
        normalization,
        species=species,
        physics=physics,
        metadata=metadata,
        aliases=aliases,
        frame=frame,
        transforms=transforms,
    )

transform_to(target)

Transform this dataset to a different coordinate reference frame.

Applies an affine transformation: translates the grid origin, scales coordinates, and rotates vector field components and pressure tensors. Scalar values pass through unchanged. When the rotation includes an axis swap or reflection, field arrays are copied to a contiguous buffer in the new axis order; pure translations and identity rotations only touch metadata.

Parameters:

Name	Type	Description	Default
target	str or FrameTransform	Target frame name (looked up in the transform registry) or a FrameTransform instance applied directly.	required

Returns:

Type	Description
FieldDataset	New dataset in the target frame.

Raises:

Type	Description
ValueError	If no transform path exists from the current frame.
KeyError	If target is a string and no transforms are registered.

Examples:

A -90° rotation about \(z\) sends \(\hat{x}\) to \(-\hat{y}\), and swaps the two grid axes with it:

>>> import numpy as np
>>> from pypic.coordinates import FrameTransform
>>> from pypic.units import Normalization
>>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0))
>>> gsm = FrameTransform(
...     "simulation", "GSM",
...     rotation=((0.0, 1.0, 0.0), (-1.0, 0.0, 0.0), (0.0, 0.0, 1.0)),
... )
>>> ds = FieldDataset.from_arrays(
...     {"B_1": np.ones((4, 3, 2)),
...      "B_2": np.zeros((4, 3, 2)),
...      "B_3": np.zeros((4, 3, 2))},
...     grid, Normalization.identity(), transforms={"GSM": gsm},
... )
>>> rotated = ds.transform_to("GSM")
>>> rotated.frame, rotated.grid.dimensions
('GSM', (3, 4, 2))
>>> float(rotated["B_2"][0, 0, 0])
-1.0

Source code in src/pypic/dataset.py

def transform_to(self, target: str | FrameTransform) -> FieldDataset:
    r"""Transform this dataset to a different coordinate reference frame.

    Applies an affine transformation: translates the grid origin,
    scales coordinates, and rotates vector field components and
    pressure tensors. Scalar values pass through unchanged. When the
    rotation includes an axis swap or reflection, field arrays are
    copied to a contiguous buffer in the new axis order; pure
    translations and identity rotations only touch metadata.

    Parameters
    ----------
    target : str or FrameTransform
        Target frame name (looked up in the transform registry) or
        a `FrameTransform` instance applied directly.

    Returns
    -------
    FieldDataset
        New dataset in the target frame.

    Raises
    ------
    ValueError
        If no transform path exists from the current frame.
    KeyError
        If *target* is a string and no transforms are registered.

    Examples
    --------
    A -90° rotation about $z$ sends $\hat{x}$ to $-\hat{y}$, and
    swaps the two grid axes with it:

    >>> import numpy as np
    >>> from pypic.coordinates import FrameTransform
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0))
    >>> gsm = FrameTransform(
    ...     "simulation", "GSM",
    ...     rotation=((0.0, 1.0, 0.0), (-1.0, 0.0, 0.0), (0.0, 0.0, 1.0)),
    ... )
    >>> ds = FieldDataset.from_arrays(
    ...     {"B_1": np.ones((4, 3, 2)),
    ...      "B_2": np.zeros((4, 3, 2)),
    ...      "B_3": np.zeros((4, 3, 2))},
    ...     grid, Normalization.identity(), transforms={"GSM": gsm},
    ... )
    >>> rotated = ds.transform_to("GSM")
    >>> rotated.frame, rotated.grid.dimensions
    ('GSM', (3, 4, 2))
    >>> float(rotated["B_2"][0, 0, 0])
    -1.0
    """
    if isinstance(target, FrameTransform):
        transform = target
        target_frame = transform.target_frame
    else:
        target_frame = target
        if target_frame == self._frame:
            return self
        if not self._transforms:
            msg = (
                f"No transforms registered; cannot transform "
                f"from {self._frame!r} to {target_frame!r}"
            )
            raise KeyError(msg)
        transform = resolve_transform(self._frame, target_frame, self._transforms)
    rotation = transform.rotation_matrix

    # Determine which axes survive and derive the axis permutation
    # and sign map from the rotation matrix up front. All field
    # arrays then flow through a single reorient → construct pass.
    if self._grid.surviving_axes is not None:
        indices = list(self._grid.surviving_axes)
    else:
        indices = list(range(len(self._grid.dimensions)))
    ndim = len(indices)

    # The grid transformation (transpose + flip) only works for
    # signed permutation matrices (axis swaps and reflections).
    # General rotations would require interpolation onto a new grid.
    rotation_sub = rotation[np.ix_(indices, indices)]
    abs_rotation = np.abs(rotation_sub)
    if not (
        np.allclose(np.sum(abs_rotation, axis=1), 1.0, atol=1e-6)
        and np.allclose(np.sum(abs_rotation, axis=0), 1.0, atol=1e-6)
    ):
        # GeometryUnsupportedError subclasses NotImplementedError, so
        # existing handlers still catch it, and the server maps it to
        # 400 rather than a 500 "internal" frame.
        raise GeometryUnsupportedError(
            "transform_to() only supports axis-swap/reflection "
            "rotation matrices (signed permutations). General "
            "rotations require grid interpolation (not yet implemented)."
        )

    # Axis permutation (target_i → source local index) and sign.
    axis_permutation: list[int] = []
    axis_signs: list[float] = []
    for target_i in indices:
        row = rotation[target_i, :]
        source_orig = int(np.argmax(np.abs(row[indices])))
        axis_permutation.append(source_orig)
        axis_signs.append(float(np.sign(row[indices[source_orig]])))

    needs_transpose = axis_permutation != list(range(ndim))
    flip_axes = [i for i, s in enumerate(axis_signs) if s < 0]

    def _reorient(arr: FloatArray) -> FloatArray:
        if needs_transpose:
            arr = np.transpose(arr, axis_permutation)
        for ax in flip_axes:
            arr = np.flip(arr, axis=ax)
        if needs_transpose or flip_axes:
            arr = np.ascontiguousarray(arr)
        return arr

    # Accumulate raw numpy arrays — no intermediate DataArray wrapping.
    new_arrays: dict[str, FloatArray] = {}
    rotated: set[str] = set()

    for n1, n2, n3 in find_vector_triplets(self.field_names()):
        r1, r2, r3 = rotate_vector_components(
            self[n1], self[n2], self[n3], rotation
        )
        for name, arr in ((n1, r1), (n2, r2), (n3, r3)):
            new_arrays[name] = _reorient(arr)
            rotated.add(name)

    for p11, p22, p33, p12, p13, p23 in find_pressure_tensor_groups(
        self.field_names()
    ):
        rp = rotate_pressure_tensor(
            self[p11],
            self[p22],
            self[p33],
            self[p12],
            self[p13],
            self[p23],
            rotation,
        )
        for name, arr in zip((p11, p22, p33, p12, p13, p23), rp, strict=True):
            new_arrays[name] = _reorient(arr)
            rotated.add(name)

    for raw_name in self._ds.data_vars:
        name = str(raw_name)
        if name not in rotated:
            new_arrays[name] = _reorient(self._ds[name].values)

    # Compute new origin, spacing, dimensions from permuted source
    transform_origin = np.array(transform.origin, dtype=np.float64)
    dx_scale = transform.scale
    old_coords = self._grid.coordinate_arrays()
    new_origin_list: list[float] = []
    new_spacing_list: list[float] = []
    new_dims_list: list[int] = []
    for src_i, sign in zip(axis_permutation, axis_signs, strict=True):
        src_coords = old_coords[src_i]
        coord_first = (
            dx_scale * sign * (src_coords[0] - transform_origin[indices[src_i]])
        )
        coord_last = (
            dx_scale * sign * (src_coords[-1] - transform_origin[indices[src_i]])
        )
        dx = dx_scale * self._grid.spacing[src_i]
        new_origin_list.append(float(min(coord_first, coord_last)) - 0.5 * dx)
        new_spacing_list.append(dx)
        new_dims_list.append(self._grid.dimensions[src_i])

    new_surviving = None
    if self._grid.surviving_axes is not None:
        new_surviving = tuple(
            self._grid.surviving_axes[axis_permutation[i]] for i in range(ndim)
        )

    new_geometry = self._grid.geometry
    if transform.target_axis_names is not None:
        new_geometry = copy.replace(
            new_geometry, axis_names=transform.target_axis_names
        )

    new_grid = copy.replace(
        self._grid,
        dimensions=tuple(new_dims_list),
        origin=tuple(new_origin_list),
        spacing=tuple(new_spacing_list),
        geometry=new_geometry,
        surviving_axes=new_surviving,
    )

    # Single Dataset construction: each field is wrapped exactly once.
    new_dim_names = list(new_grid.surviving_axis_names)
    coord_arrays = new_grid.coordinate_arrays()
    coords = {new_dim_names[i]: coord_arrays[i] for i in range(len(new_dim_names))}
    new_ds = xr.Dataset(
        data_vars={
            name: (new_dim_names, arr, dict(self._ds[name].attrs))
            for name, arr in new_arrays.items()
        },
        coords=coords,
    )

    return FieldDataset(
        new_ds,
        new_grid,
        self._normalization,
        species=self._species,
        physics=self._physics,
        metadata=self._metadata,
        frame=target_frame,
        transforms=self._transforms,
    )

resolve_key(key)

Resolve a field key through aliases to the canonical name.

Returns the canonical (data-vars) name for key, which may itself be a canonical name or an alias from this dataset's alias table. Raises KeyError (with close-match suggestions) if key matches neither.

Parameters:

Name	Type	Description	Default
key	str	Canonical field name or alias.	required

Returns:

Type	Description
str	Canonical name in self._ds.data_vars.

Examples:

>>> import numpy as np
>>> from pypic.units import Normalization
>>> grid = GridInfo(
...     dimensions=(2,), spacing=(1.0,), origin=(0.0,),
...     geometry=CARTESIAN,
... )
>>> ds = FieldDataset.from_arrays(
...     {"B_1": np.ones(2)}, grid, Normalization.identity()
... )
>>> ds.resolve_key("Bx")  # Cartesian alias
'B_1'
>>> ds.resolve_key("B_1")  # canonical
'B_1'

Source code in src/pypic/dataset.py

def resolve_key(self, key: str) -> str:
    """Resolve a field key through aliases to the canonical name.

    Returns the canonical (data-vars) name for *key*, which may itself
    be a canonical name or an alias from this dataset's alias table.
    Raises `KeyError` (with close-match suggestions) if *key*
    matches neither.

    Parameters
    ----------
    key : str
        Canonical field name or alias.

    Returns
    -------
    str
        Canonical name in ``self._ds.data_vars``.

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(
    ...     dimensions=(2,), spacing=(1.0,), origin=(0.0,),
    ...     geometry=CARTESIAN,
    ... )
    >>> ds = FieldDataset.from_arrays(
    ...     {"B_1": np.ones(2)}, grid, Normalization.identity()
    ... )
    >>> ds.resolve_key("Bx")  # Cartesian alias
    'B_1'
    >>> ds.resolve_key("B_1")  # canonical
    'B_1'
    """
    if key in self._ds.data_vars:
        return key
    canonical = self._aliases.get(key)
    if canonical is not None:
        return canonical
    available = sorted(self._ds.data_vars, key=str)
    alias_keys = sorted(self._aliases)
    msg = f"Field {key!r} not found. Available: {available}. Aliases: {alias_keys}."
    import difflib

    candidates = [str(v) for v in self._ds.data_vars] + list(self._aliases)
    suggestions = difflib.get_close_matches(key, candidates, n=3, cutoff=0.5)
    if suggestions:
        msg += f" Did you mean: {suggestions}?"
    raise UnknownFieldError(msg)

__getitem__(key)

Return field data as a NumPy array (zero-copy when possible).

Parameters:

Name	Type	Description	Default
key	str	Canonical field name or alias.	required

Returns:

Type	Description
NDArray

Source code in src/pypic/dataset.py

def __getitem__(self, key: str) -> FloatArray:
    """Return field data as a NumPy array (zero-copy when possible).

    Parameters
    ----------
    key : str
        Canonical field name or alias.

    Returns
    -------
    NDArray
    """
    resolved = self.resolve_key(key)
    return self._ds[resolved].values

has_field(key)

Check whether a field exists (canonical or alias).

Parameters:

Name	Type	Description	Default
key	str	Field name to check.	required

Returns:

Type	Description
bool

Examples:

>>> import numpy as np
>>> from pypic.units import Normalization
>>> grid = GridInfo(
...     dimensions=(2,), spacing=(1.0,), origin=(0.0,),
...     geometry=CARTESIAN,
... )
>>> ds = FieldDataset.from_arrays(
...     {"B_1": np.array([1.0, 2.0])}, grid, Normalization.identity(),
... )
>>> ds.has_field("B_1"), ds.has_field("Bx"), ds.has_field("rho")
(True, True, False)

Source code in src/pypic/dataset.py

def has_field(self, key: str) -> bool:
    """Check whether a field exists (canonical or alias).

    Parameters
    ----------
    key : str
        Field name to check.

    Returns
    -------
    bool

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(
    ...     dimensions=(2,), spacing=(1.0,), origin=(0.0,),
    ...     geometry=CARTESIAN,
    ... )
    >>> ds = FieldDataset.from_arrays(
    ...     {"B_1": np.array([1.0, 2.0])}, grid, Normalization.identity(),
    ... )
    >>> ds.has_field("B_1"), ds.has_field("Bx"), ds.has_field("rho")
    (True, True, False)
    """
    return key in self._ds.data_vars or key in self._aliases

field_names()

Return canonical field names (no aliases).

Returns:

Type	Description
list[str]

Examples:

>>> import numpy as np
>>> from pypic.units import Normalization
>>> grid = GridInfo(
...     dimensions=(2,), spacing=(1.0,), origin=(0.0,),
...     geometry=CARTESIAN,
... )
>>> ds = FieldDataset.from_arrays(
...     {"B_1": np.array([1.0, 2.0]), "rho_c": np.array([0.5, 0.5])},
...     grid, Normalization.identity(),
... )
>>> sorted(ds.field_names())
['B_1', 'rho_c']

Source code in src/pypic/dataset.py

def field_names(self) -> list[str]:
    """Return canonical field names (no aliases).

    Returns
    -------
    list[str]

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(
    ...     dimensions=(2,), spacing=(1.0,), origin=(0.0,),
    ...     geometry=CARTESIAN,
    ... )
    >>> ds = FieldDataset.from_arrays(
    ...     {"B_1": np.array([1.0, 2.0]), "rho_c": np.array([0.5, 0.5])},
    ...     grid, Normalization.identity(),
    ... )
    >>> sorted(ds.field_names())
    ['B_1', 'rho_c']
    """
    return list(self._ds.data_vars)  # type: ignore[arg-type]  # xarray types Hashable, always str

select_fields(names)

Return a new FieldDataset containing only the specified fields.

Parameters:

Name	Type	Description	Default
names	Iterable[str]	Field names to keep (canonical or alias).	required

Returns:

Type	Description
FieldDataset

Raises:

Type	Description
KeyError	If any name is not found.

Examples:

>>> import numpy as np
>>> from pypic.units import Normalization
>>> grid = GridInfo(
...     dimensions=(2,), spacing=(1.0,), origin=(0.0,),
...     geometry=CARTESIAN,
... )
>>> ds = FieldDataset.from_arrays(
...     {"B_1": np.array([1.0, 2.0]), "B_2": np.array([3.0, 4.0]),
...      "rho_c": np.array([0.5, 0.5])},
...     grid, Normalization.identity(),
... )
>>> sub = ds.select_fields(["rho_c", "B_1"])
>>> sub.field_names()
['rho_c', 'B_1']

Source code in src/pypic/dataset.py

def select_fields(self, names: Iterable[str]) -> FieldDataset:
    """Return a new FieldDataset containing only the specified fields.

    Parameters
    ----------
    names : Iterable[str]
        Field names to keep (canonical or alias).

    Returns
    -------
    FieldDataset

    Raises
    ------
    KeyError
        If any name is not found.

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(
    ...     dimensions=(2,), spacing=(1.0,), origin=(0.0,),
    ...     geometry=CARTESIAN,
    ... )
    >>> ds = FieldDataset.from_arrays(
    ...     {"B_1": np.array([1.0, 2.0]), "B_2": np.array([3.0, 4.0]),
    ...      "rho_c": np.array([0.5, 0.5])},
    ...     grid, Normalization.identity(),
    ... )
    >>> sub = ds.select_fields(["rho_c", "B_1"])
    >>> sub.field_names()
    ['rho_c', 'B_1']
    """
    # Dict-as-ordered-set: preserves caller insertion order while
    # deduplicating on the resolved canonical name. Alphabetical
    # sorting hides both request order and dataset order from the
    # caller, so we keep the order the caller asked for.
    resolved: dict[str, None] = {}
    for name in names:
        resolved[self.resolve_key(name)] = None

    new_ds = self._ds[list(resolved)]
    return FieldDataset(
        new_ds,
        self._grid,
        self._normalization,
        species=self._species,
        physics=self._physics,
        metadata=self._metadata,
        aliases={k: v for k, v in self._aliases.items() if v in resolved},
        frame=self._frame,
        transforms=self._transforms,
    )

sel(indexers=None, *, method=None, **kwargs)

Label-based selection, returning a new FieldDataset.

Accepts indexers as a dict (useful for Unicode dim names like "θ" that can't be keyword arguments) and/or **kwargs.

Parameters:

Name	Type	Description	Default
indexers	dict[str, Any] | None	Dimension-name → label mapping.	None
method	str | None	Passed to xr.Dataset.sel (e.g. "nearest").	None
**kwargs	Any	Additional dimension selections.	{}

Returns:

Type	Description
FieldDataset

Examples:

Cell centres sit at half-integer positions, so a label that falls between them needs method="nearest":

>>> import numpy as np
>>> from pypic.units import Normalization
>>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0))
>>> ds = FieldDataset.from_arrays(
...     {"B_1": np.arange(24.0).reshape(4, 3, 2)},
...     grid, Normalization.identity(),
... )
>>> midplane = ds.sel({"z": 1.4}, method="nearest")
>>> midplane.grid.dimensions, midplane.grid.surviving_axis_names
((4, 3), ('x', 'y'))

Source code in src/pypic/dataset.py

def sel(
    self,
    indexers: dict[str, Any] | None = None,
    *,
    method: str | None = None,
    **kwargs: Any,  # noqa: ANN401 — xarray passthrough
) -> FieldDataset:
    """Label-based selection, returning a new FieldDataset.

    Accepts ``indexers`` as a dict (useful for Unicode dim names like
    ``"θ"`` that can't be keyword arguments) and/or ``**kwargs``.

    Parameters
    ----------
    indexers : dict[str, Any] | None
        Dimension-name → label mapping.
    method : str | None
        Passed to ``xr.Dataset.sel`` (e.g. ``"nearest"``).
    **kwargs
        Additional dimension selections.

    Returns
    -------
    FieldDataset

    Examples
    --------
    Cell centres sit at half-integer positions, so a label that
    falls between them needs ``method="nearest"``:

    >>> import numpy as np
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0))
    >>> ds = FieldDataset.from_arrays(
    ...     {"B_1": np.arange(24.0).reshape(4, 3, 2)},
    ...     grid, Normalization.identity(),
    ... )
    >>> midplane = ds.sel({"z": 1.4}, method="nearest")
    >>> midplane.grid.dimensions, midplane.grid.surviving_axis_names
    ((4, 3), ('x', 'y'))
    """
    merged = dict(indexers) if indexers else {}
    merged.update(kwargs)
    return self._wrap_sliced(self._ds.sel(merged, method=method))

compute(name)

Compute a derived quantity by name. Returns code units.

Parameters:

Name	Type	Description	Default
name	str	Derived quantity name (e.g. "|B|", "beta", "v_A"). See pypic.compute.available_quantities() for the full list.	required

Returns:

Type	Description
FloatArray	Computed array in code units.

Examples:

>>> import numpy as np
>>> from pypic.units import Normalization
>>> grid = GridInfo(
...     dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0),
...     geometry=CARTESIAN,
... )
>>> ds = FieldDataset.from_arrays(
...     {"B_1": np.full((4,3,2), 3.0),
...      "B_2": np.full((4,3,2), 4.0),
...      "B_3": np.zeros((4,3,2))},
...     grid, Normalization.identity(),
... )
>>> ds.compute("|B|")[0, 0, 0]
np.float64(5.0)

Source code in src/pypic/dataset.py

def compute(self, name: str) -> FloatArray:
    """Compute a derived quantity by name. Returns code units.

    Parameters
    ----------
    name : str
        Derived quantity name (e.g. ``"|B|"``, ``"beta"``, ``"v_A"``).
        See ``pypic.compute.available_quantities()`` for the full list.

    Returns
    -------
    FloatArray
        Computed array in code units.

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(
    ...     dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0),
    ...     geometry=CARTESIAN,
    ... )
    >>> ds = FieldDataset.from_arrays(
    ...     {"B_1": np.full((4,3,2), 3.0),
    ...      "B_2": np.full((4,3,2), 4.0),
    ...      "B_3": np.zeros((4,3,2))},
    ...     grid, Normalization.identity(),
    ... )
    >>> ds.compute("|B|")[0, 0, 0]
    np.float64(5.0)
    """
    from pypic.compute import compute_field  # layered above dataset

    return compute_field(name, self)

with_field(name, data, quantity_type=None, *, long_name='', latex='')

Return a new FieldDataset with an additional custom field.

The field's quantity_type is stored in xarray attrs, so in_si(), field_info(), and unit conversion work without global register_field() calls.

When quantity_type is None, metadata is looked up from the field registry automatically. For unregistered fields, provide quantity_type explicitly.

Parameters:

Name	Type	Description	Default
name	str	Field name.	required
data	FloatArray	Array matching the grid dimensions.	required
quantity_type	QuantityType | str | None	Physical quantity type (e.g. QuantityType.VELOCITY). When None, auto-filled from the field registry.	None
long_name	str	Human-readable label for plot titles.	''
latex	str	LaTeX symbol for plot labels.	''

Returns:

Type	Description
FieldDataset	New dataset with the field added.

Raises:

Type	Description
ValueError	If quantity_type is not recognized, or is None and the field name is not in the registry.

Source code in src/pypic/dataset.py

def with_field(
    self,
    name: str,
    data: FloatArray,
    quantity_type: QuantityType | str | None = None,
    *,
    long_name: str = "",
    latex: str = "",
) -> FieldDataset:
    """Return a new FieldDataset with an additional custom field.

    The field's ``quantity_type`` is stored in xarray attrs, so
    ``in_si()``, ``field_info()``, and unit conversion work without
    global ``register_field()`` calls.

    When *quantity_type* is ``None``, metadata is looked up from the
    field registry automatically. For unregistered fields, provide
    *quantity_type* explicitly.

    Parameters
    ----------
    name : str
        Field name.
    data : FloatArray
        Array matching the grid dimensions.
    quantity_type : QuantityType | str | None
        Physical quantity type (e.g. ``QuantityType.VELOCITY``).
        When ``None``, auto-filled from the field registry.
    long_name : str
        Human-readable label for plot titles.
    latex : str
        LaTeX symbol for plot labels.

    Returns
    -------
    FieldDataset
        New dataset with the field added.

    Raises
    ------
    ValueError
        If *quantity_type* is not recognized, or is ``None`` and
        the field name is not in the registry.
    """
    expected = tuple(self._grid.dimensions)
    if data.shape != expected:
        msg = f"Array shape {data.shape} doesn't match grid dimensions {expected}"
        raise ValueError(msg)

    info_ud: tuple[int, int, int, int, int, int, int] | None = None
    if quantity_type is None:
        try:
            info = _field_info(name, axis_names=self._grid.geometry.axis_names)
        except KeyError:
            msg = f"Unknown field {name!r}; provide quantity_type explicitly"
            raise ValueError(msg) from None
        qt = info.quantity_type
        if not long_name:
            long_name = info.long_name
        if not latex:
            latex = info.latex
        info_ud = info.unit_dimension
    else:
        qt = str(quantity_type)

    if qt not in _QUANTITY_UNITS:
        valid = sorted(_QUANTITY_UNITS)
        msg = f"Unknown quantity_type {qt!r}. Valid: {valid}"
        raise ValueError(msg)

    si_unit = _QUANTITY_UNITS[qt]
    dim_names = list(self._grid.surviving_axis_names)
    da = xr.DataArray(data=data, dims=dim_names)
    da.attrs["quantity_type"] = qt
    da.attrs["si_unit"] = si_unit
    da.attrs["units"] = "normalized"
    da.attrs["long_name"] = long_name
    da.attrs["latex"] = latex
    da.attrs["unit_dimension"] = list(info_ud or quantity_dimension(qt))

    new_ds = self._ds.assign({name: da})
    return FieldDataset(
        new_ds,
        self._grid,
        self._normalization,
        species=self._species,
        physics=self._physics,
        metadata=self._metadata,
        aliases=dict(self._aliases),
        frame=self._frame,
        transforms=self._transforms,
    )

with_derived(*names)

Return a new dataset with derived fields computed and stored.

For each name, computes the field (if not already present), looks up metadata from the field registry, and attaches it with full metadata. Fields already in the dataset are skipped.

When a vector-component recipe is encountered (e.g. "S_1" from Poynting flux), all sibling components ("S_2", "S_3") are computed from a single function call and stored together.

Parameters:

Name	Type	Description	Default
*names	str	Derived quantity names (e.g. "|B|", "beta").	()

Returns:

Type	Description
FieldDataset	New dataset with the computed fields attached.

Examples:

>>> import numpy as np
>>> from pypic.units import Normalization
>>> grid = GridInfo(
...     dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0),
...     geometry=CARTESIAN,
... )
>>> ds = FieldDataset.from_arrays(
...     {"B_1": np.full((4,3,2), 3.0),
...      "B_2": np.full((4,3,2), 4.0),
...      "B_3": np.zeros((4,3,2))},
...     grid, Normalization.identity(),
... )
>>> ds = ds.with_derived("|B|")
>>> ds.has_field("|B|")
True
>>> ds["|B|"][0, 0, 0]
np.float64(5.0)

Source code in src/pypic/dataset.py

def with_derived(self, *names: str) -> FieldDataset:
    """Return a new dataset with derived fields computed and stored.

    For each name, computes the field (if not already present),
    looks up metadata from the field registry, and attaches it
    with full metadata. Fields already in the dataset are skipped.

    When a vector-component recipe is encountered (e.g. ``"S_1"``
    from Poynting flux), all sibling components (``"S_2"``, ``"S_3"``)
    are computed from a single function call and stored together.

    Parameters
    ----------
    *names : str
        Derived quantity names (e.g. ``"|B|"``, ``"beta"``).

    Returns
    -------
    FieldDataset
        New dataset with the computed fields attached.

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(
    ...     dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0),
    ...     geometry=CARTESIAN,
    ... )
    >>> ds = FieldDataset.from_arrays(
    ...     {"B_1": np.full((4,3,2), 3.0),
    ...      "B_2": np.full((4,3,2), 4.0),
    ...      "B_3": np.zeros((4,3,2))},
    ...     grid, Normalization.identity(),
    ... )
    >>> ds = ds.with_derived("|B|")
    >>> ds.has_field("|B|")
    True
    >>> ds["|B|"][0, 0, 0]
    np.float64(5.0)
    """
    from pypic.compute import compute_with_siblings  # layered above dataset

    result = self
    for name in names:
        if result.has_field(name):
            continue
        for field, data in compute_with_siblings(name, result).items():
            if not result.has_field(field):
                result = result.with_field(field, data)
    return result

field_info(name)

Return metadata for a field or derived quantity.

Checks xarray DataArray attrs first (set by with_field() or from_arrays()), then falls back to the global registry.

Parameters:

Name	Type	Description	Default
name	str	Field or derived quantity name (canonical or alias).	required

Returns:

Type	Description
FieldInfo

Source code in src/pypic/dataset.py

def field_info(self, name: str) -> FieldInfo:
    """Return metadata for a field or derived quantity.

    Checks xarray DataArray attrs first (set by ``with_field()``
    or ``from_arrays()``), then falls back to the global registry.

    Parameters
    ----------
    name : str
        Field or derived quantity name (canonical or alias).

    Returns
    -------
    FieldInfo
    """
    if self.has_field(name):
        resolved = self.resolve_key(name)
        attrs = self._ds[resolved].attrs
        qt = attrs.get("quantity_type")
        if qt is not None:
            return FieldInfo(
                quantity_type=qt,
                long_name=attrs.get("long_name", ""),
                si_unit=attrs.get("si_unit", ""),
                latex=attrs.get("latex", ""),
            )
    return _field_info(name, axis_names=self._grid.geometry.axis_names)

in_si(name)

Return a field or derived quantity in SI units.

Checks xarray DataArray attrs first (set by with_field() or from_arrays()), then falls back to the global registry.

For fields that have been reduced via pypic.reductions.reduce with reduction="integrate", the attrs["reduction"]["length_axes"] provenance stamp records how many length-dimension factors the integration added; in_si multiplies the registry SI factor by normalization.length_ref ** length_axes so column densities, line-of-sight integrals, etc. come out in the right SI units (e.g. column density m\ :sup:-2 instead of m\ :sup:-3). Weighted integrate does not stamp length_axes because the length factor cancels in ∫ f w dx / ∫ w dx.

Parameters:

Name	Type	Description	Default
name	str	Field or derived quantity name.	required

Returns:

Type	Description
FloatArray	Values in SI units.

Raises:

Type	Description
ValueError	If no SI conversion is registered for name — a field carrying no quantity_type attr whose name the global registry also cannot resolve.

Examples:

>>> import numpy as np
>>> from pypic.units import Normalization
>>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0))
>>> norm = Normalization.mhd_standard(6.371e6, 1.67e-17, 5.0e-9)
>>> ds = FieldDataset.from_arrays(
...     {"B_1": np.full((4, 3, 2), 2.0)}, grid, norm,
... )
>>> float(ds.in_si("B_1")[0, 0, 0])  # 2 code units at B_ref = 5 nT
1e-08

Source code in src/pypic/dataset.py

def in_si(self, name: str) -> FloatArray:
    r"""Return a field or derived quantity in SI units.

    Checks xarray DataArray attrs first (set by ``with_field()``
    or ``from_arrays()``), then falls back to the global registry.

    For fields that have been reduced via
    [`pypic.reductions.reduce`][pypic.reductions.reduce] with
    ``reduction="integrate"``,
    the ``attrs["reduction"]["length_axes"]`` provenance stamp
    records how many length-dimension factors the integration
    added; ``in_si`` multiplies the registry SI factor by
    ``normalization.length_ref ** length_axes`` so column
    densities, line-of-sight integrals, etc. come out in the
    right SI units (e.g. column density m\ :sup:`-2` instead of
    m\ :sup:`-3`).  Weighted ``integrate`` does not stamp
    ``length_axes`` because the length factor cancels in
    ``∫ f w dx / ∫ w dx``.

    Parameters
    ----------
    name : str
        Field or derived quantity name.

    Returns
    -------
    FloatArray
        Values in SI units.

    Raises
    ------
    ValueError
        If no SI conversion is registered for *name* — a field
        carrying no ``quantity_type`` attr whose name the global
        registry also cannot resolve.

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0))
    >>> norm = Normalization.mhd_standard(6.371e6, 1.67e-17, 5.0e-9)
    >>> ds = FieldDataset.from_arrays(
    ...     {"B_1": np.full((4, 3, 2), 2.0)}, grid, norm,
    ... )
    >>> float(ds.in_si("B_1")[0, 0, 0])  # 2 code units at B_ref = 5 nT
    1e-08
    """
    from pypic.compute import compute_field, field_si_factor  # above dataset

    length_axes = 0
    factor: float | None = None
    if self.has_field(name):
        resolved = self.resolve_key(name)
        data = self._ds[resolved].values
        reduction_attr = self._ds[resolved].attrs.get("reduction") or {}
        length_axes = int(reduction_attr.get("length_axes", 0))
        qt = self._ds[resolved].attrs.get("quantity_type")
        if qt is not None:
            factor = self._normalization.si_factor(qt)
    else:
        data = compute_field(name, self)
    if factor is None:
        # Either a derived quantity, or a stored field with no
        # ``quantity_type`` attr — what
        # ``from_arrays(..., strict_fields=False)`` and
        # ``open_virtual`` produce. Both resolve through the registry.
        factor = field_si_factor(name, self._normalization)
    if length_axes:
        factor *= self._normalization.length_ref**length_axes
    return data if factor == 1.0 else data * factor

in_units(name, unit_str)

Return a field or derived quantity in display units.

Parameters:

Name	Type	Description	Default
name	str	Field or derived quantity name.	required
unit_str	str	Target unit (e.g. "nT", "km/s", "cm^-3").	required

Returns:

Type	Description
FloatArray	Values in the requested units.

Notes

Temperature is stored in energy units (J), so use in_units(name, "eV") (or "keV", "MeV") for the plasma working unit, or in_units(name, "K") for the Boltzmann-factor-converted form. Magnetic field defaults to T in SI; in_units(name, "nT") is the space-physics idiom. See _DISPLAY_UNITS in pypic.units for the full vocabulary.

Examples:

>>> import numpy as np
>>> from pypic.units import Normalization
>>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0))
>>> norm = Normalization.mhd_standard(6.371e6, 1.67e-17, 5.0e-9)
>>> ds = FieldDataset.from_arrays(
...     {"B_1": np.full((4, 3, 2), 2.0)}, grid, norm,
... )
>>> float(ds.in_units("B_1", "nT")[0, 0, 0])
10.0

Source code in src/pypic/dataset.py

def in_units(self, name: str, unit_str: str) -> FloatArray:
    """Return a field or derived quantity in display units.

    Parameters
    ----------
    name : str
        Field or derived quantity name.
    unit_str : str
        Target unit (e.g. ``"nT"``, ``"km/s"``, ``"cm^-3"``).

    Returns
    -------
    FloatArray
        Values in the requested units.

    Notes
    -----
    Temperature is stored in **energy units** (J), so use
    ``in_units(name, "eV")`` (or ``"keV"``, ``"MeV"``) for the
    plasma working unit, or ``in_units(name, "K")`` for the
    Boltzmann-factor-converted form. Magnetic field defaults to
    T in SI; ``in_units(name, "nT")`` is the space-physics
    idiom. See ``_DISPLAY_UNITS`` in ``pypic.units`` for the
    full vocabulary.

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0))
    >>> norm = Normalization.mhd_standard(6.371e6, 1.67e-17, 5.0e-9)
    >>> ds = FieldDataset.from_arrays(
    ...     {"B_1": np.full((4, 3, 2), 2.0)}, grid, norm,
    ... )
    >>> float(ds.in_units("B_1", "nT")[0, 0, 0])
    10.0
    """
    from pypic.compute import display_unit_factor  # layered above dataset

    return self.in_si(name) / display_unit_factor(unit_str)

isel(indexers=None, **kwargs)

Integer-index selection, returning a new FieldDataset.

Accepts indexers as a dict and/or **kwargs.

Parameters:

Name	Type	Description	Default
indexers	dict[str, Any] | None	Dimension-name → integer index or slice mapping.	None
**kwargs	Any	Additional dimension selections.	{}

Returns:

Type	Description
FieldDataset

Examples:

>>> import numpy as np
>>> from pypic.units import Normalization
>>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0))
>>> ds = FieldDataset.from_arrays(
...     {"B_1": np.arange(24.0).reshape(4, 3, 2)},
...     grid, Normalization.identity(),
... )
>>> face = ds.isel(x=0)
>>> face.grid.dimensions, face.grid.surviving_axis_names
((3, 2), ('y', 'z'))

Source code in src/pypic/dataset.py

def isel(
    self,
    indexers: dict[str, Any] | None = None,
    **kwargs: Any,  # noqa: ANN401 — xarray passthrough
) -> FieldDataset:
    """Integer-index selection, returning a new FieldDataset.

    Accepts ``indexers`` as a dict and/or ``**kwargs``.

    Parameters
    ----------
    indexers : dict[str, Any] | None
        Dimension-name → integer index or slice mapping.
    **kwargs
        Additional dimension selections.

    Returns
    -------
    FieldDataset

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0))
    >>> ds = FieldDataset.from_arrays(
    ...     {"B_1": np.arange(24.0).reshape(4, 3, 2)},
    ...     grid, Normalization.identity(),
    ... )
    >>> face = ds.isel(x=0)
    >>> face.grid.dimensions, face.grid.surviving_axis_names
    ((3, 2), ('y', 'z'))
    """
    merged = dict(indexers) if indexers else {}
    merged.update(kwargs)
    return self._wrap_sliced(self._ds.isel(merged))

where(cond, other=np.nan)

Mask fields where cond is False.

Returns a new FieldDataset with the same grid shape. Points where cond is False are set to other (default NaN). Useful for spatial masks (spherical cutouts, boundary regions) without reducing dimensions.

Parameters:

Name	Type	Description	Default
cond	BoolArray	Boolean array with shape matching the grid dimensions. True keeps the value, False replaces with other.	required
other	float	Fill value for masked points (default NaN).	nan

Returns:

Type	Description
FieldDataset

Examples:

>>> import numpy as np
>>> from pypic.units import Normalization
>>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0))
>>> ds = FieldDataset.from_arrays(
...     {"B_1": np.arange(24.0).reshape(4, 3, 2)},
...     grid, Normalization.identity(),
... )
>>> x, _, _ = np.meshgrid(*grid.coordinate_arrays(), indexing="ij")
>>> inner = ds.where(x < 2.0)
>>> inner["B_1"].shape, int(np.isnan(inner["B_1"]).sum())
((4, 3, 2), 12)

Source code in src/pypic/dataset.py

def where(self, cond: BoolArray, other: float = np.nan) -> FieldDataset:
    r"""Mask fields where *cond* is ``False``.

    Returns a new `FieldDataset` with the same grid shape.
    Points where *cond* is ``False`` are set to *other* (default
    ``NaN``).  Useful for spatial masks (spherical cutouts, boundary
    regions) without reducing dimensions.

    Parameters
    ----------
    cond : BoolArray
        Boolean array with shape matching the grid dimensions.
        ``True`` keeps the value, ``False`` replaces with *other*.
    other : float
        Fill value for masked points (default ``NaN``).

    Returns
    -------
    FieldDataset

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0))
    >>> ds = FieldDataset.from_arrays(
    ...     {"B_1": np.arange(24.0).reshape(4, 3, 2)},
    ...     grid, Normalization.identity(),
    ... )
    >>> x, _, _ = np.meshgrid(*grid.coordinate_arrays(), indexing="ij")
    >>> inner = ds.where(x < 2.0)
    >>> inner["B_1"].shape, int(np.isnan(inner["B_1"]).sum())
    ((4, 3, 2), 12)
    """
    axis_names = list(self._grid.surviving_axis_names)
    mask_da = xr.DataArray(cond, dims=axis_names)
    return self._wrap_sliced(self._ds.where(mask_da, other=other))

reduce(axis, **kwargs)

Reduce this dataset along one or more axes.

See pypic.reduce. Convenience method equivalent to pypic.reduce(self, axis, **kwargs). Enables fluent chaining: ds.where(mask).reduce("z", reduction="integrate").

Parameters:

Name	Type	Description	Default
axis	str or tuple of str	Surviving-axis name(s) to reduce away.	required
**kwargs	Any	Forwarded to pypic.reductions.reduce: reduction, selection, fields, nan_policy.	{}

Returns:

Type	Description
FieldDataset	With axis (or every name in the tuple) removed from the grid.

Examples:

Integrating a unit field over two cells of unit spacing gives the trapezoidal path length between the two cell centres:

>>> import numpy as np
>>> from pypic.units import Normalization
>>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0))
>>> ds = FieldDataset.from_arrays(
...     {"B_1": np.ones((4, 3, 2))}, grid, Normalization.identity(),
... )
>>> column = ds.reduce("z", reduction="integrate")
>>> column.grid.dimensions, float(column["B_1"][0, 0])
((4, 3), 1.0)

Source code in src/pypic/dataset.py

def reduce(
    self,
    axis: str | tuple[str, ...],
    **kwargs: Any,  # noqa: ANN401
) -> FieldDataset:
    r"""Reduce this dataset along one or more axes.

    See [`pypic.reduce`][pypic.reduce]. Convenience method equivalent to
    ``pypic.reduce(self, axis, **kwargs)``.  Enables fluent
    chaining: ``ds.where(mask).reduce("z", reduction="integrate")``.

    Parameters
    ----------
    axis : str or tuple of str
        Surviving-axis name(s) to reduce away.
    **kwargs
        Forwarded to [`pypic.reductions.reduce`][pypic.reductions.reduce]:
        ``reduction``, ``selection``, ``fields``, ``nan_policy``.

    Returns
    -------
    FieldDataset
        With *axis* (or every name in the tuple) removed from the grid.

    Examples
    --------
    Integrating a unit field over two cells of unit spacing gives the
    trapezoidal path length between the two cell centres:

    >>> import numpy as np
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0))
    >>> ds = FieldDataset.from_arrays(
    ...     {"B_1": np.ones((4, 3, 2))}, grid, Normalization.identity(),
    ... )
    >>> column = ds.reduce("z", reduction="integrate")
    >>> column.grid.dimensions, float(column["B_1"][0, 0])
    ((4, 3), 1.0)
    """
    from pypic.reductions import reduce as _reduce  # layered above dataset

    return _reduce(self, axis, **kwargs)

GridInfo

grid

GridInfo: the structured-grid metadata every FieldDataset carries.

GridInfo dataclass

Structured grid metadata for 1D/2D/3D simulation domains.

Parameters:

Name	Type	Description	Default
dimensions	tuple[int, ...]	Number of cells along each axis.	required
spacing	tuple[float, ...]	Cell size along each axis in code units.	required
origin	tuple[float, ...]	Lower-left corner coordinate of the domain.	()
geometry	CoordinateGeometry	Coordinate system (Cartesian, spherical, cylindrical).	CARTESIAN
dt	float | None	Timestep size in code units, if known.	None
boundary	tuple[str, ...] | None	Boundary condition per axis (e.g. ("periodic", "open", "periodic")).	None

Examples:

>>> grid = GridInfo(
...     dimensions=(4,), spacing=(0.5,), origin=(0.0,),
...     geometry=CARTESIAN,
... )
>>> grid.coordinate_arrays()[0]
array([0.25, 0.75, 1.25, 1.75])

Source code in src/pypic/grid.py

@dataclass(frozen=True, slots=True)
class GridInfo:
    r"""Structured grid metadata for 1D/2D/3D simulation domains.

    Parameters
    ----------
    dimensions : tuple[int, ...]
        Number of cells along each axis.
    spacing : tuple[float, ...]
        Cell size along each axis in code units.
    origin : tuple[float, ...]
        Lower-left corner coordinate of the domain.
    geometry : CoordinateGeometry
        Coordinate system (Cartesian, spherical, cylindrical).
    dt : float | None
        Timestep size in code units, if known.
    boundary : tuple[str, ...] | None
        Boundary condition per axis (e.g. ``("periodic", "open", "periodic")``).

    Examples
    --------
    >>> grid = GridInfo(
    ...     dimensions=(4,), spacing=(0.5,), origin=(0.0,),
    ...     geometry=CARTESIAN,
    ... )
    >>> grid.coordinate_arrays()[0]
    array([0.25, 0.75, 1.25, 1.75])
    """

    dimensions: tuple[int, ...]
    spacing: tuple[float, ...]
    origin: tuple[float, ...] = ()
    geometry: CoordinateGeometry = CARTESIAN
    dt: float | None = None
    boundary: tuple[str, ...] | None = None
    surviving_axes: tuple[int, ...] | None = None

    def __post_init__(self) -> None:
        ndim = len(self.dimensions)
        if not self.origin:
            object.__setattr__(self, "origin", (0.0,) * ndim)
        if len(self.spacing) != ndim or len(self.origin) != ndim:
            msg = (
                f"Length mismatch: dimensions({ndim}), "
                f"spacing({len(self.spacing)}), origin({len(self.origin)})"
            )
            raise ValueError(msg)
        # Without this, surviving_axis_names truncates to the geometry's
        # three axes and the mismatch surfaces two calls later, inside
        # a zip() in FieldDataset.from_arrays.
        max_ndim = len(self.geometry.axis_names)
        if ndim > max_ndim:
            msg = (
                f"GridInfo holds at most {max_ndim} dimensions, got {ndim}. "
                f"pypic's containers are 1D/2D/3D structured grids; a "
                f"higher-dimensional phase space (gyrokinetic 5D, continuum "
                f"Vlasov 6D) is described by the [phase_space] section, "
                f"which reaches SimulationConfig.phase_space as typed "
                f"metadata, but no container holds its distribution function."
            )
            raise ValueError(msg)
        for i, d in enumerate(self.dimensions):
            if d <= 0:
                raise ValueError(f"dimensions[{i}] must be > 0, got {d}")
        for i, s in enumerate(self.spacing):
            if s <= 0:
                raise ValueError(f"spacing[{i}] must be > 0, got {s}")
        if self.dt is not None and self.dt <= 0:
            raise ValueError(f"dt must be > 0, got {self.dt}")
        if self.boundary is not None and len(self.boundary) != ndim:
            msg = (
                f"boundary length ({len(self.boundary)}) "
                f"must match dimensions length ({ndim})"
            )
            raise ValueError(msg)
        if self.surviving_axes is not None and len(self.surviving_axes) != ndim:
            msg = (
                f"surviving_axes length ({len(self.surviving_axes)}) "
                f"must match dimensions length ({ndim})"
            )
            raise ValueError(msg)

    @property
    def surviving_axis_names(self) -> tuple[str, ...]:
        """Axis names for the current dimensions.

        After slicing, returns only the names of axes that survived
        (e.g. ``("x", "z")`` after removing the y-axis). When no
        slicing has occurred, returns the first *ndim* names from
        the geometry.

        Examples
        --------
        >>> grid = GridInfo(
        ...     dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0),
        ...     geometry=CARTESIAN,
        ... )
        >>> grid.surviving_axis_names
        ('x', 'y', 'z')
        >>> import copy
        >>> sliced = copy.replace(
        ...     grid, dimensions=(4, 2), spacing=(1.0, 1.0),
        ...     origin=(0.0, 0.0), surviving_axes=(0, 2),
        ... )
        >>> sliced.surviving_axis_names
        ('x', 'z')
        """
        if self.surviving_axes is not None:
            return tuple(self.geometry.axis_names[i] for i in self.surviving_axes)
        return self.geometry.axis_names[: len(self.dimensions)]

    def coordinate_arrays(self) -> tuple[FloatArray, ...]:
        r"""Cell-centered coordinate arrays for each axis.

        Returns
        -------
        tuple[FloatArray, ...]
            One 1-D array per axis: ``origin[i] + (arange(n) + 0.5) * dx[i]``.

        Examples
        --------
        >>> grid = GridInfo(
        ...     dimensions=(3, 2), spacing=(1.0, 2.0), origin=(0.0, 0.0),
        ...     geometry=CARTESIAN,
        ... )
        >>> x, y = grid.coordinate_arrays()
        >>> x
        array([0.5, 1.5, 2.5])
        >>> y
        array([1., 3.])
        """
        return tuple(
            self.origin[i] + (np.arange(self.dimensions[i]) + 0.5) * self.spacing[i]
            for i in range(len(self.dimensions))
        )

surviving_axis_names property

Axis names for the current dimensions.

After slicing, returns only the names of axes that survived (e.g. ("x", "z") after removing the y-axis). When no slicing has occurred, returns the first ndim names from the geometry.

Examples:

>>> grid = GridInfo(
...     dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0),
...     geometry=CARTESIAN,
... )
>>> grid.surviving_axis_names
('x', 'y', 'z')
>>> import copy
>>> sliced = copy.replace(
...     grid, dimensions=(4, 2), spacing=(1.0, 1.0),
...     origin=(0.0, 0.0), surviving_axes=(0, 2),
... )
>>> sliced.surviving_axis_names
('x', 'z')

coordinate_arrays()

Cell-centered coordinate arrays for each axis.

Returns:

Type	Description
tuple[FloatArray, ...]	One 1-D array per axis: origin[i] + (arange(n) + 0.5) * dx[i].

Examples:

>>> grid = GridInfo(
...     dimensions=(3, 2), spacing=(1.0, 2.0), origin=(0.0, 0.0),
...     geometry=CARTESIAN,
... )
>>> x, y = grid.coordinate_arrays()
>>> x
array([0.5, 1.5, 2.5])
>>> y
array([1., 3.])

Source code in src/pypic/grid.py

def coordinate_arrays(self) -> tuple[FloatArray, ...]:
    r"""Cell-centered coordinate arrays for each axis.

    Returns
    -------
    tuple[FloatArray, ...]
        One 1-D array per axis: ``origin[i] + (arange(n) + 0.5) * dx[i]``.

    Examples
    --------
    >>> grid = GridInfo(
    ...     dimensions=(3, 2), spacing=(1.0, 2.0), origin=(0.0, 0.0),
    ...     geometry=CARTESIAN,
    ... )
    >>> x, y = grid.coordinate_arrays()
    >>> x
    array([0.5, 1.5, 2.5])
    >>> y
    array([1., 3.])
    """
    return tuple(
        self.origin[i] + (np.arange(self.dimensions[i]) + 0.5) * self.spacing[i]
        for i in range(len(self.dimensions))
    )

Simulation, particle, and tabular containers

containers

Immutable data containers: SimulationConfig, TabularData, ParticleData.

StaggerInfo dataclass

Provenance record of the original grid stagger convention.

Readers destagger to co-located grids on load. StaggerInfo documents what the grid looked like before destaggering — purely informational, never used in computation or operators.

Parameters:

Name	Type	Description	Default
convention	str	Overall grid type: "node" (all fields on vertices), "cell" (all fields at cell centers), or "staggered" (Yee mesh — B on faces, E on edges, etc.).	required
field_locations	Mapping[str, str] | None	Per-field-group stagger locations, e.g. {"B": "face", "E": "edge"}. Only meaningful for the "staggered" convention; None otherwise. Frozen to MappingProxyType after construction.	None
position	Mapping[str, tuple[float, ...]] | None	Per-component stagger offsets in [0.0, 1.0), one tuple per canonical field component (e.g. {"B_1": (0.5, 0.0, 0.0), "E_1": (0.0, 0.5, 0.5)}). Adopts the openPMD ED-PIC position semantics so Yee-mesh PIC, BATSRUS face-centered B, and any future co-located write-out can describe their native stagger losslessly even after the reader has destaggered. Only meaningful for the "staggered" convention; None otherwise. Each tuple is converted to a tuple[float, ...] on construction so the whole structure is hashable.	None
interpolation_order	int | None	Order of interpolation used during destaggering (1 = linear, 2 = quadratic). None if no destaggering was performed.	None
notes	str | None	Free-text provenance (e.g. "Yee mesh, B on faces").	None

Examples:

>>> si = StaggerInfo(convention="node")
>>> si.convention
'node'
>>> si.field_locations is None
True

>>> si = StaggerInfo(
...     convention="staggered",
...     field_locations={"B": "face", "E": "edge"},
...     notes="Yee mesh",
... )
>>> si.field_locations["B"]
'face'

>>> si = StaggerInfo(
...     convention="staggered",
...     position={"B_1": (0.5, 0.0, 0.0), "E_1": (0.0, 0.5, 0.5)},
... )
>>> si.position["B_1"]
(0.5, 0.0, 0.0)

Source code in src/pypic/containers.py

@dataclass(frozen=True, slots=True)
class StaggerInfo:
    r"""Provenance record of the original grid stagger convention.

    Readers destagger to co-located grids on load.  ``StaggerInfo``
    documents what the grid looked like *before* destaggering — purely
    informational, never used in computation or operators.

    Parameters
    ----------
    convention : str
        Overall grid type: ``"node"`` (all fields on vertices),
        ``"cell"`` (all fields at cell centers), or ``"staggered"``
        (Yee mesh — B on faces, E on edges, etc.).
    field_locations : Mapping[str, str] | None
        Per-field-group stagger locations, e.g.
        ``{"B": "face", "E": "edge"}``.  Only meaningful for the
        ``"staggered"`` convention; ``None`` otherwise.  Frozen to
        ``MappingProxyType`` after construction.
    position : Mapping[str, tuple[float, ...]] | None
        Per-component stagger offsets in ``[0.0, 1.0)``, one tuple per
        canonical field component (e.g. ``{"B_1": (0.5, 0.0, 0.0),
        "E_1": (0.0, 0.5, 0.5)}``).  Adopts the openPMD ED-PIC
        ``position`` semantics so Yee-mesh PIC, BATSRUS face-centered B,
        and any future co-located write-out can describe their native
        stagger losslessly even after the reader has destaggered.
        Only meaningful for the ``"staggered"`` convention; ``None``
        otherwise.  Each tuple is converted to a ``tuple[float, ...]``
        on construction so the whole structure is hashable.
    interpolation_order : int | None
        Order of interpolation used during destaggering (1 = linear,
        2 = quadratic).  ``None`` if no destaggering was performed.
    notes : str | None
        Free-text provenance (e.g. ``"Yee mesh, B on faces"``).

    Examples
    --------
    >>> si = StaggerInfo(convention="node")
    >>> si.convention
    'node'
    >>> si.field_locations is None
    True

    >>> si = StaggerInfo(
    ...     convention="staggered",
    ...     field_locations={"B": "face", "E": "edge"},
    ...     notes="Yee mesh",
    ... )
    >>> si.field_locations["B"]
    'face'

    >>> si = StaggerInfo(
    ...     convention="staggered",
    ...     position={"B_1": (0.5, 0.0, 0.0), "E_1": (0.0, 0.5, 0.5)},
    ... )
    >>> si.position["B_1"]
    (0.5, 0.0, 0.0)
    """

    convention: str
    field_locations: Mapping[str, str] | None = None
    position: Mapping[str, tuple[float, ...]] | None = None
    interpolation_order: int | None = None
    notes: str | None = None

    def __post_init__(self) -> None:
        if self.field_locations is not None:
            object.__setattr__(
                self,
                "field_locations",
                MappingProxyType(dict(self.field_locations)),
            )
        if self.position is not None:
            normalized: dict[str, tuple[float, ...]] = {}
            for name, offsets in self.position.items():
                offset_tuple = tuple(float(x) for x in offsets)
                if not all(0.0 <= x < 1.0 for x in offset_tuple):
                    raise ValueError(
                        f"StaggerInfo.position[{name!r}] = {offset_tuple} — "
                        f"each offset must be in [0.0, 1.0)"
                    )
                normalized[name] = offset_tuple
            object.__setattr__(self, "position", MappingProxyType(normalized))

SimulationConfig dataclass

Parsed simulation configuration from a TOML config file.

Parameters:

Name	Type	Description	Default
model_name	str	Human-readable name for the simulation run.	required
model_type	Literal['PIC', 'MHD', 'hybrid', 'vlasov', 'gyrokinetic']	Simulation type identifier — uppercase for fluid/PIC families, lowercase for kinetic continuum codes. Matches [model].type in the schema and the Pydantic Model.type Literal. The vlasov / gyrokinetic values are accepted ahead of the readers that consume them, so a conforming document from such a code validates today.	required
grid	GridInfo	Grid metadata (includes coordinate geometry).	required
normalization	Normalization	Unit system.	required
species	tuple[SpeciesInfo, ...]	Species definitions (tuple for immutability).	()
physics	PhysicsParams	Physics parameters (frozen dataclass).	PhysicsParams()
frame	str	Reference frame label (e.g. "GSM", "simulation").	'simulation'
transforms	Mapping[str, FrameTransform]	Validated [coordinates.transforms] entries, keyed by target-frame name. Chains resolve breadth-first from frame.	dict()
initial_conditions	InitialConditions | None	Validated [initial_conditions] object from the v2.0 schema, or None when the section is absent.	None
output	Output | None	Validated [output] umbrella object from the v2.0 schema, or None when the section is absent.	None
bodies	tuple[Body, ...]	Validated [[bodies]] registry — planets, stars, coils. Drivers and initial conditions reference these by name.	()
drivers	tuple[Driver, ...]	Validated [[drivers]] registry — magnetograms, solar-wind inflows, coupled models, and other ongoing external input.	()
restart	Restart | None	Validated [restart] continuation pointer, or None.	None
run	Run | None	Validated [run] provenance object — identifier, UTC epoch, published references, authors, DOI, license, funding, embargo, ensemble, resource accounting. None only when this config was assembled by hand outside the schema path.	None
probes	tuple[Probe, ...]	Validated [[probes]] entries; empty when no probes declared.	()
collisions	tuple[Collision, ...]	Validated [[collisions]] entries; empty when no collision pairs declared.	()
phase_space	PhaseSpace | None	Validated [phase_space] block for >3D kinetic codes. Continuum-Vlasov sparse-block storage knobs live under phase_space.storage.	None
metadata	Mapping[str, Any]	Free-form annotations from readers (stagger, scaling, version, description, ...). Schema-typed sections live on dedicated attributes above rather than as opaque dict entries here.	dict()

Examples:

>>> from pypic.units import Normalization, SpeciesInfo
>>> cfg = SimulationConfig(
...     model_name="test", model_type="PIC",
...     grid=GridInfo(
...         dimensions=(4,), spacing=(1.0,), origin=(0.0,),
...         geometry=CARTESIAN,
...     ),
...     normalization=Normalization.identity(),
...     species=(SpeciesInfo(name="e", charge=-1.0, mass=1.0),),
... )
>>> cfg.model_name
'test'

Source code in src/pypic/containers.py

@dataclass(frozen=True, slots=True)
class SimulationConfig:
    """Parsed simulation configuration from a TOML config file.

    Parameters
    ----------
    model_name : str
        Human-readable name for the simulation run.
    model_type : Literal["PIC", "MHD", "hybrid", "vlasov", "gyrokinetic"]
        Simulation type identifier — uppercase for fluid/PIC families,
        lowercase for kinetic continuum codes. Matches ``[model].type``
        in the schema and the Pydantic ``Model.type`` Literal. The
        ``vlasov`` / ``gyrokinetic`` values are accepted ahead of the
        readers that consume them, so a conforming document from such a
        code validates today.
    grid : GridInfo
        Grid metadata (includes coordinate geometry).
    normalization : Normalization
        Unit system.
    species : tuple[SpeciesInfo, ...]
        Species definitions (tuple for immutability).
    physics : PhysicsParams
        Physics parameters (frozen dataclass).
    frame : str
        Reference frame label (e.g. ``"GSM"``, ``"simulation"``).
    transforms : Mapping[str, FrameTransform]
        Validated ``[coordinates.transforms]`` entries, keyed by
        target-frame name. Chains resolve breadth-first from *frame*.
    initial_conditions : InitialConditions | None
        Validated ``[initial_conditions]`` object from the v2.0 schema,
        or ``None`` when the section is absent.
    output : Output | None
        Validated ``[output]`` umbrella object from the v2.0 schema, or
        ``None`` when the section is absent.
    bodies : tuple[Body, ...]
        Validated ``[[bodies]]`` registry — planets, stars, coils.
        Drivers and initial conditions reference these by name.
    drivers : tuple[Driver, ...]
        Validated ``[[drivers]]`` registry — magnetograms, solar-wind
        inflows, coupled models, and other ongoing external input.
    restart : Restart | None
        Validated ``[restart]`` continuation pointer, or ``None``.
    run : Run | None
        Validated ``[run]`` provenance object — identifier, UTC epoch,
        published references, authors, DOI, license, funding, embargo,
        ensemble, resource accounting. ``None`` only when this config
        was assembled by hand outside the schema path.
    probes : tuple[Probe, ...]
        Validated ``[[probes]]`` entries; empty when no probes declared.
    collisions : tuple[Collision, ...]
        Validated ``[[collisions]]`` entries; empty when no collision
        pairs declared.
    phase_space : PhaseSpace | None
        Validated ``[phase_space]`` block for >3D kinetic codes.
        Continuum-Vlasov sparse-block storage knobs live under
        ``phase_space.storage``.
    metadata : Mapping[str, Any]
        Free-form annotations from readers (stagger, scaling, version,
        description, ...). Schema-typed sections live on dedicated
        attributes above rather than as opaque dict entries here.

    Examples
    --------
    >>> from pypic.units import Normalization, SpeciesInfo
    >>> cfg = SimulationConfig(
    ...     model_name="test", model_type="PIC",
    ...     grid=GridInfo(
    ...         dimensions=(4,), spacing=(1.0,), origin=(0.0,),
    ...         geometry=CARTESIAN,
    ...     ),
    ...     normalization=Normalization.identity(),
    ...     species=(SpeciesInfo(name="e", charge=-1.0, mass=1.0),),
    ... )
    >>> cfg.model_name
    'test'
    """

    model_name: str
    model_type: ModelType
    grid: GridInfo
    normalization: Normalization
    species: tuple[SpeciesInfo, ...] = ()
    physics: PhysicsParams = field(default_factory=PhysicsParams)
    frame: str = "simulation"
    transforms: Mapping[str, FrameTransform] = field(default_factory=dict)
    initial_conditions: InitialConditions | None = None
    output: Output | None = None
    bodies: tuple[Body, ...] = ()
    drivers: tuple[Driver, ...] = ()
    restart: Restart | None = None
    run: Run | None = None
    probes: tuple[Probe, ...] = ()
    collisions: tuple[Collision, ...] = ()
    phase_space: PhaseSpace | None = None
    metadata: Mapping[str, Any] = field(default_factory=dict)  # frozen in __post_init__

    def __post_init__(self) -> None:
        # Wrap mutable dicts in read-only proxies to enforce true immutability.
        # Callers pass plain dicts; frozen assignment uses object.__setattr__.
        object.__setattr__(
            self,
            "transforms",
            MappingProxyType(dict(self.transforms)),
        )
        object.__setattr__(self, "metadata", MappingProxyType(dict(self.metadata)))

TabularData dataclass

Generic columnar container for auxiliary time-series data.

Stores named 1-D arrays sharing a common length, with optional index column designation. Used for conserved quantities, solver diagnostics, virtual satellite probes, etc.

Parameters:

Name	Type	Description	Default
name	str	Dataset label (e.g. "conserved_quantities").	required
columns	Mapping[str, FloatArray]	Column name → 1-D array mapping. All arrays must have the same length.	required
index_column	str | None	Which column serves as the index (e.g. "cycle"). None means row-indexed.	None
metadata	Mapping[str, Any]	Source info (reader name, file path, etc.).	dict()

Examples:

>>> import numpy as np
>>> tab = TabularData(
...     name="diagnostics",
...     columns={"cycle": np.array([0.0, 1.0, 2.0]),
...              "energy": np.array([1.0, 0.9, 0.8])},
...     index_column="cycle",
... )
>>> tab["energy"]
array([1. , 0.9, 0.8])
>>> len(tab)
3
>>> "cycle" in tab
True
>>> tab.column_names
['cycle', 'energy']

Source code in src/pypic/containers.py

@dataclass(frozen=True, slots=True)
class TabularData:
    r"""Generic columnar container for auxiliary time-series data.

    Stores named 1-D arrays sharing a common length, with optional
    index column designation.  Used for conserved quantities, solver
    diagnostics, virtual satellite probes, etc.

    Parameters
    ----------
    name : str
        Dataset label (e.g. ``"conserved_quantities"``).
    columns : Mapping[str, FloatArray]
        Column name → 1-D array mapping.  All arrays must have
        the same length.
    index_column : str | None
        Which column serves as the index (e.g. ``"cycle"``).
        ``None`` means row-indexed.
    metadata : Mapping[str, Any]
        Source info (reader name, file path, etc.).

    Examples
    --------
    >>> import numpy as np
    >>> tab = TabularData(
    ...     name="diagnostics",
    ...     columns={"cycle": np.array([0.0, 1.0, 2.0]),
    ...              "energy": np.array([1.0, 0.9, 0.8])},
    ...     index_column="cycle",
    ... )
    >>> tab["energy"]
    array([1. , 0.9, 0.8])
    >>> len(tab)
    3
    >>> "cycle" in tab
    True
    >>> tab.column_names
    ['cycle', 'energy']
    """

    name: str
    columns: Mapping[str, FloatArray]  # frozen at runtime via __post_init__
    index_column: str | None = None
    metadata: Mapping[str, Any] = field(default_factory=dict)  # frozen at runtime

    def __post_init__(self) -> None:
        # Validate before freezing
        if self.columns:
            lengths = {k: len(v) for k, v in self.columns.items()}
            unique_lengths = set(lengths.values())
            if len(unique_lengths) > 1:
                msg = f"All columns must have equal length, got {lengths}"
                raise ValueError(msg)
        if self.index_column is not None and self.index_column not in self.columns:
            msg = (
                f"index_column {self.index_column!r} not found "
                f"in columns: {sorted(self.columns)}"
            )
            raise ValueError(msg)
        # Wrap mutable dicts in read-only proxies
        object.__setattr__(self, "columns", MappingProxyType(dict(self.columns)))
        object.__setattr__(self, "metadata", MappingProxyType(dict(self.metadata)))

    def __getitem__(self, key: str) -> FloatArray:
        """Return a column by name.

        Parameters
        ----------
        key : str
            Column name.

        Returns
        -------
        FloatArray

        Raises
        ------
        KeyError
            If *key* is not a column name.
        """
        try:
            return self.columns[key]
        except KeyError:
            msg = f"Column {key!r} not found. Available: {sorted(self.columns)}"
            raise KeyError(msg) from None

    def __contains__(self, key: object) -> bool:
        """Check whether *key* is a column name."""
        return key in self.columns

    def __len__(self) -> int:
        """Return the number of rows (common array length)."""
        if not self.columns:
            return 0
        return len(next(iter(self.columns.values())))

    @property
    def column_names(self) -> list[str]:
        """Sorted list of column names."""
        return sorted(self.columns)

    @property
    def index(self) -> FloatArray:
        """Index array: the designated index column, or ``np.arange(len)``."""
        if self.index_column is not None:
            return self.columns[self.index_column]
        return np.arange(len(self), dtype=np.float64)

column_names property

Sorted list of column names.

index property

Index array: the designated index column, or np.arange(len).

__getitem__(key)

Return a column by name.

Parameters:

Name	Type	Description	Default
key	str	Column name.	required

Returns:

Type	Description
FloatArray

Raises:

Type	Description
KeyError	If key is not a column name.

Source code in src/pypic/containers.py

def __getitem__(self, key: str) -> FloatArray:
    """Return a column by name.

    Parameters
    ----------
    key : str
        Column name.

    Returns
    -------
    FloatArray

    Raises
    ------
    KeyError
        If *key* is not a column name.
    """
    try:
        return self.columns[key]
    except KeyError:
        msg = f"Column {key!r} not found. Available: {sorted(self.columns)}"
        raise KeyError(msg) from None

__contains__(key)

Check whether key is a column name.

Source code in src/pypic/containers.py

def __contains__(self, key: object) -> bool:
    """Check whether *key* is a column name."""
    return key in self.columns

__len__()

Return the number of rows (common array length).

Source code in src/pypic/containers.py

def __len__(self) -> int:
    """Return the number of rows (common array length)."""
    if not self.columns:
        return 0
    return len(next(iter(self.columns.values())))

ParticleData dataclass

Container for particle data from a single species at one timestep.

Canonical per-macroparticle representation: every macroparticle carries a weight (number of physical particles it represents); the species as a whole carries scalar species_charge and species_mass. The per-macroparticle charge and mass that enter moments and the equations of motion are derived on demand via macro_charge and macro_mass. This shape is code-agnostic — readers for combined- storage codes (iPIC3D, OSIRIS) split the native q_s × w column into weight + scalars on load; separate-storage codes (VPIC, WarpX, Smilei, PIConGPU, TRISTAN-MP) populate the same fields directly.

Parameters:

Name	Type	Description	Default
species_index	int	Zero-based species index.	required
species_name	str	Human-readable species name (e.g. "electrons").	required
position	FloatArray | None	Particle positions, shape (N, 3). None if not loaded.	required
velocity	FloatArray | None	Particle velocities, shape (N, 3). None if not loaded.	required
n_particles	int	Total particle count.	required
id	IntArray | None	Integer particle tracking IDs, shape (N,). None if not available or not requested.	None
weight	FloatArray | None	Per-macroparticle weight \(w\) (number of physical particles per macroparticle), shape (N,), float64. Required for macro_charge and macro_mass.	None
species_charge	float | None	Scalar species charge \(q_s\) in code units (e.g. -1.0 for electrons in iPIC3D normalization).	None
species_mass	float | None	Scalar species mass \(m_s\) in code units.	None
metadata	Mapping[str, Any]	Source info (file path, format, etc.).	required

Examples:

>>> import numpy as np
>>> pcl = ParticleData(
...     species_index=0, species_name="electrons",
...     position=np.zeros((10, 3)),
...     velocity=np.ones((10, 3)),
...     n_particles=10, metadata={},
...     weight=np.ones(10), species_charge=-1.0, species_mass=1.0,
... )
>>> pcl.x.shape
(10,)
>>> len(pcl)
10

Source code in src/pypic/containers.py

@dataclass(frozen=True, slots=True)
class ParticleData:
    r"""Container for particle data from a single species at one timestep.

    Canonical per-macroparticle representation: every macroparticle carries
    a ``weight`` (number of physical particles it represents); the species
    as a whole carries scalar ``species_charge`` and ``species_mass``. The
    per-macroparticle charge and mass that enter moments and the equations
    of motion are derived on demand via `macro_charge` and
    `macro_mass`. This shape is code-agnostic — readers for combined-
    storage codes (iPIC3D, OSIRIS) split the native ``q_s × w`` column into
    ``weight`` + scalars on load; separate-storage codes (VPIC, WarpX,
    Smilei, PIConGPU, TRISTAN-MP) populate the same fields directly.

    Parameters
    ----------
    species_index : int
        Zero-based species index.
    species_name : str
        Human-readable species name (e.g. ``"electrons"``).
    position : FloatArray | None
        Particle positions, shape ``(N, 3)``. ``None`` if not loaded.
    velocity : FloatArray | None
        Particle velocities, shape ``(N, 3)``. ``None`` if not loaded.
    n_particles : int
        Total particle count.
    id : IntArray | None
        Integer particle tracking IDs, shape ``(N,)``. ``None`` if not
        available or not requested.
    weight : FloatArray | None
        Per-macroparticle weight $w$ (number of physical particles per
        macroparticle), shape ``(N,)``, float64. Required for
        `macro_charge` and `macro_mass`.
    species_charge : float | None
        Scalar species charge $q_s$ in code units (e.g. ``-1.0`` for
        electrons in iPIC3D normalization).
    species_mass : float | None
        Scalar species mass $m_s$ in code units.
    metadata : Mapping[str, Any]
        Source info (file path, format, etc.).

    Examples
    --------
    >>> import numpy as np
    >>> pcl = ParticleData(
    ...     species_index=0, species_name="electrons",
    ...     position=np.zeros((10, 3)),
    ...     velocity=np.ones((10, 3)),
    ...     n_particles=10, metadata={},
    ...     weight=np.ones(10), species_charge=-1.0, species_mass=1.0,
    ... )
    >>> pcl.x.shape
    (10,)
    >>> len(pcl)
    10
    """

    species_index: int
    species_name: str
    position: FloatArray | None
    velocity: FloatArray | None
    n_particles: int
    metadata: Mapping[str, Any]  # frozen at runtime via __post_init__
    id: IntArray | None = None
    weight: FloatArray | None = None
    species_charge: float | None = None
    species_mass: float | None = None

    def __post_init__(self) -> None:
        if self.position is None and self.velocity is None:
            msg = "At least one of position or velocity must be provided"
            raise ValueError(msg)
        n = self.n_particles
        _check_particle_shape(self.position, "position", (n, 3))
        _check_particle_shape(self.velocity, "velocity", (n, 3))
        _check_particle_shape(self.id, "id", (n,))
        _check_particle_shape(self.weight, "weight", (n,))
        if self.weight is not None and self.weight.dtype != np.float64:
            msg = f"weight must be float64, got {self.weight.dtype}"
            raise ValueError(msg)
        object.__setattr__(self, "metadata", MappingProxyType(dict(self.metadata)))

    def _component(self, attr: str, idx: int) -> FloatArray:
        """Return one column of *attr*, raising if it was not loaded."""
        arr: FloatArray | None = getattr(self, attr)
        if arr is None:
            msg = f"{attr} was not loaded"
            raise ValueError(msg)
        return arr[:, idx]

    @property
    def x(self) -> FloatArray:
        """X positions (view into ``position[:, 0]``)."""
        return self._component("position", 0)

    @property
    def y(self) -> FloatArray:
        """Y positions (view into ``position[:, 1]``)."""
        return self._component("position", 1)

    @property
    def z(self) -> FloatArray:
        """Z positions (view into ``position[:, 2]``)."""
        return self._component("position", 2)

    @property
    def vx(self) -> FloatArray:
        """X velocities (view into ``velocity[:, 0]``)."""
        return self._component("velocity", 0)

    @property
    def vy(self) -> FloatArray:
        """Y velocities (view into ``velocity[:, 1]``)."""
        return self._component("velocity", 1)

    @property
    def vz(self) -> FloatArray:
        """Z velocities (view into ``velocity[:, 2]``)."""
        return self._component("velocity", 2)

    @property
    def macro_charge(self) -> FloatArray:
        r"""Per-macroparticle charge $q_s w$."""
        if self.species_charge is None or self.weight is None:
            msg = "Cannot compute macro_charge: need both 'species_charge' and 'weight'"
            raise ValueError(msg)
        return self.species_charge * self.weight

    @property
    def macro_mass(self) -> FloatArray:
        r"""Per-macroparticle mass $m_s w$."""
        if self.species_mass is None or self.weight is None:
            msg = "Cannot compute macro_mass: need both 'species_mass' and 'weight'"
            raise ValueError(msg)
        return self.species_mass * self.weight

    def __len__(self) -> int:
        return self.n_particles

    def __repr__(self) -> str:
        loaded = []
        if self.position is not None:
            loaded.append("position")
        if self.velocity is not None:
            loaded.append("velocity")
        if self.id is not None:
            loaded.append("id")
        if self.weight is not None:
            loaded.append("weight")
        if self.species_charge is not None:
            loaded.append(f"species_charge={self.species_charge:g}")
        if self.species_mass is not None:
            loaded.append(f"species_mass={self.species_mass:g}")
        return (
            f"ParticleData({self.species_name!r}, "
            f"n={self.n_particles:,}, "
            f"loaded=[{', '.join(loaded)}])"
        )

x property

X positions (view into position[:, 0]).

y property

Y positions (view into position[:, 1]).

z property

Z positions (view into position[:, 2]).

vx property

X velocities (view into velocity[:, 0]).

vy property

Y velocities (view into velocity[:, 1]).

vz property

Z velocities (view into velocity[:, 2]).

macro_charge property

Per-macroparticle charge \(q_s w\).

macro_mass property

Per-macroparticle mass \(m_s w\).