Coordinates & Operators¶

Coordinate geometry definitions and geometry-aware discrete differential operators (divergence, curl, gradient).

`coordinates` ¶

Coordinate geometry definitions and discrete differential operators.

`CoordinateGeometry` `dataclass` ¶

Coordinate system with axis metadata and metric scale factors.

Parameters:

Name	Type	Description	Default
`type`	`GeometryType`	The coordinate geometry type.	required
`axis_names`	`tuple[str, str, str]`	Human-readable axis labels, e.g. `("x", "y", "z")`.	required
`axis_units`	`tuple[str, str, str]`	Dimension kind per axis: `"length"` or `"angle"`.	required

Examples:

>>> CARTESIAN.axis_names
('x', 'y', 'z')
>>> SPHERICAL.axis_names
('r', 'θ', 'φ')
>>> CYLINDRICAL.axis_units
('length', 'angle', 'length')

Source code in src/pypic/coordinates/geometry.py

@dataclass(frozen=True, slots=True)
class CoordinateGeometry:
    r"""Coordinate system with axis metadata and metric scale factors.

    Parameters
    ----------
    type : GeometryType
        The coordinate geometry type.
    axis_names : tuple[str, str, str]
        Human-readable axis labels, e.g. ``("x", "y", "z")``.
    axis_units : tuple[str, str, str]
        Dimension kind per axis: ``"length"`` or ``"angle"``.

    Examples
    --------
    >>> CARTESIAN.axis_names
    ('x', 'y', 'z')
    >>> SPHERICAL.axis_names
    ('r', 'θ', 'φ')
    >>> CYLINDRICAL.axis_units
    ('length', 'angle', 'length')
    """

    type: GeometryType
    axis_names: tuple[str, str, str]
    axis_units: tuple[str, str, str]

    def metric_factors(self, x1: Numeric, x2: Numeric, x3: Numeric) -> ScaleFactors:
        r"""Compute metric scale factors $(h_1, h_2, h_3)$ for this geometry.

        The line element is $ds^2 = h_1^2 dx_1^2 + h_2^2 dx_2^2 + h_3^2 dx_3^2$.

        Parameters
        ----------
        x1 : Numeric
            First coordinate: $x$ (Cartesian), $r$ (spherical/cylindrical).
        x2 : Numeric
            Second coordinate: $y$ (Cartesian), $θ$ (spherical), $φ$ (cylindrical).
        x3 : Numeric
            Third coordinate: $z$ (Cartesian), $φ$ (spherical), $z$ (cylindrical).

        Returns
        -------
        tuple[ScaleFactor, ScaleFactor, ScaleFactor]
            Scale factors $(h_1, h_2, h_3)$. Constant factors are returned as
            ``np.float64(1.0)`` which broadcasts with any array shape.

        Examples
        --------
        >>> CARTESIAN.metric_factors(1.0, 2.0, 3.0)
        (np.float64(1.0), np.float64(1.0), np.float64(1.0))

        >>> h1, h2, h3 = SPHERICAL.metric_factors(2.0, np.pi / 2, 0.0)
        >>> float(h1), float(h2), float(h3)
        (1.0, 2.0, 2.0)

        >>> h1, h2, h3 = CYLINDRICAL.metric_factors(3.0, 0.0, 1.0)
        >>> float(h1), float(h2), float(h3)
        (1.0, 3.0, 1.0)
        """
        _one = np.float64(1.0)
        match self.type:
            case GeometryType.CARTESIAN:
                return (_one, _one, _one)
            case GeometryType.SPHERICAL:
                r = np.asarray(x1, dtype=np.float64)
                theta = np.asarray(x2, dtype=np.float64)
                return (_one, r, r * np.sin(theta))  # ds² = dr² + r²dθ² + r²sin²θ dφ²
            case GeometryType.CYLINDRICAL:
                r = np.asarray(x1, dtype=np.float64)
                return (_one, r, _one)  # ds² = dr² + r²dφ² + dz²
            case _ as unreachable:
                assert_never(unreachable)

`metric_factors(x1, x2, x3)` ¶

Compute metric scale factors \((h_1, h_2, h_3)\) for this geometry.

The line element is \(ds^2 = h_1^2 dx_1^2 + h_2^2 dx_2^2 + h_3^2 dx_3^2\).

Parameters:

Name	Type	Description	Default
`x1`	`Numeric`	First coordinate: \(x\) (Cartesian), \(r\) (spherical/cylindrical).	required
`x2`	`Numeric`	Second coordinate: \(y\) (Cartesian), \(θ\) (spherical), \(φ\) (cylindrical).	required
`x3`	`Numeric`	Third coordinate: \(z\) (Cartesian), \(φ\) (spherical), \(z\) (cylindrical).	required

Returns:

Type	Description
`tuple[ScaleFactor, ScaleFactor, ScaleFactor]`	Scale factors \((h_1, h_2, h_3)\). Constant factors are returned as `np.float64(1.0)` which broadcasts with any array shape.

Examples:

>>> CARTESIAN.metric_factors(1.0, 2.0, 3.0)
(np.float64(1.0), np.float64(1.0), np.float64(1.0))

>>> h1, h2, h3 = SPHERICAL.metric_factors(2.0, np.pi / 2, 0.0)
>>> float(h1), float(h2), float(h3)
(1.0, 2.0, 2.0)

>>> h1, h2, h3 = CYLINDRICAL.metric_factors(3.0, 0.0, 1.0)
>>> float(h1), float(h2), float(h3)
(1.0, 3.0, 1.0)

Source code in src/pypic/coordinates/geometry.py

