# plasma_plots.analysis

*module*

Numerical diagnostics returning values and labeled arrays, without rendering.

Fits (growth and damping rates, convergence orders, dispersion branches), norms and errors,
volume integrals and field energies on mapped domains, vector calculus on the logical grid,
velocity moments and orbit diagnostics. Most functions are also available as accessor methods,
e.g. ``array.plasma.analysis.error(...)``.

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L1-L1)

## plasma_plots.analysis.HELICITIES

*attribute* · *module attribute*

```python
HELICITIES = {'QA': (1, 0), 'QP': (0, 1), 'QH': (1, None)}
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L2282-L2282)

## plasma_plots.analysis.ORBIT_CLASSES

*attribute* · *module attribute*

```python
ORBIT_CLASSES = {0: 'passing', 1: 'trapped', -1: 'lost'}
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L627-L627)

## plasma_plots.analysis.VELOCITY_DIMS

*attribute* · *module attribute*

```python
VELOCITY_DIMS = ('v1', 'v2', 'v3')
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L449-L449)

## plasma_plots.analysis.BranchFit

*class* · *dataclass*

```python
class BranchFit
```

A dispersion branch fitted as ``omega = velocity * k``, from [`fit_dispersion_branches()`][plasma_plots.analysis.fit_dispersion_branches].

**Attributes**

- `velocity` (`float`) — The fitted slope, the phase velocity of the branch.
- `k` (`numpy.ndarray`) — The wavenumbers of the ridge points used in the fit.
- `omega` (`numpy.ndarray`) — The angular frequencies of the ridge points, one per ``k``.

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L695-L711)

### plasma_plots.analysis.BranchFit.k

*attribute* · *instance attribute*

```python
k: np.ndarray
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L710-L710)

### plasma_plots.analysis.BranchFit.omega

*attribute* · *instance attribute*

```python
omega: np.ndarray
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L711-L711)

### plasma_plots.analysis.BranchFit.velocity

*attribute* · *instance attribute*

```python
velocity: float
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L709-L709)

## plasma_plots.analysis.ConvergenceFit

*class* · *dataclass*

```python
class ConvergenceFit
```

A power law ``error = constant * size**order``, from [`convergence_order()`][plasma_plots.analysis.convergence_order].

**Attributes**

- `order` (`float`) — The fitted exponent: negative when the error shrinks as the size grows (e.g. points per cell), positive when it shrinks with the size (e.g. ``dt``).
- `constant` (`float`) — The fitted prefactor.
- `sizes` (`numpy.ndarray`) — The sizes of the valid (finite, positive) samples used in the fit.
- `fitted` (`numpy.ndarray`) — The fitted errors at ``sizes``, ready to plot over the data.

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L67-L87)

### plasma_plots.analysis.ConvergenceFit.constant

*attribute* · *instance attribute*

```python
constant: float
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L85-L85)

### plasma_plots.analysis.ConvergenceFit.fitted

*attribute* · *instance attribute*

```python
fitted: np.ndarray
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L87-L87)

### plasma_plots.analysis.ConvergenceFit.order

*attribute* · *instance attribute*

```python
order: float
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L84-L84)

### plasma_plots.analysis.ConvergenceFit.sizes

*attribute* · *instance attribute*

```python
sizes: np.ndarray
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L86-L86)

## plasma_plots.analysis.FitResult

*class* · *dataclass*

```python
class FitResult
```

An exponential fit ``exp(rate t + intercept)``, from [`growth_rate()`][plasma_plots.analysis.growth_rate] or [`damping_rate()`][plasma_plots.analysis.damping_rate].

**Attributes**

- `rate` (`float`) — The fitted rate: positive for growth, negative for damping. With ``GrowthFit.amplitude_from_quadratic`` it is the rate of the amplitude.
- `intercept` (`float`) — The fitted intercept of the logarithm (of the amplitude, with ``GrowthFit.amplitude_from_quadratic``).
- `time` (`numpy.ndarray`) — The times of the samples that were used in the fit.
- `fitted` (`numpy.ndarray`) — The fitted curve at ``time``, in the units of the fitted series (squared again with ``GrowthFit.amplitude_from_quadratic``), ready to plot over the data.

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L42-L64)

### plasma_plots.analysis.FitResult.fitted

*attribute* · *instance attribute*

```python
fitted: np.ndarray
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L64-L64)

### plasma_plots.analysis.FitResult.intercept

*attribute* · *instance attribute*

```python
intercept: float
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L62-L62)

### plasma_plots.analysis.FitResult.rate

*attribute* · *instance attribute*

```python
rate: float
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L61-L61)

### plasma_plots.analysis.FitResult.time

*attribute* · *instance attribute*

```python
time: np.ndarray
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L63-L63)

## plasma_plots.analysis.GrowthFit

*class* · *dataclass*

```python
class GrowthFit
```

Configuration for an exponential growth-rate fit.

Passed to [`growth_rate()`][plasma_plots.analysis.growth_rate] and [`damping_rate()`][plasma_plots.analysis.damping_rate].

**Attributes**

- `window` (`(float or None, float or None)`) — The time interval ``(start, end)`` of the samples that are fitted; ``None`` means the first or last time. The two ends may be given in either order. Default: every sample.
- `amplitude_from_quadratic` (`bool`) — Whether the series is quadratic in an amplitude (e.g. an energy). The fit then uses the square root of the samples, so that ``rate`` is the growth rate of the amplitude, half that of the series itself. Default: ``False``.

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L21-L39)

### plasma_plots.analysis.GrowthFit.amplitude_from_quadratic

*attribute* · *class attribute* · *instance attribute*

```python
amplitude_from_quadratic: bool = False
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L39-L39)

### plasma_plots.analysis.GrowthFit.window

*attribute* · *class attribute* · *instance attribute*

```python
window: tuple[float | None, float | None] = (None, None)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L38-L38)

## plasma_plots.analysis.OscillationFit

*class* · *dataclass*

```python
class OscillationFit
```

The frequency of an oscillating time series, from [`oscillation_frequency()`][plasma_plots.analysis.oscillation_frequency].

**Attributes**

- `omega` (`float`) — The angular frequency ``2π / period``.
- `period` (`float`) — The period, from a straight-line fit through the times of successive zero crossings (half periods apart) or peaks (a period apart).
- `times` (`numpy.ndarray`) — The times of the crossings or peaks that were used.
- `method` (`str`) — ``"zero_crossings"`` or ``"peaks"``.

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L258-L278)

### plasma_plots.analysis.OscillationFit.method

*attribute* · *instance attribute*

```python
method: str
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L278-L278)

### plasma_plots.analysis.OscillationFit.omega

*attribute* · *instance attribute*

```python
omega: float
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L275-L275)

### plasma_plots.analysis.OscillationFit.period

*attribute* · *instance attribute*

```python
period: float
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L276-L276)

### plasma_plots.analysis.OscillationFit.times

*attribute* · *instance attribute*

```python
times: np.ndarray
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L277-L277)

## plasma_plots.analysis.boozer_spectrum

*function*

```python
def boozer_spectrum(data: xr.DataArray, *, top: int | None = None, angles: str = 'boozer') -> xr.DataArray
```

The Boozer harmonics ``B_mn`` of a quantity on flux surfaces, real amplitudes over the radius.

The Fourier amplitudes of [`plasma_plots.spectral.mode_spectrum()`][plasma_plots.spectral.mode_spectrum] over the Boozer angles,
combined into the real amplitudes of [`plasma_plots.spectral.mode_amplitudes()`][plasma_plots.spectral.mode_amplitudes]: a field
``Σ B_mn cos(m θ_B + n ζ_B)`` gives ``B_mn`` at ``(m, n)``, with ``n`` the full-torus mode
number (a multiple of nfp). The usual stellarator convention ``cos(m θ_B − n ζ_B)`` has the
opposite sign of ``n``. In Boozer angles, the spectrum of ``|B|`` is what the guiding-center
drifts see: a quasi-symmetric field has a single helicity in it (see
[`quasisymmetry_error()`][plasma_plots.analysis.quasisymmetry_error]).

**Parameters**

- `data` (`xarray.DataArray`) — The quantity, e.g. GVEC's ``mod_B`` on a Boozer grid (``from_gvec`` of ``state.evaluate_sfl(..., sfl="boozer")``), over ``(rho, theta_B, zeta_B)`` and any other dimensions.
- `top` (`int`) (default: `None`) — Keep only the ``top`` harmonics with the largest peak amplitude over the radius, strongest first. Default: all.
- `angles` (`(boozer, any)`) (default: `"boozer"`) — Require the Boozer angles ``theta_B``, ``zeta_B`` (default), or take the harmonics in whatever angles the field has (not a Boozer spectrum then; e.g. for a comparison).

**Returns**

- (`xarray.DataArray`) — The amplitudes over ``mode`` (with the coordinates ``m``, ``n`` and a label such as ``"(1, -5)"``) and the remaining dimensions, e.g. ``rho``; named ``B_mn``.

**Raises**

- `ValueError` — If the field is not over the Boozer angles (unless ``angles="any"``), or does not sample
them over full periods.