def metric_factors(self, x1: Numeric, x2: Numeric, x3: Numeric) -> ScaleFactors:
    r"""Compute metric scale factors $(h_1, h_2, h_3)$ for this geometry.

    The line element is $ds^2 = h_1^2 dx_1^2 + h_2^2 dx_2^2 + h_3^2 dx_3^2$.

    Parameters
    ----------
    x1 : Numeric
        First coordinate: $x$ (Cartesian), $r$ (spherical/cylindrical).
    x2 : Numeric
        Second coordinate: $y$ (Cartesian), $θ$ (spherical), $φ$ (cylindrical).
    x3 : Numeric
        Third coordinate: $z$ (Cartesian), $φ$ (spherical), $z$ (cylindrical).

    Returns
    -------
    tuple[ScaleFactor, ScaleFactor, ScaleFactor]
        Scale factors $(h_1, h_2, h_3)$. Constant factors are returned as
        ``np.float64(1.0)`` which broadcasts with any array shape.

    Examples
    --------
    >>> CARTESIAN.metric_factors(1.0, 2.0, 3.0)
    (np.float64(1.0), np.float64(1.0), np.float64(1.0))

    >>> h1, h2, h3 = SPHERICAL.metric_factors(2.0, np.pi / 2, 0.0)
    >>> float(h1), float(h2), float(h3)
    (1.0, 2.0, 2.0)

    >>> h1, h2, h3 = CYLINDRICAL.metric_factors(3.0, 0.0, 1.0)
    >>> float(h1), float(h2), float(h3)
    (1.0, 3.0, 1.0)
    """
    _one = np.float64(1.0)
    match self.type:
        case GeometryType.CARTESIAN:
            return (_one, _one, _one)
        case GeometryType.SPHERICAL:
            r = np.asarray(x1, dtype=np.float64)
            theta = np.asarray(x2, dtype=np.float64)
            return (_one, r, r * np.sin(theta))  # ds² = dr² + r²dθ² + r²sin²θ dφ²
        case GeometryType.CYLINDRICAL:
            r = np.asarray(x1, dtype=np.float64)
            return (_one, r, _one)  # ds² = dr² + r²dφ² + dz²
        case _ as unreachable:
            assert_never(unreachable)

`GeometryType` ¶

Bases: StrEnum

Coordinate geometry type.

Values match the [coordinates] geometry key in schema.md.

Examples:

>>> GeometryType.CARTESIAN
<GeometryType.CARTESIAN: 'cartesian'>
>>> GeometryType("spherical")
<GeometryType.SPHERICAL: 'spherical'>

Source code in src/pypic/coordinates/geometry.py

class GeometryType(StrEnum):
    r"""Coordinate geometry type.

    Values match the ``[coordinates] geometry`` key in schema.md.

    Examples
    --------
    >>> GeometryType.CARTESIAN
    <GeometryType.CARTESIAN: 'cartesian'>
    >>> GeometryType("spherical")
    <GeometryType.SPHERICAL: 'spherical'>
    """

    CARTESIAN = "cartesian"
    SPHERICAL = "spherical"
    CYLINDRICAL = "cylindrical"

`FrameTransform` `dataclass` ¶

Affine transformation between coordinate reference frames.

Transforms a point \(\mathbf{x}\) from the source frame to the target:

\[\mathbf{x}_{target} = s \cdot R \cdot (\mathbf{x}_{source} - \mathbf{o})\]

Parameters:

Name	Type	Description	Default
`source_frame`	`str`	Name of the source frame (e.g. `"simulation"`).	required
`target_frame`	`str`	Name of the target frame (e.g. `"GSM"`).	required
`origin`	`tuple[float, float, float]`	Source-frame coordinates of the target-frame origin. Subtracted before rotation.	`(0.0, 0.0, 0.0)`
`rotation`	`Rotation3x3`	3×3 orthogonal rotation matrix as nested tuples.	`_IDENTITY_3X3`
`scale`	`float`	Converts code length units to target-frame units. `scale=0.25` with target in R_E means 1 d_i = 0.25 R_E. Default 1.0 (code and target use the same length unit). Auto-computed from `physical_extent` when available.	`1.0`
`target_axis_names`	`tuple[str, str, str] \| None`	Axis names in the target frame. `None` keeps source names.	`None`

Examples:

>>> t = FrameTransform("sim", "GSM", origin=(52.0, 26.0, 64.0))
>>> t.source_frame
'sim'
>>> t.is_identity
False
>>> FrameTransform("a", "a").is_identity
True

Source code in src/pypic/coordinates/transforms.py

@dataclass(frozen=True, slots=True)
class FrameTransform:
    r"""Affine transformation between coordinate reference frames.

    Transforms a point $\mathbf{x}$ from the source frame to the target:

    $$\mathbf{x}_{target} = s \cdot R \cdot (\mathbf{x}_{source} - \mathbf{o})$$

    Parameters
    ----------
    source_frame : str
        Name of the source frame (e.g. ``"simulation"``).
    target_frame : str
        Name of the target frame (e.g. ``"GSM"``).
    origin : tuple[float, float, float]
        Source-frame coordinates of the target-frame origin.
        Subtracted before rotation.
    rotation : Rotation3x3
        3×3 orthogonal rotation matrix as nested tuples.
    scale : float
        Converts code length units to target-frame units.
        ``scale=0.25`` with target in R_E means 1 d_i = 0.25 R_E.
        Default 1.0 (code and target use the same length unit).
        Auto-computed from ``physical_extent`` when available.
    target_axis_names : tuple[str, str, str] | None
        Axis names in the target frame. ``None`` keeps source names.

    Examples
    --------
    >>> t = FrameTransform("sim", "GSM", origin=(52.0, 26.0, 64.0))
    >>> t.source_frame
    'sim'
    >>> t.is_identity
    False
    >>> FrameTransform("a", "a").is_identity
    True
    """

    source_frame: str
    target_frame: str
    origin: tuple[float, float, float] = (0.0, 0.0, 0.0)
    rotation: Rotation3x3 = _IDENTITY_3X3
    scale: float = 1.0
    target_axis_names: tuple[str, str, str] | None = None

    def __post_init__(self) -> None:
        if len(self.rotation) != 3 or any(len(row) != 3 for row in self.rotation):
            msg = f"rotation must be 3×3, got shape {len(self.rotation)}×..."
            raise ValueError(msg)
        if self.scale <= 0:
            raise ValueError(f"scale must be > 0, got {self.scale}")
        r = self.rotation_matrix
        err = np.max(np.abs(r.T @ r - np.eye(3)))
        if err > 1e-6:
            raise ValueError(
                f"rotation matrix is not orthogonal (max |R^T R - I| = {err:.2e})"
            )

    @property
    def rotation_matrix(self) -> FloatArray:
        """Return the rotation as a (3, 3) NumPy array.

        Examples
        --------
        >>> FrameTransform("a", "b").rotation_matrix.shape
        (3, 3)
        """
        return np.array(self.rotation, dtype=np.float64)

    @property
    def is_identity(self) -> bool:
        """True if this transform is a no-op.

        Examples
        --------
        >>> FrameTransform("a", "a").is_identity
        True
        >>> FrameTransform("a", "b", origin=(1.0, 0.0, 0.0)).is_identity
        False
        """
        return (
            self.origin == (0.0, 0.0, 0.0)
            and self.rotation == _IDENTITY_3X3
            and self.scale == 1.0
        )

    def inverse(self) -> FrameTransform:
        r"""Return the inverse transform (target → source).

        For $\mathbf{x}_t = s R (\mathbf{x}_s - \mathbf{o})$, the inverse
        is $\mathbf{x}_s = R^T \mathbf{x}_t / s + \mathbf{o}$.

        Examples
        --------
        >>> t = FrameTransform("a", "b", origin=(1.0, 2.0, 3.0), scale=2.0)
        >>> inv = t.inverse()
        >>> inv.source_frame, inv.target_frame
        ('b', 'a')
        """
        r = self.rotation_matrix
        o = np.array(self.origin, dtype=np.float64)
        # origin_inv = -s R o (from inverting the affine map)
        return FrameTransform(
            source_frame=self.target_frame,
            target_frame=self.source_frame,
            origin=_to_vec3(-self.scale * r @ o),
            rotation=_to_rotation(r.T),
            scale=1.0 / self.scale,
            target_axis_names=None,
        )

`rotation_matrix` `property` ¶

Return the rotation as a (3, 3) NumPy array.

Examples:

>>> FrameTransform("a", "b").rotation_matrix.shape
(3, 3)

`is_identity` `property` ¶

True if this transform is a no-op.

Examples:

>>> FrameTransform("a", "a").is_identity
True
>>> FrameTransform("a", "b", origin=(1.0, 0.0, 0.0)).is_identity
False

`inverse()` ¶

Return the inverse transform (target → source).

For \(\mathbf{x}_t = s R (\mathbf{x}_s - \mathbf{o})\), the inverse is \(\mathbf{x}_s = R^T \mathbf{x}_t / s + \mathbf{o}\).

Examples:

>>> t = FrameTransform("a", "b", origin=(1.0, 2.0, 3.0), scale=2.0)
>>> inv = t.inverse()
>>> inv.source_frame, inv.target_frame
('b', 'a')

Source code in src/pypic/coordinates/transforms.py

def inverse(self) -> FrameTransform:
    r"""Return the inverse transform (target → source).

    For $\mathbf{x}_t = s R (\mathbf{x}_s - \mathbf{o})$, the inverse
    is $\mathbf{x}_s = R^T \mathbf{x}_t / s + \mathbf{o}$.

    Examples
    --------
    >>> t = FrameTransform("a", "b", origin=(1.0, 2.0, 3.0), scale=2.0)
    >>> inv = t.inverse()
    >>> inv.source_frame, inv.target_frame
    ('b', 'a')
    """
    r = self.rotation_matrix
    o = np.array(self.origin, dtype=np.float64)
    # origin_inv = -s R o (from inverting the affine map)
    return FrameTransform(
        source_frame=self.target_frame,
        target_frame=self.source_frame,
        origin=_to_vec3(-self.scale * r @ o),
        rotation=_to_rotation(r.T),
        scale=1.0 / self.scale,
        target_axis_names=None,
    )

`curl(f1, f2, f3, d1, d2, d3=None, *, geometry=GeometryType.CARTESIAN)` ¶

Compute the curl of a vector field.

\[(\nabla \times \mathbf{F})_1 = \frac{\partial F_3}{\partial y} - \frac{\partial F_2}{\partial z}\]

\[(\nabla \times \mathbf{F})_2 = \frac{\partial F_1}{\partial z} - \frac{\partial F_3}{\partial x}\]

\[(\nabla \times \mathbf{F})_3 = \frac{\partial F_2}{\partial x} - \frac{\partial F_1}{\partial y}\]

Parameters:

Name	Type	Description	Default
`f1`	`NDArray`	First component of the vector field, shape `(nx, ny[, nz])`.	required
`f2`	`NDArray`	Second component of the vector field, same shape as f1.	required
`f3`	`NDArray`	Third component of the vector field, same shape as f1.	required
`d1`	`float`	Grid spacing along the first axis.	required
`d2`	`float`	Grid spacing along the second axis.	required
`d3`	`float or None`	Grid spacing along the third axis, or `None` for 2D data, where \(\partial/\partial x_3 \equiv 0\) drops every term that differentiates along it.	`None`
`geometry`	`GeometryType`	Coordinate geometry. Only `CARTESIAN` is currently supported.	`CARTESIAN`

Returns:

Type	Description
`tuple[NDArray, NDArray, NDArray]`	Curl components `(curl_1, curl_2, curl_3)`, each with the same shape as the input arrays.

Raises:

Type	Description
`GeometryUnsupportedError`	If `geometry` is spherical or cylindrical.

Examples:

>>> import numpy as np
>>> f = np.ones((4, 4, 4))
>>> c1, c2, c3 = curl(f, f, f, 1.0, 1.0, 1.0)
>>> np.max(np.abs(c1))
np.float64(0.0)