> **See Also**
>
> [`quasisymmetry_error()`][plasma_plots.analysis.quasisymmetry_error] : The symmetry-breaking part of the spectrum.
> [`plasma_plots.spectral_plots.plot_boozer_spectrum()`][plasma_plots.spectral_plots.plot_boozer_spectrum] : The harmonics over the radius.

**Examples**

```pycon
>>> B_mn = boozer_spectrum(boozer.mod_B, top=8)
>>> B_mn.sel(mode="(0, 0)")  # the mean field on each surface
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L2297-L2360)

## plasma_plots.analysis.bounce_period

*function*

```python
def bounce_period(orbits: xr.Dataset, *, v_par: str = 'v_par') -> xr.DataArray
```

The bounce period of each trapped marker.

Twice the mean time between reversals of its parallel velocity; the reversal times are
interpolated linearly between the saved samples.

**Parameters**

- `orbits` (`xarray.Dataset`) — The orbits product, with variables over ``(t, marker)``.
- `v_par` (`str`) (default: `'v_par'`) — The name of the parallel-velocity variable. Default: ``"v_par"``.

**Returns**

- (`xarray.DataArray`) — The period over ``marker``, named ``bounce_period``; NaN for markers with fewer than two reversals, e.g. passing ones.

> **See Also**
>
> [`classify_orbits()`][plasma_plots.analysis.classify_orbits]

**Examples**

```pycon
>>> periods = bounce_period(orbits)
>>> periods.where(classify_orbits(orbits) == 1).mean()
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L2227-L2276)

## plasma_plots.analysis.classify_orbits

*function*

```python
def classify_orbits(orbits: xr.Dataset, *, v_par: str = 'v_par') -> xr.DataArray
```

Classify each marker of an orbits product as passing (0), trapped (1) or lost (-1).

The same criteria as Struphy's ``post_process_orbit_classification``: a marker is trapped if
its parallel velocity ``v_par`` ever has the opposite sign to its initial one, and lost if at
any saved time every quantity is zero (how Struphy stores a marker that has left the domain).
Lost takes precedence over trapped.

**Parameters**

- `orbits` (`xarray.Dataset`) — The orbits product, with variables over ``(t, marker)``.
- `v_par` (`str`) (default: `'v_par'`) — The name of the parallel-velocity variable. Default: ``"v_par"``.

**Returns**

- (`xarray.DataArray`) — A ``(marker,)`` array of integer codes named ``classification``; the names are in ``attrs["flag_meanings"]`` and in [`ORBIT_CLASSES`][plasma_plots.analysis.ORBIT_CLASSES].

**Raises**

- `ValueError` — If ``orbits`` has no variable ``v_par``, or it is not over ``(t, marker)``.

> **See Also**
>
> [`bounce_period()`][plasma_plots.analysis.bounce_period], [`orbit_invariants()`][plasma_plots.analysis.orbit_invariants]

**Examples**

```pycon
>>> codes = classify_orbits(orbits)
>>> trapped = orbits.sel(marker=codes == 1)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L630-L692)

## plasma_plots.analysis.convergence_order

*function*

```python
def convergence_order(sizes, errors) -> ConvergenceFit | None
```

Fit ``error = constant * size**order`` in log-log space.

Only valid samples, where both the size and the error are finite and positive, are used.

**Parameters**

- `sizes` (`array_like of float`) — Typically a resolution (points per cell, coarser to finer) or a step size (``dt``).
- `errors` (`array_like of float`) — The corresponding, necessarily positive, error norms, e.g. from [`error()`][plasma_plots.analysis.error].

**Returns**

- (`ConvergenceFit or None`) — The fit. ``order`` is negative when the error shrinks as ``sizes`` grows (e.g. more points per cell), and positive when it shrinks as ``sizes`` shrinks (e.g. a smaller ``dt``). ``None`` with fewer than two valid (finite, positive) samples.

**Examples**

```pycon
>>> errors = [
...     run.evaluate("T").plasma.analysis.error(exact).isel(t=-1)
...     for run in runs
... ]
>>> convergence_order([16, 32, 64], errors).order
-2.01
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L90-L126)

## plasma_plots.analysis.critical_points

*function*

```python
def critical_points(data: xr.DataArray, *, refine: bool = True) -> xr.Dataset
```

The O-points (extrema) and X-points (saddles) of a 2-D flux function.

The critical points are where the gradient vanishes: every grid cell whose corners carry
both signs of each gradient component (central differences, wrapping around a periodic
direction) is a candidate, and the zero of the bilinear interpolant of the gradient inside
it, found by Newton's method, is the point (a candidate whose zero lies outside the cell is
dropped). The Hessian of second differences, interpolated to the point, tells the kind:
``det H > 0`` an O-point (a maximum when its trace is negative, else a minimum),
``det H < 0`` an X-point. A time dimension is kept: the points are found at every time.

**Parameters**

- `data` (`xarray.DataArray`) — The flux function, e.g. [`flux_function()`][plasma_plots.analysis.flux_function] of a 2-D magnetic field, over two logical directions (a third with a single point is fine) and optionally ``t``. A direction wraps around when its coordinate has a ``period`` attribute or its physical coordinates close (see [`plasma_plots.arrays.periodicity()`][plasma_plots.arrays.periodicity]); a periodic box without either is treated as bounded, and loses the points on its edges.
- `refine` (`bool`) (default: `True`) — Locate the zero within the cell by Newton's method (the cell's centre otherwise). Default: True.

**Returns**

- (`xarray.Dataset`) — Over ``point`` (and ``t``): the logical coordinates of each point (named as in the field), ``X`` and ``Y`` when the field has physical coordinates (interpolated), ``value`` (the flux there), ``kind`` (``"O"`` or ``"X"``) and ``sign`` (``+1`` for a maximum, ``−1`` for a minimum, ``0`` for a saddle). O-points first, then X-points, each sorted by value; padded with NaN (and ``""``) when the number of points varies in time. The attrs name the ``plane`` and the ``periods`` of its directions (NaN when bounded).

**Raises**

- `ValueError` — If the field is not over exactly two logical directions with more than one point.

> **See Also**
>
> [`reconnected_flux()`][plasma_plots.analysis.reconnected_flux] : The flux between an O- and an X-point over time.
> [`plasma_plots.plotting.plot_critical_points()`][plasma_plots.plotting.plot_critical_points] : The points over the flux contours.

**Examples**

```pycon
>>> points = critical_points(flux_function(B))
>>> points.isel(t=-1).to_dataframe()
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L2506-L2732)

## plasma_plots.analysis.curl

*function*

```python
def curl(vector: xr.DataArray, *, components: str = 'cartesian', domain=None) -> xr.DataArray
```

The curl of a vector field on a mapped domain, in Cartesian components ``(x, y, z)``.

Derivatives as for [`divergence()`][plasma_plots.analysis.divergence]. E.g. the current ``J = ∇ × B``, or the vorticity of a
flow; for a 2-D field in the ``x``-``y`` plane only the ``z`` component is non-zero.

**Parameters**

- `vector` (`xarray.DataArray`) — The vector field, with a ``component`` dimension of size 3 and the dimensions ``eta1``, ``eta2``, ``eta3``.
- `components` (`(cartesian, contravariant)`) (default: `"cartesian"`) — What the components are: ``"cartesian"`` (default) ``(x, y, z)``, or ``"contravariant"`` components, which are pushed forward first.
- `domain` (`struphy domain`) (default: `None`) — The mapping (``out.domain``), for the exact Jacobian. Default: from the ``X``, ``Y``, ``Z`` coordinates.

**Returns**

- (`xarray.DataArray`) — The curl, with a ``component`` dimension ``(x, y, z)`` (coordinates 0, 1, 2) first, named ``curl_<name>``.

**Raises**

- `ValueError` — As for [`divergence()`][plasma_plots.analysis.divergence].

> **See Also**
>
> [`gradient()`][plasma_plots.analysis.gradient], [`divergence()`][plasma_plots.analysis.divergence]

**Examples**

```pycon
>>> J = curl(B)  # the current, in Cartesian components
>>> vorticity = curl(u).sel(component=2)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L1854-L1906)

## plasma_plots.analysis.cylindrical_components

*function*

```python
def cylindrical_components(vector: xr.DataArray) -> xr.DataArray
```

Cartesian components rotated to cylindrical ones ``(R, phi, Z)`` about the ``Z`` axis.

The rotation is by the toroidal angle ``φ = atan2(Y, X)`` at every point of the field.

**Parameters**

- `vector` (`xarray.DataArray`) — The field in Cartesian components, with a ``component`` dimension of size 3 and the ``X``, ``Y``, ``Z`` coordinates.

**Returns**

- (`xarray.DataArray`) — The components ``v_R``, ``v_φ``, ``v_Z``, with ``component`` coordinates ``"R"``, ``"phi"``, ``"Z"`` and the dimensions of ``vector``.

**Raises**

- `ValueError` — If ``vector`` has no ``component`` dimension of size 3.

> **See Also**
>
> [`toroidal_components()`][plasma_plots.analysis.toroidal_components]

**Examples**