Source code in src/pypic/coordinates/operators.py

def curl(
    f1: FloatArray,
    f2: FloatArray,
    f3: FloatArray,
    d1: float,
    d2: float,
    d3: float | None = None,
    *,
    geometry: GeometryType = GeometryType.CARTESIAN,
) -> tuple[FloatArray, FloatArray, FloatArray]:
    r"""Compute the curl of a vector field.

    $$(\nabla \times \mathbf{F})_1
    = \frac{\partial F_3}{\partial y} - \frac{\partial F_2}{\partial z}$$

    $$(\nabla \times \mathbf{F})_2
    = \frac{\partial F_1}{\partial z} - \frac{\partial F_3}{\partial x}$$

    $$(\nabla \times \mathbf{F})_3
    = \frac{\partial F_2}{\partial x} - \frac{\partial F_1}{\partial y}$$

    Parameters
    ----------
    f1 : NDArray
        First component of the vector field, shape ``(nx, ny[, nz])``.
    f2 : NDArray
        Second component of the vector field, same shape as *f1*.
    f3 : NDArray
        Third component of the vector field, same shape as *f1*.
    d1 : float
        Grid spacing along the first axis.
    d2 : float
        Grid spacing along the second axis.
    d3 : float or None
        Grid spacing along the third axis, or ``None`` for 2D data,
        where $\partial/\partial x_3 \equiv 0$ drops every term that
        differentiates along it.
    geometry : GeometryType
        Coordinate geometry. Only ``CARTESIAN`` is currently supported.

    Returns
    -------
    tuple[NDArray, NDArray, NDArray]
        Curl components ``(curl_1, curl_2, curl_3)``, each with the same
        shape as the input arrays.

    Raises
    ------
    GeometryUnsupportedError
        If ``geometry`` is spherical or cylindrical.

    Examples
    --------
    >>> import numpy as np
    >>> f = np.ones((4, 4, 4))
    >>> c1, c2, c3 = curl(f, f, f, 1.0, 1.0, 1.0)
    >>> np.max(np.abs(c1))
    np.float64(0.0)
    """
    _require_cartesian(geometry, "curl")
    _require_positive_spacing(d1, d2, d3)
    df3_d2: FloatArray = np.gradient(f3, d2, axis=1)
    df3_d1: FloatArray = np.gradient(f3, d1, axis=0)
    curl_3: FloatArray = np.gradient(f2, d1, axis=0) - np.gradient(f1, d2, axis=1)
    if d3 is None:
        return (df3_d2, -df3_d1, curl_3)
    curl_1: FloatArray = df3_d2 - np.gradient(f2, d3, axis=2)
    curl_2: FloatArray = np.gradient(f1, d3, axis=2) - df3_d1
    return (curl_1, curl_2, curl_3)

`divergence(f1, f2, f3, d1, d2, d3=None, *, geometry=GeometryType.CARTESIAN)` ¶

Compute the divergence of a vector field.

\[\nabla \cdot \mathbf{F} = \frac{\partial F_1}{\partial x} + \frac{\partial F_2}{\partial y} + \frac{\partial F_3}{\partial z}\]

Parameters:

Name	Type	Description	Default
`f1`	`NDArray`	First component of the vector field, shape `(nx, ny[, nz])`.	required
`f2`	`NDArray`	Second component of the vector field, same shape as f1.	required
`f3`	`NDArray`	Third component of the vector field, same shape as f1.	required
`d1`	`float`	Grid spacing along the first axis.	required
`d2`	`float`	Grid spacing along the second axis.	required
`d3`	`float or None`	Grid spacing along the third axis, or `None` for 2D data, where \(\partial/\partial x_3 \equiv 0\) drops the third term.	`None`
`geometry`	`GeometryType`	Coordinate geometry. Only `CARTESIAN` is currently supported.	`CARTESIAN`

Returns:

Type	Description
`NDArray`	Divergence field, same shape as the input arrays.

Raises:

Type	Description
`GeometryUnsupportedError`	If `geometry` is spherical or cylindrical.

Examples:

>>> import numpy as np
>>> f = np.ones((4, 4, 4))
>>> np.max(np.abs(divergence(f, f, f, 1.0, 1.0, 1.0)))
np.float64(0.0)

Source code in src/pypic/coordinates/operators.py

def divergence(
    f1: FloatArray,
    f2: FloatArray,
    f3: FloatArray,
    d1: float,
    d2: float,
    d3: float | None = None,
    *,
    geometry: GeometryType = GeometryType.CARTESIAN,
) -> FloatArray:
    r"""Compute the divergence of a vector field.

    $$\nabla \cdot \mathbf{F} = \frac{\partial F_1}{\partial x}
    + \frac{\partial F_2}{\partial y}
    + \frac{\partial F_3}{\partial z}$$

    Parameters
    ----------
    f1 : NDArray
        First component of the vector field, shape ``(nx, ny[, nz])``.
    f2 : NDArray
        Second component of the vector field, same shape as *f1*.
    f3 : NDArray
        Third component of the vector field, same shape as *f1*.
    d1 : float
        Grid spacing along the first axis.
    d2 : float
        Grid spacing along the second axis.
    d3 : float or None
        Grid spacing along the third axis, or ``None`` for 2D data,
        where $\partial/\partial x_3 \equiv 0$ drops the third term.
    geometry : GeometryType
        Coordinate geometry. Only ``CARTESIAN`` is currently supported.

    Returns
    -------
    NDArray
        Divergence field, same shape as the input arrays.

    Raises
    ------
    GeometryUnsupportedError
        If ``geometry`` is spherical or cylindrical.

    Examples
    --------
    >>> import numpy as np
    >>> f = np.ones((4, 4, 4))
    >>> np.max(np.abs(divergence(f, f, f, 1.0, 1.0, 1.0)))
    np.float64(0.0)
    """
    _require_cartesian(geometry, "divergence")
    _require_positive_spacing(d1, d2, d3)
    df1_d1: FloatArray = np.gradient(f1, d1, axis=0)
    df2_d2: FloatArray = np.gradient(f2, d2, axis=1)
    if d3 is None:
        result_2d: FloatArray = df1_d1 + df2_d2
        return result_2d
    df3_d3: FloatArray = np.gradient(f3, d3, axis=2)
    result: FloatArray = df1_d1 + df2_d2 + df3_d3
    return result

`gradient(f, d1, d2, d3=None, *, geometry=GeometryType.CARTESIAN)` ¶

Compute the gradient of a scalar field.

\[(\nabla f)_i = \frac{\partial f}{\partial x_i} \quad \text{for } i = 1, 2, 3\]

Parameters:

Name	Type	Description	Default
`f`	`NDArray`	Scalar field, shape `(nx, ny[, nz])`.	required
`d1`	`float`	Grid spacing along the first axis.	required
`d2`	`float`	Grid spacing along the second axis.	required
`d3`	`float or None`	Grid spacing along the third axis, or `None` for 2D data, where the third component is identically zero.	`None`
`geometry`	`GeometryType`	Coordinate geometry. Only `CARTESIAN` is currently supported.	`CARTESIAN`

Returns:

Type	Description
`tuple[NDArray, NDArray, NDArray]`	Gradient components `(df_d1, df_d2, df_d3)`, each with the same shape as the input array.

Raises:

Type	Description
`GeometryUnsupportedError`	If `geometry` is spherical or cylindrical.

Examples:

>>> import numpy as np
>>> f = np.ones((4, 4, 4))
>>> g1, g2, g3 = gradient(f, 1.0, 1.0, 1.0)
>>> np.max(np.abs(g1))
np.float64(0.0)

Source code in src/pypic/coordinates/operators.py

def gradient(
    f: FloatArray,
    d1: float,
    d2: float,
    d3: float | None = None,
    *,
    geometry: GeometryType = GeometryType.CARTESIAN,
) -> tuple[FloatArray, FloatArray, FloatArray]:
    r"""Compute the gradient of a scalar field.

    $$(\nabla f)_i = \frac{\partial f}{\partial x_i}
    \quad \text{for } i = 1, 2, 3$$

    Parameters
    ----------
    f : NDArray
        Scalar field, shape ``(nx, ny[, nz])``.
    d1 : float
        Grid spacing along the first axis.
    d2 : float
        Grid spacing along the second axis.
    d3 : float or None
        Grid spacing along the third axis, or ``None`` for 2D data,
        where the third component is identically zero.
    geometry : GeometryType
        Coordinate geometry. Only ``CARTESIAN`` is currently supported.

    Returns
    -------
    tuple[NDArray, NDArray, NDArray]
        Gradient components ``(df_d1, df_d2, df_d3)``, each with the same
        shape as the input array.

    Raises
    ------
    GeometryUnsupportedError
        If ``geometry`` is spherical or cylindrical.

    Examples
    --------
    >>> import numpy as np
    >>> f = np.ones((4, 4, 4))
    >>> g1, g2, g3 = gradient(f, 1.0, 1.0, 1.0)
    >>> np.max(np.abs(g1))
    np.float64(0.0)
    """
    _require_cartesian(geometry, "gradient")
    _require_positive_spacing(d1, d2, d3)
    df_d1: FloatArray = np.gradient(f, d1, axis=0)
    df_d2: FloatArray = np.gradient(f, d2, axis=1)
    df_d3: FloatArray = np.zeros_like(f) if d3 is None else np.gradient(f, d3, axis=2)
    return (df_d1, df_d2, df_d3)

`compose_transforms(first, second)` ¶

Compose two transforms: apply first then second.

If first maps A→B and second maps B→C, the result maps A→C.

The composition follows from:

\[\mathbf{x}_C = s_2 R_2 (s_1 R_1 (\mathbf{x}_A - \mathbf{o}_1) - \mathbf{o}_2)\]

Parameters:

Name	Type	Description	Default
`first`	`FrameTransform`	Transform applied first (A→B).	required
`second`	`FrameTransform`	Transform applied second (B→C).	required

Returns:

Type	Description
`FrameTransform`	Composed transform (A→C).

Raises:

Type	Description
`ValueError`	If `first.target_frame != second.source_frame`.

Examples:

>>> a_to_b = FrameTransform("A", "B", origin=(1.0, 0.0, 0.0))
>>> b_to_c = FrameTransform("B", "C", origin=(0.0, 2.0, 0.0))
>>> a_to_c = compose_transforms(a_to_b, b_to_c)
>>> a_to_c.source_frame, a_to_c.target_frame
('A', 'C')

Source code in src/pypic/coordinates/transforms.py