```pycon
>>> E_R = cylindrical_components(-gradient(phi)).sel(component="R")
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L2005-L2057)

## plasma_plots.analysis.damping_rate

*function*

```python
def damping_rate(data: xr.DataArray, fit: GrowthFit | None = None) -> FitResult | None
```

Fit ``exp(rate*t + intercept)`` to the envelope of an oscillating time series.

Use this for signals such as the field energy in Landau damping, where [`growth_rate()`][plasma_plots.analysis.growth_rate] on
the raw series would fit the oscillation. The envelope is the series' local maxima, see
[`envelope()`][plasma_plots.analysis.envelope].

**Parameters**

- `data` (`xarray.DataArray`) — The oscillating time series, with ``t`` as its only dimension.
- `fit` (`GrowthFit`) (default: `None`) — ``fit.window`` restricts the peaks that are used; with ``fit.amplitude_from_quadratic`` the rate of the amplitude is returned. Default: ``GrowthFit()``.

**Returns**

- (`FitResult or None`) — The fit to the peaks; the rate is negative for damping. ``None`` with fewer than two finite, positive peaks within the window.

**Raises**

- `ValueError` — If ``data`` has dimensions other than ``t``.

> **See Also**
>
> [`growth_rate()`][plasma_plots.analysis.growth_rate], [`envelope()`][plasma_plots.analysis.envelope]

**Examples**

```pycon
>>> damping_rate(
...     out.scalars["en_E"], GrowthFit(amplitude_from_quadratic=True)
... ).rate
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L219-L255)

## plasma_plots.analysis.divergence

*function*

```python
def divergence(vector: xr.DataArray, *, components: str = 'cartesian', domain=None) -> xr.DataArray
```

The divergence ``sum_a d v_a / d x_a`` of a vector field on a mapped domain.

Derivatives as in [`gradient()`][plasma_plots.analysis.gradient]: spectral around periodic angles; a flat direction (2-D
run) is left out. E.g. a ``div B`` check of an MHD run.

**Parameters**

- `vector` (`xarray.DataArray`) — The vector field, with a ``component`` dimension of size 3 and the dimensions ``eta1``, ``eta2``, ``eta3``.
- `components` (`(cartesian, contravariant)`) (default: `"cartesian"`) — What the components are, as for the 3-D views: ``"cartesian"`` (default) ``(x, y, z)``, or ``"contravariant"`` components, which are pushed forward first.
- `domain` (`struphy domain`) (default: `None`) — The mapping (``out.domain``), for the exact Jacobian. Default: from the ``X``, ``Y``, ``Z`` coordinates.

**Returns**

- (`xarray.DataArray`) — The divergence, over the field's dimensions without ``component``, named ``div_<name>``.

**Raises**

- `ValueError` — If ``vector`` has no ``component`` dimension of size 3, ``components`` is unknown, or as
for [`gradient()`][plasma_plots.analysis.gradient].

> **See Also**
>
> [`gradient()`][plasma_plots.analysis.gradient], [`curl()`][plasma_plots.analysis.curl]

**Examples**

```pycon
>>> div_B = divergence(B)  # should stay at round-off
>>> norm(div_B).plasma.plot.timeseries()
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L1802-L1851)

## plasma_plots.analysis.drift

*function*

```python
def drift(data: xr.DataArray, *, ref=None) -> xr.DataArray
```

Signed deviation from an explicit reference or the first time sample.

**Parameters**

- `data` (`xarray.DataArray`) — The array, with a ``t`` dimension.
- `ref` (`float or array_like or xarray.DataArray`) (default: `None`) — The reference, broadcast against ``data``. Default: ``data`` at the first time sample.

**Returns**

- (`xarray.DataArray`) — ``data - ref``, with the attributes of ``data`` and the label ``"... drift"``.

> **See Also**
>
> [`relative_error()`][plasma_plots.analysis.relative_error] : The absolute relative deviation.

**Examples**

```pycon
>>> drift(out.scalars["en_tot"]).plasma.plot.timeseries()
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L418-L446)

## plasma_plots.analysis.envelope

*function*

```python
def envelope(data: xr.DataArray) -> xr.DataArray
```

Local maxima of a time series: the interior samples not smaller than their neighbours.

A sample is a peak when it is larger than the previous sample and not smaller than the next;
the first and last samples never are.

**Parameters**

- `data` (`xarray.DataArray`) — The time series, with ``t`` as its only dimension.

**Returns**

- (`xarray.DataArray`) — The peaks of ``data``, a selection along ``t`` with the attributes and coordinates kept.

**Raises**

- `ValueError` — If ``data`` has dimensions other than ``t``.

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L189-L216)

## plasma_plots.analysis.error

*function*

```python
def error(data: xr.DataArray, exact, *, norm: str = 'rms', relative: bool = False, dims=None, weighted: bool = False, domain=None, args=None) -> xr.DataArray
```

The error of ``data`` against an exact solution, as a function of the other dimensions.

Unweighted, the norms are means over the grid points (``l1``, and ``rms`` = ``l2`` its root
mean square), i.e. integrals over the logical unit cube. With ``weighted`` the norms over the
logical grid are integrals over the physical volume (with ``|√g|``, from ``domain`` or the
``X``, ``Y``, ``Z`` coordinates; see [`volume_integral()`][plasma_plots.analysis.volume_integral]), and ``rms`` is divided by that
volume: the proper L2 error on a mapped domain, where a plain mean over grid points is not.

**Parameters**

- `data` (`xarray.DataArray`) — The numerical solution.
- `exact` (`(callable, array_like, number or xarray.DataArray)`) — The exact solution: an array (aligned with ``data``, or broadcast to it) or a function of its coordinates, evaluated with [`evaluate_on()`][plasma_plots.analysis.evaluate_on], e.g. ``lambda x, y, z, t: ...``.
- `norm` (`(rms, max, l1, l2, pointwise)`) (default: `"rms"`) — ``"pointwise"`` is the difference ``data − exact`` itself; ``"max"`` the largest absolute difference; ``"l1"`` the mean (or, weighted, the integral) of ``|data − exact|``; ``"l2"`` the square root of the mean (or integral) of ``|data − exact|²``; ``"rms"`` (default) as ``"l2"``, divided by the volume when weighted.
- `relative` (`bool`) (default: `False`) — Divide by the same norm of the exact solution; for ``"pointwise"``, by the largest ``|exact|`` over the whole array. Default: ``False``.
- `dims` (`str or sequence of str`) (default: `None`) — The dimensions the norm is taken over. Default: every dimension but ``t``. Weighted norms need exactly ``eta1``, ``eta2``, ``eta3``.
- `weighted` (`bool`) (default: `False`) — Integrate over the physical volume instead of averaging over the grid points (no effect on ``"max"`` and ``"pointwise"``). Default: ``False``.
- `domain` (`struphy domain`) (default: `None`) — The mapping (``out.domain``), for the exact ``|√g|`` of weighted norms. Default: from the ``X``, ``Y``, ``Z`` coordinates.
- `args` (`sequence of str`) (default: `None`) — The coordinates passed to a callable ``exact``, as in [`evaluate_on()`][plasma_plots.analysis.evaluate_on]. Default: ``X``, ``Y``, ``Z`` (or the logical dimensions), then ``t``.

**Returns**

- (`xarray.DataArray`) — The error as a function of the dimensions not in ``dims`` (typically ``t``); for ``"pointwise"``, an array like ``data``. Labeled e.g. ``"relative rms error of ..."``.

**Raises**

- `ValueError` — If ``norm`` is unknown, or a weighted norm is not over ``eta1``, ``eta2``, ``eta3``.

> **See Also**
>
> [`convergence_order()`][plasma_plots.analysis.convergence_order] : The order of a series of errors.

**Examples**

```pycon
>>> error(T, exact, relative=True)  # relative RMS error over time
>>> error(T, exact, norm="max")  # largest pointwise error
>>> # √∫|u − u_exact|² dV
>>> error(u, lambda x, y, z, t: np.sin(x - t), norm="l2", weighted=True)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L1569-L1696)

## plasma_plots.analysis.evaluate_on

*function*

```python
def evaluate_on(data: xr.DataArray, function, args=None) -> xr.DataArray
```

``function`` evaluated on the coordinates of ``data``, broadcast to its shape.

E.g. an exact solution ``exact(x, y, z, t)`` on the points of a field, or ``exact(eta1, t)``
on a 1-D profile.

**Parameters**

- `data` (`xarray.DataArray`) — The array whose coordinates and shape are used.
- `function` (`callable`) — Called with the coordinates named by ``args``, positionally, as xarray.DataArrays; returns an array (or number) broadcastable to ``data``.
- `args` (`sequence of str`) (default: `None`) — The coordinates passed positionally, in order. Default: the physical ``X``, ``Y``, ``Z`` (if attached, else the logical dimensions ``eta1``... that ``data`` has), followed by ``t`` if ``data`` has it.

**Returns**

- (`xarray.DataArray`) — The values, broadcast to the shape of ``data`` with its dimensions first.

**Raises**

- `ValueError` — If one of ``args`` is not a coordinate of ``data``.

**Examples**

```pycon
>>> exact = evaluate_on(u, lambda x, y, z, t: np.sin(x - t))
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L1520-L1566)

## plasma_plots.analysis.field_energy

*function*

```python
def field_energy(data: xr.DataArray, *, form: int | str | None = None, weight=None, domain=None, normalization: float = 1.0, quadrature=None) -> xr.DataArray
```

The quadratic energy ``α ½ ∫ w ωᵀ A ω dη``, as a function of time.

``form`` says what ``data`` holds, and sets the metric factor ``A`` (as struphy's mass
matrices do, so that e.g. LinearMHD's ``en_U`` is ``field_energy(u, form=2, weight=n0)``):

========================  ===========================================  ==============
``form``                  data                                         ``A``
========================  ===========================================  ==============
``None`` (default)        a function, or Cartesian vector components    ``|√g|``
``0``                     0-form                                        ``|√g|``
``1``                     1-form components                             ``G⁻¹ |√g|``
``2``                     2-form components                             ``G / |√g|``
``3``                     3-form                                        ``1 / |√g|``
``"v"``                   contravariant vector components               ``G |√g|``
========================  ===========================================  ==============

Here ``G = Jᵀ J`` is the metric of the mapping and ``|√g|`` its Jacobian determinant. The
energy of a filtered field, e.g. from [`filter_time()`][plasma_plots.spectral.filter_time], measures
how much of the energy is in that mode. A spline field squared is integrated exactly only
with enough points per element: evaluate it at Gauss points for accurate energies.

**Parameters**

- `data` (`xarray.DataArray`) — The field, with the dimensions ``eta1``, ``eta2``, ``eta3``; vectors have a ``component`` dimension of size 3.
- `form` (`(None, 0, 1, 2, 3, v)`) (default: `None`) — What ``data`` holds, see the table. Default: ``None``, a function or Cartesian vector components.
- `weight` (`array_like`) (default: `None`) — The weight ``w`` over ``(eta1, eta2, eta3)`` (e.g. a background density ``n0``). Non-finite values (e.g. ``1/p0`` where ``p0`` vanishes on the boundary) are left out.
- `domain` (`struphy domain`) (default: `None`) — The mapping (``out.domain``), which gives the exact ``|√g|`` and ``G``. Without one, they come from the numerical Jacobian of the ``X``, ``Y``, ``Z`` coordinates (see [`plasma_plots.arrays.mapping_jacobian()`][plasma_plots.arrays.mapping_jacobian]).
- `normalization` (`float`) (default: `1.0`) — The prefactor ``α``. Default: 1.
- `quadrature` (`dict`) (default: `None`) — Explicit quadrature weights per logical dimension, as for [`volume_integral()`][plasma_plots.analysis.volume_integral].

**Returns**

- (`xarray.DataArray`) — The energy, as a function of the dimensions other than ``eta1``, ``eta2``, ``eta3`` and ``component`` (typically ``t``), labeled ``"energy of ..."``.

**Raises**

- `ValueError` — If ``form`` is not one of the above or does not match ``data`` (a vector form for a
scalar or the reverse), a vector does not have 3 components, or a logical dimension is
missing.

> **See Also**
>
> [`volume_integral()`][plasma_plots.analysis.volume_integral]

**Examples**

```pycon
>>> # ½ uᵀ M2n u, as LinearMHD's en_U
>>> field_energy(u_2form, form=2, weight=n0, domain=out.domain)
>>> field_energy(E)  # ½ ∫ |E|² dV
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L1275-L1401)

## plasma_plots.analysis.fit_dispersion_branches

*function*

```python
def fit_dispersion_branches(spectrum: xr.DataArray, *, n_branches: int, k_range: tuple[float, float] | None = None, noise_level: float = 0.5, order: int = 10) -> list[BranchFit]
```

Fit ``omega = v * k`` to each of ``n_branches`` straight ridges in an (omega, k) power spectrum.

The spectrum is as returned by [`power_spectrum()`][plasma_plots.analysis.power_spectrum], and there is no theoretical curve to
guide the search -- useful when several linear wave branches are excited at once and there is
nothing to search around yet. (For a single, possibly non-linear branch with a known
theoretical curve to guide the search instead, take the frequency of maximum power in a window
of ``spectrum`` around that curve at each k, rather than this blind approach.)

Only non-negative ``omega`` and ``k`` are scanned (a real signal's spectrum is symmetric under
``(k, omega) -> (-k, -omega)``, so every branch already appears on both sides of ``k = 0``). At
each remaining k, the local maxima of the spectrum along omega are found; a k column
contributes to the fit only where it has exactly ``n_branches`` maxima above ``noise_level``
times that column's own peak power, taken in increasing-omega order.

**Parameters**

- `spectrum` (`xarray.DataArray`) — The power spectrum, with dimensions ``omega`` and ``k``.
- `n_branches` (`int`) — The number of branches, at least 1.
- `k_range` (`(float, float)`) (default: `None`) — The interval of non-negative k's that are scanned. Default: ``(k.max() / 8, k.max() / 2)``, which in practice skips both the low-k region where branches have not yet separated, and the folded Nyquist edge.
- `noise_level` (`float`) (default: `0.5`) — Maxima count only above this fraction of the column's peak power. Default: 0.5.
- `order` (`int`) (default: `10`) — A local maximum must exceed every one of its ``order`` neighbors on both sides along omega (out-of-range neighbors are clipped to the edge sample, as in ``scipy.signal.argrelextrema``). Default: 10.

**Returns**

- (`list of BranchFit`) — One fit per branch, in increasing-omega order at the low-k end of ``k_range``; ``.velocity`` is the fitted slope, ``.k``/``.omega`` the ridge points used.

**Raises**

- `ValueError` — If ``spectrum`` lacks the ``omega`` or ``k`` dimension, ``n_branches`` is not positive, or
no k in ``k_range`` has exactly ``n_branches`` peaks above the noise level.

> **See Also**
>
> [`power_spectrum()`][plasma_plots.analysis.power_spectrum]

**Examples**

```pycon
>>> spectrum = power_spectrum(b.isel(eta2=0, eta3=0, component=0))
>>> [fit.velocity for fit in fit_dispersion_branches(spectrum, n_branches=2)]
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L730-L831)

## plasma_plots.analysis.flux_function

*function*

```python
def flux_function(vector: xr.DataArray) -> xr.DataArray
```

The flux function (or stream function) ``A`` of a 2-D, divergence-free in-plane field.

For a magnetic field ``B_x = ∂A/∂y``, ``B_y = −∂A/∂x`` (so ``B = ∇A × ẑ``); its contour
lines are the field lines, e.g. drawn over the current with a slice's ``contours_of``. For a
velocity it is the stream function. ``A`` is integrated along the grid lines (trapezoidal
rule), ``A(x, y) = ∫ B_x(x₀, y′) dy′ − ∫ B_y(x′, y) dx′``, and shifted to zero mean.

**Parameters**

- `vector` (`xarray.DataArray`) — The field in Cartesian components, with a ``component`` dimension of size 3 (only ``x`` and ``y`` are used). It must lie on a Cartesian grid in the ``x``-``y`` plane: exactly two logical directions with more than one point, physical ``X`` depending on one of them and ``Y`` on the other, e.g. a slab or periodic box.

**Returns**

- (`xarray.DataArray`) — ``A``, named ``flux_function``, as a function of every dimension of ``vector`` except ``component``.

**Raises**

- `ValueError` — If ``vector`` is not a 3-component field on such a grid.

**Examples**

```pycon
>>> A = flux_function(B.isel(t=-1))
>>> J = curl(B).isel(t=-1, component=2)
>>> J.plasma.plot.slice(overlays={"contours_of": A})
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L1909-L1994)

## plasma_plots.analysis.gradient

*function*

```python
def gradient(data: xr.DataArray, *, domain=None) -> xr.DataArray
```

The Cartesian gradient of a scalar field on the logical grid.

``∇f = J⁻ᵀ ∂f/∂η``, with the Jacobian ``J`` of the mapping from a struphy ``domain`` (exact)
or, without one, from the attached ``X``, ``Y``, ``Z`` coordinates. The logical derivatives
are spectral around periodic directions and second order elsewhere (see
[`plasma_plots.arrays.logical_derivative()`][plasma_plots.arrays.logical_derivative]). A direction with a single point (a 2-D run)
is left out: the result is then the gradient within the plane (with the pseudo-inverse of
the remaining Jacobian columns). For example, ``E = -gradient(phi)``, or the E × B velocity
``ẑ × ∇φ`` of a 2-D drift-wave model.

**Parameters**

- `data` (`xarray.DataArray`) — The scalar field, with the dimensions ``eta1``, ``eta2``, ``eta3`` (at least one with more than one point) and no ``component`` dimension.
- `domain` (`struphy domain`) (default: `None`) — The mapping (``out.domain``), for the exact Jacobian. Default: differentiate the ``X``, ``Y``, ``Z`` coordinates numerically.

**Returns**

- (`xarray.DataArray`) — The gradient, with a ``component`` dimension ``(x, y, z)`` (coordinates 0, 1, 2) in front of the field's own dimensions (e.g. ``t`` is kept), named ``grad_<name>``.

**Raises**

- `ValueError` — If a logical dimension is missing, ``data`` has a ``component`` dimension, no direction
has more than one point, or there is neither a ``domain`` nor ``X``, ``Y``, ``Z``.

> **See Also**
>
> [`divergence()`][plasma_plots.analysis.divergence], [`curl()`][plasma_plots.analysis.curl]

**Examples**

```pycon
>>> E = -gradient(phi)  # Cartesian components (x, y, z), over time
>>> E = -gradient(phi, domain=out.domain)  # with the exact Jacobian
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L1404-L1511)