def compose_transforms(first: FrameTransform, second: FrameTransform) -> FrameTransform:
    r"""Compose two transforms: apply *first* then *second*.

    If *first* maps A→B and *second* maps B→C, the result maps A→C.

    The composition follows from:

    $$\mathbf{x}_C = s_2 R_2 (s_1 R_1 (\mathbf{x}_A - \mathbf{o}_1) - \mathbf{o}_2)$$

    Parameters
    ----------
    first : FrameTransform
        Transform applied first (A→B).
    second : FrameTransform
        Transform applied second (B→C).

    Returns
    -------
    FrameTransform
        Composed transform (A→C).

    Raises
    ------
    ValueError
        If ``first.target_frame != second.source_frame``.

    Examples
    --------
    >>> a_to_b = FrameTransform("A", "B", origin=(1.0, 0.0, 0.0))
    >>> b_to_c = FrameTransform("B", "C", origin=(0.0, 2.0, 0.0))
    >>> a_to_c = compose_transforms(a_to_b, b_to_c)
    >>> a_to_c.source_frame, a_to_c.target_frame
    ('A', 'C')
    """
    if first.target_frame != second.source_frame:
        msg = (
            f"Cannot compose: first maps to {first.target_frame!r} "
            f"but second starts from {second.source_frame!r}"
        )
        raise ValueError(msg)

    r1 = first.rotation_matrix
    r2 = second.rotation_matrix
    o1 = np.array(first.origin, dtype=np.float64)
    o2 = np.array(second.origin, dtype=np.float64)

    # x_C = s2 R2 (s1 R1 (x_A - o1) - o2)
    #      = s2 s1 R2 R1 (x_A - o1) - s2 R2 o2
    #      = s_c R_c (x_A - o_c)
    # where R_c = R2 R1, s_c = s1 s2
    # s_c R_c (x_A - o_c) = s_c R_c x_A - s_c R_c o_c
    # must equal: s_c R_c x_A - s_c R_c o1 - s2 R2 o2
    # => o_c = o1 + (1/s1) R1^T o2
    combined_rotation = r2 @ r1
    combined_scale = first.scale * second.scale
    combined_origin = o1 + (1.0 / first.scale) * (r1.T @ o2)

    return FrameTransform(
        source_frame=first.source_frame,
        target_frame=second.target_frame,
        origin=_to_vec3(combined_origin),
        rotation=_to_rotation(combined_rotation),
        scale=combined_scale,
        target_axis_names=second.target_axis_names or first.target_axis_names,
    )

`find_pressure_tensor_groups(field_names)` ¶

Group pressure tensor fields into complete symmetric tensors.

Returns (P_11, P_22, P_33, P_12, P_13, P_23) tuples for each complete set. Handles per-species tensors (P_s0_11 etc.).

Parameters:

Name	Type	Description	Default
`field_names`	`Iterable[str]`	All field names in the dataset.	required

Returns:

Type	Description
`list[tuple[str, str, str, str, str, str]]`

Examples:

>>> fields = ["P_11", "P_22", "P_33", "P_12", "P_13", "P_23"]
>>> find_pressure_tensor_groups(fields)
[('P_11', 'P_22', 'P_33', 'P_12', 'P_13', 'P_23')]

Source code in src/pypic/coordinates/transforms.py

def find_pressure_tensor_groups(
    field_names: Iterable[str],
) -> list[tuple[str, str, str, str, str, str]]:
    """Group pressure tensor fields into complete symmetric tensors.

    Returns ``(P_11, P_22, P_33, P_12, P_13, P_23)`` tuples for each complete
    set. Handles per-species tensors (``P_s0_11`` etc.).

    Parameters
    ----------
    field_names : Iterable[str]
        All field names in the dataset.

    Returns
    -------
    list[tuple[str, str, str, str, str, str]]

    Examples
    --------
    >>> fields = ["P_11", "P_22", "P_33", "P_12", "P_13", "P_23"]
    >>> find_pressure_tensor_groups(fields)
    [('P_11', 'P_22', 'P_33', 'P_12', 'P_13', 'P_23')]
    """
    groups: dict[str, dict[str, str]] = {}
    for name in field_names:
        m = _PRESSURE_RE.match(name)
        if m:
            species, i, j = m.groups()
            key = f"P_s{species}" if species else "P"
            groups.setdefault(key, {})[f"{i}{j}"] = name
    required = {"11", "22", "33", "12", "13", "23"}
    return [
        (g["11"], g["22"], g["33"], g["12"], g["13"], g["23"])
        for g in groups.values()
        if required <= set(g)
    ]

`find_vector_triplets(field_names)` ¶

Group field names into vector triplets needing rotation.

Returns a list of (name1, name2, name3) tuples for each complete vector field, stored or derived, that the field registry knows as a vector: total fields (B_1, B_2, B_3), per-species fields (J_s0_1, J_s0_2, J_s0_3) and derived vectors (E_prime_1, ...).

Parameters:

Name	Type	Description	Default
`field_names`	`Iterable[str]`	All field names in the dataset.	required

Returns:

Type	Description
`list[tuple[str, str, str]]`	Complete vector triplets.

Examples:

>>> find_vector_triplets(["B_1", "B_2", "B_3", "rho_c"])
[('B_1', 'B_2', 'B_3')]
>>> find_vector_triplets(["J_s0_1", "J_s0_2", "J_s0_3", "J_s1_1"])
[('J_s0_1', 'J_s0_2', 'J_s0_3')]

Source code in src/pypic/coordinates/transforms.py

def find_vector_triplets(
    field_names: Iterable[str],
) -> list[tuple[str, str, str]]:
    """Group field names into vector triplets needing rotation.

    Returns a list of ``(name1, name2, name3)`` tuples for each complete
    vector field, stored or derived, that the field registry knows as a
    vector: total fields (``B_1, B_2, B_3``), per-species fields
    (``J_s0_1, J_s0_2, J_s0_3``) and derived vectors (``E_prime_1, ...``).

    Parameters
    ----------
    field_names : Iterable[str]
        All field names in the dataset.

    Returns
    -------
    list[tuple[str, str, str]]
        Complete vector triplets.

    Examples
    --------
    >>> find_vector_triplets(["B_1", "B_2", "B_3", "rho_c"])
    [('B_1', 'B_2', 'B_3')]
    >>> find_vector_triplets(["J_s0_1", "J_s0_2", "J_s0_3", "J_s1_1"])
    [('J_s0_1', 'J_s0_2', 'J_s0_3')]
    """
    groups: dict[str, dict[int, str]] = {}
    for name in field_names:
        parts = vector_component(name)
        if parts is not None:
            base, component = parts
            groups.setdefault(base, {})[component] = name
    return [(g[1], g[2], g[3]) for g in groups.values() if len(g) == 3]

`identity_transform(frame)` ¶

Return a no-op transform within a single frame.

Parameters:

Name	Type	Description	Default
`frame`	`str`	Frame name used as both source and target.	required

Returns:

Type	Description
`FrameTransform`	Transform with no translation, rotation, or scaling, for which `is_identity` is `True`.

Examples:

>>> identity_transform("sim").is_identity
True

Source code in src/pypic/coordinates/transforms.py

def identity_transform(frame: str) -> FrameTransform:
    """Return a no-op transform within a single frame.

    Parameters
    ----------
    frame : str
        Frame name used as both source and target.

    Returns
    -------
    FrameTransform
        Transform with no translation, rotation, or scaling, for which
        ``is_identity`` is ``True``.

    Examples
    --------
    >>> identity_transform("sim").is_identity
    True
    """
    return FrameTransform(source_frame=frame, target_frame=frame)

`resolve_transform(source_frame, target_frame, transforms)` ¶

Find or compose a transform from source_frame to target_frame.

Tries direct lookup first, then searches for a one-hop chain (A→B→C) through intermediate frames. Also checks inverse transforms (if A→B is registered, B→A is available via inverse).

Parameters:

Name	Type	Description	Default
`source_frame`	`str`	Current frame name.	required
`target_frame`	`str`	Desired frame name.	required
`transforms`	`dict[str, FrameTransform]`	Registry mapping target frame names to transforms.	required

Returns:

Type	Description
`FrameTransform`

Raises:

Type	Description
`ValueError`	If no transform path exists.

Examples:

>>> t = FrameTransform("sim", "GSM", origin=(1.0, 0.0, 0.0))
>>> resolve_transform("sim", "GSM", {"GSM": t}).target_frame
'GSM'

Source code in src/pypic/coordinates/transforms.py

def resolve_transform(
    source_frame: str,
    target_frame: str,
    transforms: dict[str, FrameTransform],
) -> FrameTransform:
    """Find or compose a transform from *source_frame* to *target_frame*.

    Tries direct lookup first, then searches for a one-hop chain
    (A→B→C) through intermediate frames. Also checks inverse
    transforms (if A→B is registered, B→A is available via inverse).

    Parameters
    ----------
    source_frame : str
        Current frame name.
    target_frame : str
        Desired frame name.
    transforms : dict[str, FrameTransform]
        Registry mapping target frame names to transforms.

    Returns
    -------
    FrameTransform

    Raises
    ------
    ValueError
        If no transform path exists.

    Examples
    --------
    >>> t = FrameTransform("sim", "GSM", origin=(1.0, 0.0, 0.0))
    >>> resolve_transform("sim", "GSM", {"GSM": t}).target_frame
    'GSM'
    """
    if source_frame == target_frame:
        return identity_transform(source_frame)

    # Build a full map of all available directed edges
    edges: dict[tuple[str, str], FrameTransform] = {}
    for t in transforms.values():
        edges[(t.source_frame, t.target_frame)] = t
        edges[(t.target_frame, t.source_frame)] = t.inverse()

    # Direct lookup
    if (source_frame, target_frame) in edges:
        return edges[(source_frame, target_frame)]

    # One-hop chain: source → intermediate → target
    for (src, mid), first in list(edges.items()):
        if src != source_frame:
            continue
        if (mid, target_frame) in edges:
            return compose_transforms(first, edges[(mid, target_frame)])

    available = sorted(
        {t.source_frame for t in transforms.values()}
        | {t.target_frame for t in transforms.values()}
    )
    msg = (
        f"No transform path from {source_frame!r} to {target_frame!r}. "
        f"Available frames: {available}"
    )
    raise ValueError(msg)

`rotate_pressure_tensor(p11, p22, p33, p12, p13, p23, rotation)` ¶

Rotate a symmetric pressure tensor by a rotation matrix.