## plasma_plots.analysis.growth_rate

*function*

```python
def growth_rate(data: xr.DataArray, fit: GrowthFit | None = None) -> FitResult | None
```

Fit ``exp(rate*t + intercept)`` to a time series, using only finite, positive samples.

The fit is a straight line through the logarithm of the samples within ``fit.window``.

**Parameters**

- `data` (`xarray.DataArray`) — The time series, with ``t`` as its only dimension (e.g. a field energy).
- `fit` (`GrowthFit`) (default: `None`) — The time window and whether the series is quadratic in the amplitude. Default: ``GrowthFit()``, every sample, the series itself.

**Returns**

- (`FitResult or None`) — The rate, the intercept, the times used and the fitted curve. ``None`` with fewer than two samples, or fewer than two finite, positive samples within the window.

**Raises**

- `ValueError` — If ``data`` has dimensions other than ``t``.

> **See Also**
>
> [`damping_rate()`][plasma_plots.analysis.damping_rate] : The same fit to the envelope of an oscillating series.

**Examples**

```pycon
>>> energy = out.scalars["en_E"]
>>> growth_rate(
...     energy, GrowthFit(window=(0.0, 5.0), amplitude_from_quadratic=True)
... ).rate
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L129-L186)

## plasma_plots.analysis.loss_map

*function*

```python
def loss_map(markers: xr.Dataset, *, x: str = 'v_par', y: str | None = None, t=0, absB=None) -> xr.Dataset
```

Each marker's initial phase-space position, whether it is lost, and when.

``x`` and ``y`` are variables of the Dataset, or ``"energy"``, ``"pitch"`` and ``"speed"``
from [`orbit_invariants()`][plasma_plots.analysis.orbit_invariants] (``absB`` needed for the first two).

**Parameters**

- `markers` (`xarray.Dataset`) — An orbits product over ``(t, marker)``.
- `x` (`str`) (default: `'v_par'`) — The first quantity. Default: ``"v_par"``.
- `y` (`str`) (default: `None`) — The second quantity. Default: ``"mu"``, or ``"v_perp"`` without ``mu``.
- `t` (`int or float`) (default: `0`) — The time of the plotted values: an integer position (default ``0``, the initial one) or a float nearest value.
- `absB` (`callable`) (default: `None`) — ``|B|(x, y, z)``, for ``"energy"`` and ``"pitch"`` (see [`orbit_invariants()`][plasma_plots.analysis.orbit_invariants]).

**Returns**

- (`xarray.Dataset`) — Over ``marker``: ``x`` and ``y`` (named as the quantities), ``lost`` (bool) and ``loss_time`` (the first time a marker is gone; NaN while confined).

**Raises**

- `ValueError` — If a quantity is unknown, or needs ``absB``.

> **See Also**
>
> [`lost_fraction()`][plasma_plots.analysis.lost_fraction] : The losses over time.
> [`plasma_plots.plotting.plot_loss_map()`][plasma_plots.plotting.plot_loss_map] : The plot.

**Examples**

```pycon
>>> losses = loss_map(orbits, x="energy", y="pitch", absB=absB)
>>> losses.where(losses.lost, drop=True)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L3072-L3161)

## plasma_plots.analysis.lost_fraction

*function*

```python
def lost_fraction(markers: xr.Dataset, *, weight: str | None = None) -> xr.DataArray
```

The fraction of markers that have left the domain, over time.

A marker is lost from the first time every saved quantity is zero (how Struphy stores it),
so the fraction never decreases. With ``weight`` the fraction is of the initial weights,
i.e. of the particles the markers represent.

**Parameters**

- `markers` (`xarray.Dataset`) — A marker Dataset over ``(t, marker)``, e.g. an orbits product.
- `weight` (`str`) (default: `None`) — Weigh each marker by its initial value of this variable, e.g. ``"weight"``. Default: count the markers.

**Returns**

- (`xarray.DataArray`) — The fraction over ``t``, between 0 and 1, named ``lost_fraction``.

> **See Also**
>
> [`loss_map()`][plasma_plots.analysis.loss_map] : Which markers are lost, and when.
> [`classify_orbits()`][plasma_plots.analysis.classify_orbits] : Lost markers as a class.

**Examples**

```pycon
>>> lost_fraction(orbits).plasma.plot.timeseries(logy=False)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L3024-L3069)

## plasma_plots.analysis.marker_density

*function*

```python
def marker_density(markers: xr.Dataset, *, dims=('eta1'), bins=32, weight: str | None = None, ranges=None) -> xr.DataArray
```

Bin the markers over position variables: the number (or weight) of markers per unit volume.

Without ``weight`` this is the sampling density, where the markers are (``s0`` in Struphy's
terms, the marker loading); with it the density the markers represent (the physical density
of a full-f run, the perturbation of a δf run). Markers that have left the domain are left
out.

**Parameters**

- `markers` (`xarray.Dataset`) — A marker Dataset with position variables over ``(t, marker)`` (or ``marker``).
- `dims` (`str or sequence of str`) (default: `('eta1')`) — The position variables to bin over, e.g. ``("eta1",)`` or ``("x", "y")``. Default: ``("eta1",)``.
- `bins` (`int or sequence of int`) (default: `32`) — The number of bins, for all or for each. Default: 32.
- `weight` (`str`) (default: `None`) — Weigh each marker by this variable, e.g. ``"weight"``. Default: count the markers.
- `ranges` (`dict`) (default: `None`) — ``{variable: (low, high)}`` bin ranges. Default: ``(0, 1)`` for the logical ``eta`` coordinates, the markers' extent otherwise.

**Returns**

- (`xarray.DataArray`) — The density over ``t`` and the binned variables (bin centres as coordinates, with the edges in ``attrs["edges"]``), named ``marker_density`` or ``weighted_density``.

**Raises**

- `ValueError` — If a variable is missing.

> **See Also**
>
> [`weight_statistics()`][plasma_plots.analysis.weight_statistics] : The weights' statistics.
> [`plasma_plots.plotting.plot_marker_density()`][plasma_plots.plotting.plot_marker_density] : Sampling against physical density.

**Examples**

```pycon
>>> sampling = marker_density(orbits, dims="eta1", bins=24)
>>> physical = marker_density(orbits, dims="eta1", bins=24, weight="weight")
>>> (physical / sampling).isel(t=-1).plasma.plot.lineout()
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L2913-L3021)

## plasma_plots.analysis.norm

*function*

```python
def norm(data: xr.DataArray, *, dims=None, squared: bool = False) -> xr.DataArray
```

L2 norm over ``dims`` (default: every dimension except ``t``), as a function of the rest.

The discrete norm ``√(Σ f²)``, a plain sum over the grid points: neither divided by their
number nor weighted by the volume element (see [`error()`][plasma_plots.analysis.error] and [`volume_integral()`][plasma_plots.analysis.volume_integral]
for those).

**Parameters**

- `data` (`xarray.DataArray`) — The field.
- `dims` (`str or sequence of str`) (default: `None`) — The dimensions summed over. Default: every dimension except ``t``.
- `squared` (`bool`) (default: `False`) — Return the squared norm ``Σ f²`` instead. Default: ``False``.

**Returns**

- (`xarray.DataArray`) — The norm as a function of the remaining dimensions, labeled ``"norm of ..."`` (or ``"squared norm of ..."``), keeping the run provenance attributes.

**Examples**

```pycon
>>> div_B.plasma.analysis.norm().plasma.plot.timeseries()
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L378-L415)

## plasma_plots.analysis.orbit_invariants

*function*

```python
def orbit_invariants(orbits: xr.Dataset, *, absB=None) -> xr.Dataset
```

Kinetic invariants of saved marker orbits, over ``(t, marker)``.

Computes whatever the saved quantities allow:

* ``speed`` ``|v|`` from ``v1``, ``v2``, ``v3`` (full orbits);
* with ``v_par``, ``mu``, the positions ``x``, ``y``, ``z`` and ``absB``: the guiding-centre
  ``energy`` ``v_par² / 2 + mu |B|`` and the ``pitch`` ``v_par / v``, with
  ``v = √(2 energy)``.

Their drift, e.g. [`relative_error()`][plasma_plots.analysis.relative_error] of the energy, measures the pusher's accuracy.

**Parameters**

- `orbits` (`xarray.Dataset`) — The orbits product, with variables over ``(t, marker)``.
- `absB` (`callable`) (default: `None`) — ``|B|(x, y, z)`` as a function of the physical positions (numpy arrays), e.g. ``lambda x, y, z: out.equil.absB0(*out.domain.inverse_map(x, y, z))``. Needed for the energy and the pitch.

**Returns**

- (`xarray.Dataset`) — The invariants that could be computed (``speed``, ``energy``, ``pitch``), over ``(t, marker)``. Samples where a marker has left the domain are NaN.

**Raises**

- `ValueError` — If no invariant can be computed (need ``v1``..``v3``, or ``v_par``, ``mu`` and ``absB``).

> **See Also**
>
> [`classify_orbits()`][plasma_plots.analysis.classify_orbits], [`bounce_period()`][plasma_plots.analysis.bounce_period]

**Examples**

```pycon
>>> invariants = orbit_invariants(orbits, absB=absB)
>>> relative_error(invariants.energy.isel(marker=0)).plasma.plot.timeseries()
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L2156-L2224)

## plasma_plots.analysis.oscillation_frequency

*function*

```python
def oscillation_frequency(data: xr.DataArray, *, window: tuple[float | None, float | None] = (None, None), method: str = 'zero_crossings', detrend: bool = True) -> OscillationFit | None
```

Measure the frequency of an oscillating time series from its zero crossings or its peaks.

Zero crossings are unaffected by damping or growth, and are found to a fraction of a sample by
linear interpolation; peaks are refined by a parabola through each peak and its neighbours.
The period is the slope of a straight line through the crossing (or peak) times against their
number, so a single late or early crossing hardly matters. A spectrum
([`spectral_peaks()`][plasma_plots.spectral.spectral_peaks]) resolves several frequencies at once; this
measures one, from a few periods, better than a frequency bin.

**Parameters**

- `data` (`xarray.DataArray`) — The time series, with ``t`` as its only dimension, e.g. a probe or a mode amplitude.
- `window` (`(float or None, float or None)`) (default: `(None, None)`) — The time interval ``(t0, t1)`` used; ``None`` for an open end. Default: every sample.
- `method` (`(zero_crossings, peaks)`) (default: `"zero_crossings"`) — Count the crossings of the mean (half a period apart), or the maxima (a period apart; for a signal that does not cross its mean, e.g. an energy, whose peaks are half the field's period apart). Default: ``"zero_crossings"``.
- `detrend` (`bool`) (default: `True`) — Subtract the mean over the window first, so that crossings are of the mean. Default: ``True``.

**Returns**

- (`OscillationFit or None`) — The frequency, period and the times used; ``None`` with fewer than two crossings or peaks.

**Raises**

- `ValueError` — If ``data`` has dimensions other than ``t``, or ``method`` is unknown.

> **See Also**
>
> [`damping_rate()`][plasma_plots.analysis.damping_rate] : The decay of the same oscillation.

**Examples**

```pycon
>>> oscillation_frequency(
...     phi.isel(eta1=8, eta2=0, eta3=0), window=(5.0, 40.0)
... ).omega
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L281-L375)

## plasma_plots.analysis.polar_coordinates

*function*

```python
def polar_coordinates(data: xr.DataArray, *, center=(0.0, 0.0)) -> xr.DataArray
```

Attach the polar coordinates ``r`` and ``theta`` of each point in the ``X``-``Y`` plane.

E.g. for profiles against the radius.

**Parameters**

- `data` (`xarray.DataArray`) — The array, with the ``X`` and ``Y`` coordinates.
- `center` (`(float, float)`) (default: `(0.0, 0.0)`) — The origin ``(X, Y)`` of the polar coordinates. Default: ``(0.0, 0.0)``.

**Returns**

- (`xarray.DataArray`) — ``data`` with the extra coordinates ``r`` and ``theta`` (radians, in ``(−π, π]``) of its points about ``center``.

**Examples**

```pycon
>>> n_polar = polar_coordinates(n, center=(0.0, 0.0)).isel(t=-1, eta3=0)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L2127-L2153)

## plasma_plots.analysis.power_spectrum

*function*

```python
def power_spectrum(data: xr.DataArray, *, dim: str | None = None, detrend: bool = True) -> xr.DataArray
```

The 2-D power spectrum of a ``(t, dim)`` signal, as a function of frequency and wavenumber.

A plain space-time FFT, as a function of angular frequency and wavenumber -- the basis of a
dispersion-relation plot ([`dispersion()`][plasma_plots.accessors.ArrayPlots.dispersion]),
independent of Struphy. Built on [`plasma_plots.spectral.fft()`][plasma_plots.spectral.fft], so it shares its
conventions: coefficients are divided by the sample counts, and the power sums to the mean
square of the signal.

**Parameters**

- `data` (`xarray.DataArray`) — The signal, with exactly the two dimensions ``t`` and ``dim``, each on a uniform grid.
- `dim` (`str`) (default: `None`) — The spatial dimension. Default: the sole dimension other than ``t``.
- `detrend` (`bool`) (default: `True`) — Remove the time-mean at each point of ``dim`` first, which otherwise dominates the spectrum as a spurious zero-frequency line. Default: ``True``.

**Returns**

- (`xarray.DataArray`) — The power ``|ĉ|²``, named ``power``, over ``(omega, k)``.

**Raises**

- `ValueError` — If ``dim`` is not given and ``data`` has more than one dimension besides ``t``, or
``data`` has dimensions other than ``t`` and ``dim``.

> **See Also**
>
> [`fit_dispersion_branches()`][plasma_plots.analysis.fit_dispersion_branches]

**Examples**

```pycon
>>> spectrum = power_spectrum(phi.isel(eta2=0, eta3=0), dim="eta1")
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L834-L899)

## plasma_plots.analysis.project_mode

*function*

```python
def project_mode(data: xr.DataArray, *, dim: str, number: float, kind: str = 'sin', period: float = 1.0, bin_correction: bool = False) -> xr.DataArray
```

The amplitude of one Fourier mode along a periodic dimension, as a function of the rest.

The amplitude is twice the mean over the samples of ``data`` times the basis function, e.g.
``2 ⟨f sin(…)⟩``.

**Parameters**

- `data` (`xarray.DataArray`) — The field. ``dim`` must sample one full period uniformly (a duplicate endpoint is dropped).
- `dim` (`str`) — The periodic dimension, e.g. ``"eta1"``.
- `number` (`float`) — The mode number, not a wavenumber: ``k = 2π number / L``.
- `kind` (`(sin, cos, complex)`) (default: `"sin"`) — ``"sin"`` (default) gives ``a`` of ``a sin(2π number x / period)``: ``2 ⟨f sin(…)⟩``; ``"cos"`` the cosine amplitude, and ``"complex"`` the complex amplitude ``2 ⟨f exp(−i…)⟩`` (``abs()`` the amplitude, ``np.angle()`` the phase).
- `period` (`float`) (default: `1.0`) — The period of ``dim``. Default: 1, the logical unit interval.
- `bin_correction` (`bool`) (default: `False`) — Undo the damping of binned particle data, whose bins average the mode over their width ``h``: divide by ``sinc(number h / period)``. Default: ``False``.

**Returns**

- (`xarray.DataArray`) — The amplitude as a function of the other dimensions (complex for ``"complex"``), named e.g. ``sin_1``.

**Raises**

- `ValueError` — If ``kind`` is unknown, or ``dim`` does not sample one full period uniformly.

**Examples**

```pycon
>>> amplitude = project_mode(e1, dim="eta1", number=1)  # sine amplitude over t
>>> phase = np.angle(project_mode(rho, dim="eta2", number=3, kind="complex"))
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L1699-L1775)

## plasma_plots.analysis.quadrature_weights

*function*

```python
def quadrature_weights(coordinate, *, period: float | None = None) -> np.ndarray
```

Quadrature weights for samples of a logical coordinate in ``[0, 1]``, or of an angle.