\[P'_{ij} = \sum_{k,l} R_{ik} \, R_{jl} \, P_{kl}\]

The trace \(P_{11} + P_{22} + P_{33}\) is invariant.

Parameters:

Name	Type	Description	Default
`p11`	`FloatArray`	Independent components of the symmetric tensor.	required
`p22`	`FloatArray`	Independent components of the symmetric tensor.	required
`p33`	`FloatArray`	Independent components of the symmetric tensor.	required
`p12`	`FloatArray`	Independent components of the symmetric tensor.	required
`p13`	`FloatArray`	Independent components of the symmetric tensor.	required
`p23`	`FloatArray`	Independent components of the symmetric tensor.	required
`rotation`	`FloatArray`	Shape `(3, 3)` rotation matrix.	required

Returns:

Type	Description
`tuple[FloatArray, ...]`	`(P_11', P_22', P_33', P_12', P_13', P_23')` in the rotated frame.

Examples:

>>> import numpy as np
>>> I = np.eye(3)
>>> p = rotate_pressure_tensor(
...     np.array([1.0]), np.array([2.0]), np.array([3.0]),
...     np.array([0.0]), np.array([0.0]), np.array([0.0]), I,
... )
>>> [float(x[0]) for x in p]
[1.0, 2.0, 3.0, 0.0, 0.0, 0.0]

Source code in src/pypic/coordinates/transforms.py

def rotate_pressure_tensor(
    p11: FloatArray,
    p22: FloatArray,
    p33: FloatArray,
    p12: FloatArray,
    p13: FloatArray,
    p23: FloatArray,
    rotation: FloatArray,
) -> tuple[FloatArray, FloatArray, FloatArray, FloatArray, FloatArray, FloatArray]:
    r"""Rotate a symmetric pressure tensor by a rotation matrix.

    $$P'_{ij} = \sum_{k,l} R_{ik} \, R_{jl} \, P_{kl}$$

    The trace $P_{11} + P_{22} + P_{33}$ is invariant.

    Parameters
    ----------
    p11, p22, p33, p12, p13, p23 : FloatArray
        Independent components of the symmetric tensor.
    rotation : FloatArray
        Shape ``(3, 3)`` rotation matrix.

    Returns
    -------
    tuple[FloatArray, ...]
        ``(P_11', P_22', P_33', P_12', P_13', P_23')`` in the rotated frame.

    Examples
    --------
    >>> import numpy as np
    >>> I = np.eye(3)
    >>> p = rotate_pressure_tensor(
    ...     np.array([1.0]), np.array([2.0]), np.array([3.0]),
    ...     np.array([0.0]), np.array([0.0]), np.array([0.0]), I,
    ... )
    >>> [float(x[0]) for x in p]
    [1.0, 2.0, 3.0, 0.0, 0.0, 0.0]
    """
    r = rotation
    # Build the full 3×3 symmetric tensor per grid point, then rotate.
    # P'_ij = sum_kl R_ik R_jl P_kl
    # Expand all 9 components (P is symmetric: P_21=P_12, P_31=P_13, P_32=P_23)
    p = [[p11, p12, p13], [p12, p22, p23], [p13, p23, p33]]

    def _component(i: int, j: int) -> FloatArray:
        result = r[i, 0] * r[j, 0] * p[0][0]
        for k in range(3):
            for el in range(3):
                if k == 0 and el == 0:
                    continue
                result = result + r[i, k] * r[j, el] * p[k][el]
        return cast("FloatArray", result)

    return (
        _component(0, 0),
        _component(1, 1),
        _component(2, 2),
        _component(0, 1),
        _component(0, 2),
        _component(1, 2),
    )

`rotate_vector_components(v1, v2, v3, rotation)` ¶

Rotate vector field components by a 3×3 rotation matrix.

\[v'_i = \sum_j R_{ij} \, v_j\]

Parameters:

Name	Type	Description	Default
`v1`	`FloatArray`	Vector field components (arbitrary shape, must match).	required
`v2`	`FloatArray`	Vector field components (arbitrary shape, must match).	required
`v3`	`FloatArray`	Vector field components (arbitrary shape, must match).	required
`rotation`	`FloatArray`	Shape `(3, 3)` rotation matrix.	required

Returns:

Type	Description
`tuple[FloatArray, FloatArray, FloatArray]`

Examples:

>>> import numpy as np
>>> v1, v2, v3 = np.array([1.0]), np.array([0.0]), np.array([0.0])
>>> R = np.array([[0, 0, 1], [1, 0, 0], [0, 1, 0]], dtype=float)
>>> r1, r2, r3 = rotate_vector_components(v1, v2, v3, R)
>>> float(r1[0]), float(r2[0]), float(r3[0])
(0.0, 1.0, 0.0)

Source code in src/pypic/coordinates/transforms.py

def rotate_vector_components(
    v1: FloatArray,
    v2: FloatArray,
    v3: FloatArray,
    rotation: FloatArray,
) -> tuple[FloatArray, FloatArray, FloatArray]:
    r"""Rotate vector field components by a 3×3 rotation matrix.

    $$v'_i = \sum_j R_{ij} \, v_j$$

    Parameters
    ----------
    v1, v2, v3 : FloatArray
        Vector field components (arbitrary shape, must match).
    rotation : FloatArray
        Shape ``(3, 3)`` rotation matrix.

    Returns
    -------
    tuple[FloatArray, FloatArray, FloatArray]

    Examples
    --------
    >>> import numpy as np
    >>> v1, v2, v3 = np.array([1.0]), np.array([0.0]), np.array([0.0])
    >>> R = np.array([[0, 0, 1], [1, 0, 0], [0, 1, 0]], dtype=float)
    >>> r1, r2, r3 = rotate_vector_components(v1, v2, v3, R)
    >>> float(r1[0]), float(r2[0]), float(r3[0])
    (0.0, 1.0, 0.0)
    """
    r = rotation
    v1_new = r[0, 0] * v1 + r[0, 1] * v2 + r[0, 2] * v3
    v2_new = r[1, 0] * v1 + r[1, 1] * v2 + r[1, 2] * v3
    v3_new = r[2, 0] * v1 + r[2, 1] * v2 + r[2, 2] * v3
    return v1_new, v2_new, v3_new

Coordinates & Operators¶

coordinates ¶

CoordinateGeometry dataclass ¶

metric_factors(x1, x2, x3) ¶

GeometryType ¶

FrameTransform dataclass ¶

rotation_matrix property ¶

is_identity property ¶

inverse() ¶

curl(f1, f2, f3, d1, d2, d3=None, *, geometry=GeometryType.CARTESIAN) ¶

divergence(f1, f2, f3, d1, d2, d3=None, *, geometry=GeometryType.CARTESIAN) ¶

gradient(f, d1, d2, d3=None, *, geometry=GeometryType.CARTESIAN) ¶

compose_transforms(first, second) ¶

find_pressure_tensor_groups(field_names) ¶

find_vector_triplets(field_names) ¶

identity_transform(frame) ¶

resolve_transform(source_frame, target_frame, transforms) ¶

rotate_pressure_tensor(p11, p22, p33, p12, p13, p23, rotation) ¶

rotate_vector_components(v1, v2, v3, rotation) ¶

`coordinates` ¶

`CoordinateGeometry` `dataclass` ¶

`metric_factors(x1, x2, x3)` ¶

`GeometryType` ¶

`FrameTransform` `dataclass` ¶

`rotation_matrix` `property` ¶

`is_identity` `property` ¶

`inverse()` ¶

`curl(f1, f2, f3, d1, d2, d3=None, *, geometry=GeometryType.CARTESIAN)` ¶

`divergence(f1, f2, f3, d1, d2, d3=None, *, geometry=GeometryType.CARTESIAN)` ¶

`gradient(f, d1, d2, d3=None, *, geometry=GeometryType.CARTESIAN)` ¶

`compose_transforms(first, second)` ¶

`find_pressure_tensor_groups(field_names)` ¶

`find_vector_triplets(field_names)` ¶

`identity_transform(frame)` ¶

`resolve_transform(source_frame, target_frame, transforms)` ¶

`rotate_pressure_tensor(p11, p22, p33, p12, p13, p23, rotation)` ¶

`rotate_vector_components(v1, v2, v3, rotation)` ¶