Struphy's cell centers (uniform, half a cell from each end) get the midpoint rule, which
integrates over the whole unit interval; a uniform grid over one full ``period`` (an angle,
e.g. GVEC's) the rectangle rule, exact for its Fourier modes; any other grid the trapezoidal
rule over the sampled range. A single point (a 2-D run's flat direction) has weight 1.

**Parameters**

- `coordinate` (`array_like of float`) — The sample points along one logical direction.
- `period` (`float`) (default: `None`) — The direction's period, if it is an angle (see [`plasma_plots.arrays.angle_period()`][plasma_plots.arrays.angle_period]). Default: none.

**Returns**

- (`numpy.ndarray`) — One weight per sample point.

> **See Also**
>
> [`volume_integral()`][plasma_plots.analysis.volume_integral], [`field_energy()`][plasma_plots.analysis.field_energy]

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L907-L951)

## plasma_plots.analysis.quasisymmetry_error

*function*

```python
def quasisymmetry_error(data: xr.DataArray, *, helicity='QA', angles: str = 'boozer') -> xr.DataArray
```

The quasi-symmetry error of ``|B|`` on each flux surface, from its Boozer spectrum.

A field is quasi-symmetric with helicity ``(M, N)`` when ``|B|`` depends on the Boozer
angles only through ``M θ_B − N ζ_B``: its harmonics ``(m, n)`` (in the convention of
[`boozer_spectrum()`][plasma_plots.analysis.boozer_spectrum], ``cos(m θ_B + n ζ_B)``) all satisfy ``m N + n M = 0``. The error is
the symmetry-breaking content relative to the mean field,

    f_QS(ρ) = √(Σ_breaking B_mn²) / B_00,

the common measure of how far an equilibrium is from quasi-symmetry (see e.g. DESC's and
SIMSOPT's quasi-symmetry objectives). ``"QA"`` is quasi-axisymmetry (``N = 0``: only
``n = 0`` survives), ``"QP"`` quasi-poloidal symmetry (``M = 0``), ``"QH"`` quasi-helical
symmetry with ``(1, nfp)``, tried with both signs of ``N`` (the handedness), keeping the
smaller error.

**Parameters**

- `data` (`xarray.DataArray`) — ``|B|`` on a Boozer grid, as for [`boozer_spectrum()`][plasma_plots.analysis.boozer_spectrum].
- `helicity` (`(QA, QP, QH)`) (default: `"QA"`) — The symmetry: a name, or ``(M, N)`` with ``N`` a full-torus toroidal mode number (its sign picks the handedness). Default: ``"QA"``.
- `angles` (`(boozer, any)`) (default: `"boozer"`) — As for [`boozer_spectrum()`][plasma_plots.analysis.boozer_spectrum].

**Returns**

- (`xarray.DataArray`) — ``f_QS`` over the radius (and every other non-angle dimension), named ``quasisymmetry_error``, with the ``helicity`` used in its attrs.

**Raises**

- `ValueError` — As for [`boozer_spectrum()`][plasma_plots.analysis.boozer_spectrum], or for an unknown helicity.

> **See Also**
>
> [`boozer_spectrum()`][plasma_plots.analysis.boozer_spectrum] : The harmonics themselves.

> **References**
>
> M. Landreman and E. Paul, "Magnetic fields with precise quasisymmetry for plasma
> confinement", Phys. Rev. Lett. 128, 035001 (2022).

**Examples**

```pycon
>>> f_qs = quasisymmetry_error(boozer.mod_B, helicity="QA")
>>> f_qs.plasma.plot.lineout()
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L2375-L2461)

## plasma_plots.analysis.rational_surfaces

*function*

```python
def rational_surfaces(profile: xr.DataArray, *, count: int = 4, nfp: int | None = None, max_denominator: int = 12) -> xr.DataArray
```

Where a rotational transform (or safety factor) profile takes low-order rational values.

The values are the fractions ``n/m`` with ``m ≤ max_denominator`` and ``n`` a multiple of
``nfp`` (for ι = n/m, the resonances of a stellarator with nfp field periods) that the profile
reaches; the ``count`` of lowest order (smallest ``m``, then ``|n|``) are kept. Each crossing
is found by linear interpolation between samples, so a non-monotonic profile gives several.

**Parameters**

- `profile` (`xarray.DataArray`) — A 1-D profile, e.g. ``iota`` over ``rho``.
- `count` (`int`) (default: `4`) — The number of rational values. Default: 4.
- `nfp` (`int`) (default: `None`) — The numerators are multiples of it. Default: the profile's ``nfp`` attribute (GVEC's, see [`plasma_plots.gvec.from_gvec()`][plasma_plots.gvec.from_gvec]), else 1.
- `max_denominator` (`int`) (default: `12`) — The largest ``m``. Default: 12.

**Returns**

- (`xarray.DataArray`) — The positions of the crossings (in the profile's coordinate), over a ``surface`` dimension with the coordinates ``n``, ``m`` and ``value`` (``n/m``), lowest order first; empty if the profile reaches none.

**Raises**

- `ValueError` — If ``profile`` is not one-dimensional.

> **See Also**
>
> [`plasma_plots.plotting.plot_lineout()`][plasma_plots.plotting.plot_lineout] : ``rationals=`` marks them on the profile.

**Examples**

```pycon
>>> rational_surfaces(ev.iota, count=3)
>>> ev.iota.plasma.plot.lineout(rationals=3)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L1175-L1272)

## plasma_plots.analysis.reconnected_flux

*function*

```python
def reconnected_flux(data: xr.DataArray, *, relative: bool = True, o_point=None, x_point=None) -> xr.DataArray
```

The reconnected flux over time: the flux function between an O-point and an X-point.

At each time the O- and X-points are found ([`critical_points()`][plasma_plots.analysis.critical_points]); the flux of each
O-point's island is ``|A(O) − A(X)|`` to the X-point nearest in flux (the one on its
separatrix), and the island with the most flux is taken, the dominant island or
reconnection region. ``o_point`` and ``x_point`` pin the pair to the points nearest given
positions instead (e.g. the X-point of a Harris sheet at the origin). For a tearing mode the flux grows as ``exp(γt)``
(see [`growth_rate()`][plasma_plots.analysis.growth_rate]); for the GEM challenge it is the usual reconnected-flux curve.

**Parameters**

- `data` (`xarray.DataArray`) — The flux function over ``t`` and two logical directions, e.g. [`flux_function()`][plasma_plots.analysis.flux_function] of ``B``.
- `relative` (`bool`) (default: `True`) — Subtract the value at the first time. Default: True.
- `o_point` (`(float, float)`) (default: `None`) — The logical coordinates near which to look for the O-point. Default: the dominant island's.
- `x_point` (`(float, float)`) (default: `None`) — The same for the X-point.

**Returns**

- (`xarray.DataArray`) — ``reconnected_flux`` over ``t``, labeled ``ΔΨ``, NaN at times without such a pair.

**Raises**

- `ValueError` — If ``data`` has no ``t`` dimension, or is not over two logical directions.

> **See Also**
>
> [`critical_points()`][plasma_plots.analysis.critical_points] : The points themselves.

**Examples**

```pycon
>>> flux = reconnected_flux(flux_function(B))
>>> flux.plasma.plot.timeseries(fit=(10.0, 30.0))
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L2735-L2822)

## plasma_plots.analysis.relative_error

*function*

```python
def relative_error(data: xr.DataArray, *, ref=None, skip_first=True) -> xr.DataArray
```

Absolute relative deviation from an explicit reference or first sample.

**Parameters**

- `data` (`xarray.DataArray`) — The array, with a ``t`` dimension, e.g. a conserved energy.
- `ref` (`float or array_like or xarray.DataArray`) (default: `None`) — The reference, broadcast against ``data``; must be non-zero everywhere. Default: ``data`` at the first time sample.
- `skip_first` (`bool`) (default: `True`) — Leave out the first time sample (zero against the default reference). Default: ``True``.

**Returns**

- (`xarray.DataArray`) — ``|data − ref| / |ref|``, labeled ``"relative error of ..."`` with empty units.

**Raises**

- `ValueError` — If the reference is zero anywhere.

> **See Also**
>
> [`drift()`][plasma_plots.analysis.drift] : The signed deviation.

**Examples**

```pycon
>>> relative_error(out.scalars["en_tot"]).plasma.plot.timeseries()
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L584-L624)

## plasma_plots.analysis.spatial_average

*function*

```python
def spatial_average(data: xr.DataArray, *, dims=None) -> xr.DataArray
```

Mean over the logical space dimensions, e.g. a binned f(t, eta1, v1) becomes f(t, v1).

The mean is uniform in the logical coordinates, which is the volume average on a Cartesian
domain; on a mapped domain it is not weighted by the Jacobian (use [`volume_integral()`][plasma_plots.analysis.volume_integral]
for that). Physical ``X``, ``Y``, ``Z`` coordinates that depend on the averaged dimensions
are dropped.

**Parameters**

- `data` (`xarray.DataArray`) — The array, e.g. a binned distribution function.
- `dims` (`str or sequence of str`) (default: `None`) — The dimensions averaged over. Default: every logical dimension (``eta1``, ``eta2``, ``eta3``, or GVEC's ``rho``, ``theta``, ``zeta``) that ``data`` has.

**Returns**

- (`xarray.DataArray`) — The mean over ``dims``, as a function of the remaining dimensions, with the attributes of ``data`` and the label ``"average of ..."``.

**Raises**

- `ValueError` — If ``data`` has none of the default dimensions, or lacks one of ``dims``.

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L475-L507)

## plasma_plots.analysis.surface_average

*function*

```python
def surface_average(data: xr.DataArray, *, jacobian=None, domain=None, quadrature=None) -> xr.DataArray
```

The flux-surface average ``⟨f⟩ = ∫ f √g dθ dζ / ∫ √g dθ dζ`` over the two angles.

The angles are the second and third logical dimensions (``eta2``, ``eta3``, or GVEC's
``theta``, ``zeta`` or Boozer or PEST angles; see [`plasma_plots.arrays.logical_dims()`][plasma_plots.arrays.logical_dims]),
and the average is a function of the radius and every other dimension. The angles are
integrated as in [`volume_integral()`][plasma_plots.analysis.volume_integral]: GVEC's Gauss weights where given, else the
rectangle rule over a full period. On the magnetic axis, where ``√g`` vanishes, it is the
plain mean over the angles.

**Parameters**

- `data` (`xarray.DataArray`) — The field, over the radial and both angular logical dimensions.
- `jacobian` (`xarray.DataArray`) (default: `None`) — ``√g`` on the same grid, e.g. GVEC's ``Jac`` (its absolute value is used). Default: from ``domain``, else the numerical Jacobian of the ``X``, ``Y``, ``Z`` coordinates, which needs at least two radial points.
- `domain` (`struphy domain`) (default: `None`) — The mapping, for the exact ``√g`` of a Struphy run.
- `quadrature` (`dict`) (default: `None`) — Explicit weights for the angles, as for [`volume_integral()`][plasma_plots.analysis.volume_integral].

**Returns**

- (`xarray.DataArray`) — ``⟨f⟩`` as a function of the radius (and every non-spatial dimension), labeled ``"⟨...⟩"``.

**Raises**

- `ValueError` — If an angle is missing, or ``√g`` can't be computed (no ``jacobian``, ``domain`` or
``X``, ``Y``, ``Z`` on at least two radial points).

> **See Also**
>
> [`volume_integral()`][plasma_plots.analysis.volume_integral]

**Examples**

```pycon
>>> # GVEC: ⟨|B|⟩(rho)
>>> ev.mod_B.plasma.analysis.surface_average(jacobian=ev.Jac)
>>> surface_average(p, domain=out.domain)  # Struphy
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L1095-L1172)

## plasma_plots.analysis.toroidal_components

*function*

```python
def toroidal_components(vector: xr.DataArray, *, R0: float, Z0: float = 0.0) -> xr.DataArray
```

Cartesian components rotated to the local ``(radial, poloidal, toroidal)`` directions.

The directions are those about a circular magnetic axis at major radius ``R0`` (height
``Z0``), at every point. The radial direction points away from the axis in the poloidal
plane, the poloidal one along ``θ = atan2(Z − Z0, R − R0)``, the toroidal one along ``φ``:
the natural components for waves in a tokamak or torus.

**Parameters**

- `vector` (`xarray.DataArray`) — The field in Cartesian components, with a ``component`` dimension of size 3 and the ``X``, ``Y``, ``Z`` coordinates.
- `R0` (`float`) — The major radius of the magnetic axis.
- `Z0` (`float`) (default: `0.0`) — The height of the magnetic axis. Default: 0.

**Returns**

- (`xarray.DataArray`) — The components, with ``component`` coordinates ``"radial"``, ``"poloidal"``, ``"toroidal"`` and the dimensions of ``vector``.

**Raises**

- `ValueError` — If ``vector`` has no ``component`` dimension of size 3.

> **See Also**
>
> [`cylindrical_components()`][plasma_plots.analysis.cylindrical_components]

**Examples**

```pycon
>>> local = toroidal_components(u, R0=3.0)
>>> local.sel(component="poloidal").isel(t=-1, eta3=0).plasma.plot.slice()
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L2060-L2124)

## plasma_plots.analysis.velocity_moments

*function*

```python
def velocity_moments(f: xr.DataArray, *, dims=None) -> xr.Dataset
```

Moments of a binned distribution function over its velocity dimensions.

The moments are functions of the remaining dimensions, for example ``(t, eta1)`` for an
``e1_v1`` product. The integrals are sums over the bins, weighted by the bin widths (from
``numpy.gradient`` of the velocity coordinates). The values keep the normalization of the
run; see ``Output.to_si``.

**Parameters**

- `f` (`xarray.DataArray`) — The binned distribution function. A product named ``delta_f`` has only the density, which is then the density perturbation, because its mean and variance are not defined.
- `dims` (`str or sequence of str`) (default: `None`) — The velocity dimensions integrated over, each with a coordinate of at least two bins. Default: every one of ``v1``, ``v2``, ``v3`` that ``f`` has.

**Returns**

- (`xarray.Dataset`) — The moments, over the remaining dimensions:  * ``density``: the zeroth moment, ``n = ∫ f dv``. * ``mean_<dim>``: the mean velocity ``u = ∫ v f dv / n`` along each dimension. * ``variance_<dim>``: ``∫ (v − u)² f dv / n``. In normalized units this is the   temperature over the particle mass along that direction, ``T/m``.  Where the density is not positive, the mean and variance are NaN.

**Raises**

- `ValueError` — If ``f`` has none of the default dimensions, lacks one of ``dims``, or a dimension has
fewer than two bins.

**Examples**

```pycon
>>> moments = velocity_moments(f)  # f(t, eta1, v1)
>>> moments.variance_v1.isel(t=-1).plasma.plot.lineout()
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L519-L581)

## plasma_plots.analysis.volume_integral

*function*

```python
def volume_integral(data: xr.DataArray, *, form: int = 0, weight=None, domain=None, quadrature=None, jacobian=None) -> xr.DataArray
```

``∫ w f dV`` over the logical grid, as a function of every other dimension (e.g. ``t``).

A function (``form=0``) is integrated with the volume element ``|√g| dη``; a density (a
3-form, e.g. Struphy's L2 fields, ``form=3``) as ``∫ f dη``.

**Parameters**

- `data` (`xarray.DataArray`) — The integrand, with the three logical dimensions (``eta1``, ``eta2``, ``eta3``, or GVEC's ``rho``, ``theta``, ``zeta``; see [`plasma_plots.arrays.logical_dims()`][plasma_plots.arrays.logical_dims]).
- `form` (`(0, 3)`) (default: `0`) — ``0`` (default) for a function, ``3`` for a density.
- `weight` (`array_like`) (default: `None`) — An extra weight ``w`` over the logical grid, broadcast against ``data``.
- `domain` (`struphy domain`) (default: `None`) — The mapping (``out.domain``), which gives the exact ``|√g|``. Without one, ``|√g|`` comes from the numerical Jacobian of the ``X``, ``Y``, ``Z`` coordinates (see [`plasma_plots.arrays.mapping_jacobian()`][plasma_plots.arrays.mapping_jacobian]). Only needed for ``form=0``.
- `quadrature` (`dict`) (default: `None`) — Maps logical dimensions to explicit weights, one per point (e.g. Gauss weights, see ``out.analysis.quadrature_grid()``). Directions left out take a ``<dim>_weight`` coordinate (GVEC's integration points, see [`plasma_plots.gvec.from_gvec()`][plasma_plots.gvec.from_gvec]), else [`quadrature_weights()`][plasma_plots.analysis.quadrature_weights]: the midpoint rule on Struphy's cell centers, the rectangle rule over a full period of an angle, else trapezoidal.
- `jacobian` (`xarray.DataArray`) (default: `None`) — The Jacobian determinant ``√g`` on the grid, e.g. GVEC's ``Jac`` (its absolute value is used), instead of ``domain`` or the numerical one. Only needed for ``form=0``.

**Returns**

- (`xarray.DataArray`) — The integral, as a function of the dimensions other than the logical ones, labeled ``"integral of ..."``. On GVEC's grid over one field period, multiply by ``nfp`` for the whole device.

**Raises**

- `ValueError` — If ``form`` is not 0 or 3, a logical dimension is missing, or a quadrature has the wrong
number of weights.

> **See Also**
>
> [`field_energy()`][plasma_plots.analysis.field_energy], [`quadrature_weights()`][plasma_plots.analysis.quadrature_weights]

**Examples**

```pycon
>>> total_charge = volume_integral(rho, domain=out.domain)
>>> mass = volume_integral(n3, form=3)  # a 3-form: no |√g|
>>> # GVEC
>>> volume = volume_integral(xr.ones_like(ev.mod_B), jacobian=ev.Jac) * ev.nfp
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L1016-L1092)

## plasma_plots.analysis.weight_statistics

*function*

```python
def weight_statistics(markers: xr.Dataset, *, weight: str = 'weight') -> xr.Dataset
```

Statistics of the marker weights over time, with the noise a δf (or PIC) estimate carries.

Besides the mean, spread and extremes of the weights, the relative statistical error of
the total (the zeroth moment, ``Σw``) that random marker positions would give,
``noise = √(Σw²) / |Σw|``, and Kish's effective number of markers ``(Σw)² / Σw²`` (the
number of equal-weight markers with the same noise). Markers that have left the domain
(every quantity zero, as Struphy stores them) are left out.

**Parameters**

- `markers` (`xarray.Dataset`) — A marker Dataset with the weights over ``(t, marker)`` (or ``marker``), e.g. an orbits product.
- `weight` (`str`) (default: `'weight'`) — The weight variable. Default: ``"weight"``, Struphy's.

**Returns**

- (`xarray.Dataset`) — Over ``t`` (scalars without ``t``): ``mean``, ``std``, ``min``, ``max``, ``total``, ``noise``, ``effective_markers`` and ``count`` (the markers in the domain).

**Raises**

- `ValueError` — If ``markers`` has no variable ``weight``.

> **See Also**
>
> [`marker_density()`][plasma_plots.analysis.marker_density] : Where the markers are, against what they represent.
> [`plasma_plots.plotting.plot_weight_histogram()`][plasma_plots.plotting.plot_weight_histogram] : The distribution of the weights.

**Examples**

```pycon
>>> stats = weight_statistics(orbits)
>>> stats.noise.plasma.plot.timeseries(logy=False)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/analysis.py#L2845-L2910)