plasma_plots.analysis

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(...).

Attributes

Name	Description
HELICITIES	No description.
ORBIT_CLASSES	No description.
VELOCITY_DIMS	No description.

Classes

Name	Description
BranchFit	A dispersion branch fitted as omega = velocity k, from fitdispersion_branches().
ConvergenceFit	A power law error = constant size*order, from convergence_order().
FitResult	An exponential fit exp(rate t + intercept), from growth_rate() or damping_rate().
GrowthFit	Configuration for an exponential growth-rate fit.
OscillationFit	The frequency of an oscillating time series, from oscillation_frequency().

Functions

Name	Description
boozer_spectrum	The Boozer harmonics B_mn of a quantity on flux surfaces, real amplitudes over the radius.
bounce_period	The bounce period of each trapped marker.
classify_orbits	Classify each marker of an orbits product as passing (0), trapped (1) or lost (-1).
convergence_order	Fit error = constant size*order in log-log space.
critical_points	The O-points (extrema) and X-points (saddles) of a 2-D flux function.
curl	The curl of a vector field on a mapped domain, in Cartesian components (x, y, z).
cylindrical_components	Cartesian components rotated to cylindrical ones (R, phi, Z) about the Z axis.
damping_rate	Fit exp(rate*t + intercept) to the envelope of an oscillating time series.
divergence	The divergence sum_a d v_a / d x_a of a vector field on a mapped domain.
drift	Signed deviation from an explicit reference or the first time sample.
envelope	Local maxima of a time series: the interior samples not smaller than their neighbours.
error	The error of data against an exact solution, as a function of the other dimensions.
evaluate_on	function evaluated on the coordinates of data, broadcast to its shape.
field_energy	The quadratic energy α ½ ∫ w ωᵀ A ω dη, as a function of time.
fit_dispersion_branches	Fit omega = v k to each of nbranches straight ridges in an (omega, k) power spectrum.
flux_function	The flux function (or stream function) A of a 2-D, divergence-free in-plane field.
gradient	The Cartesian gradient of a scalar field on the logical grid.
growth_rate	Fit exp(rate*t + intercept) to a time series, using only finite, positive samples.
loss_map	Each marker's initial phase-space position, whether it is lost, and when.
lost_fraction	The fraction of markers that have left the domain, over time.
marker_density	Bin the markers over position variables: the number (or weight) of markers per unit volume.
norm	L2 norm over dims (default: every dimension except t), as a function of the rest.
orbit_invariants	Kinetic invariants of saved marker orbits, over (t, marker).
oscillation_frequency	Measure the frequency of an oscillating time series from its zero crossings or its peaks.
polar_coordinates	Attach the polar coordinates r and theta of each point in the X-Y plane.
power_spectrum	The 2-D power spectrum of a (t, dim) signal, as a function of frequency and wavenumber.
project_mode	The amplitude of one Fourier mode along a periodic dimension, as a function of the rest.
quadrature_weights	Quadrature weights for samples of a logical coordinate in [0, 1], or of an angle.
quasisymmetry_error	The quasi-symmetry error of |B| on each flux surface, from its Boozer spectrum.
rational_surfaces	Where a rotational transform (or safety factor) profile takes low-order rational values.
reconnected_flux	The reconnected flux over time: the flux function between an O-point and an X-point.
relative_error	Absolute relative deviation from an explicit reference or first sample.
spatial_average	Mean over the logical space dimensions, e.g. a binned f(t, eta1, v1) becomes f(t, v1).
surface_average	The flux-surface average ⟨f⟩ = ∫ f √g dθ dζ / ∫ √g dθ dζ over the two angles.
toroidal_components	Cartesian components rotated to the local (radial, poloidal, toroidal) directions.
velocity_moments	Moments of a binned distribution function over its velocity dimensions.
volume_integral	∫ w f dV over the logical grid, as a function of every other dimension (e.g. t).
weight_statistics	Statistics of the marker weights over time, with the noise a δf (or PIC) estimate carries.

HELICITIESattributemodule attribute#

View source

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

ORBIT_CLASSESattributemodule attribute#

View source

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

VELOCITY_DIMSattributemodule attribute#

View source

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

BranchFitclassdataclass#

View source

class BranchFit(velocity: float, k: np.ndarray, omega: np.ndarray)

A dispersion branch fitted as omega = velocity * k, from fit_dispersion_branches().

Attributes

Name	Type	Description
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.

kattributeinstance attribute#

View source

k: np.ndarray

omegaattributeinstance attribute#

View source

omega: np.ndarray

velocityattributeinstance attribute#

View source

velocity: float

ConvergenceFitclassdataclass#

View source

class ConvergenceFit(order: float, constant: float, sizes: np.ndarray, fitted: np.ndarray)

A power law error = constant * size**order, from convergence_order().

Attributes

Name	Type	Description
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.

constantattributeinstance attribute#

View source

constant: float

fittedattributeinstance attribute#

View source

fitted: np.ndarray

orderattributeinstance attribute#

View source

order: float

sizesattributeinstance attribute#

View source

sizes: np.ndarray

FitResultclassdataclass#

View source

class FitResult(rate: float, intercept: float, time: np.ndarray, fitted: np.ndarray)

An exponential fit exp(rate t + intercept), from growth_rate() or damping_rate().

Attributes

Name	Type	Description
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.

fittedattributeinstance attribute#

View source

fitted: np.ndarray

interceptattributeinstance attribute#

View source

intercept: float

rateattributeinstance attribute#

View source

rate: float

timeattributeinstance attribute#

View source

time: np.ndarray

GrowthFitclassdataclass#

View source

class GrowthFit(window: tuple[float | None, float | None] = (None, None), amplitude_from_quadratic: bool = False)

Configuration for an exponential growth-rate fit.

Passed to growth_rate() and damping_rate().

Attributes

Name	Type	Description
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.

amplitude_from_quadraticattributeclass attributeinstance attribute#

View source

amplitude_from_quadratic: bool = False

windowattributeclass attributeinstance attribute#

View source

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

OscillationFitclassdataclass#

View source

class OscillationFit(omega: float, period: float, times: np.ndarray, method: str)

The frequency of an oscillating time series, from oscillation_frequency().

Attributes

Name	Type	Description
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".

methodattributeinstance attribute#

View source

method: str

omegaattributeinstance attribute#

View source

omega: float

periodattributeinstance attribute#

View source

period: float

timesattributeinstance attribute#

View source

times: np.ndarray

boozer_spectrumfunction#

View source

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() over the Boozer angles, combined into the real amplitudes of 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()).

Parameters

Name	Type	Default	Description
data	xarray.DataArray	required	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	None	Keep only the top harmonics with the largest peak amplitude over the radius, strongest first. Default: all.
angles	(boozer, any)	"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() : The symmetry-breaking part of the spectrum. plasma_plots.spectral_plots.plot_boozer_spectrum() : The harmonics over the radius.

Examples

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

bounce_periodfunction#

View source

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

Name	Type	Default	Description
orbits	xarray.Dataset	required	The orbits product, with variables over (t, marker).
v_par	str	'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()

Examples

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

classify_orbitsfunction#

View source

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

Name	Type	Default	Description
orbits	xarray.Dataset	required	The orbits product, with variables over (t, marker).
v_par	str	'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.

Raises

ValueError

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

See Also

bounce_period(), orbit_invariants()

Examples

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

convergence_orderfunction#

View source

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

Name	Type	Description
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().

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

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

critical_pointsfunction#

View source

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

Name	Type	Default	Description
data	xarray.DataArray	required	The flux function, e.g. 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()); a periodic box without either is treated as bounded, and loses the points on its edges.
refine	bool	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() : The flux between an O- and an X-point over time. plasma_plots.plotting.plot_critical_points() : The points over the flux contours.

Examples

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

curlfunction#

View source

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(). 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

Name	Type	Default	Description
vector	xarray.DataArray	required	The vector field, with a component dimension of size 3 and the dimensions eta1, eta2, eta3.
components	(cartesian, contravariant)	"cartesian"	What the components are: "cartesian" (default) (x, y, z), or "contravariant" components, which are pushed forward first.
domain	struphy domain	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().

See Also

gradient(), divergence()

Examples

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

cylindrical_componentsfunction#

View source

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

Name	Type	Description
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()

Examples

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

damping_ratefunction#

View source

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() on the raw series would fit the oscillation. The envelope is the series’ local maxima, see envelope().

Parameters

Name	Type	Default	Description
data	xarray.DataArray	required	The oscillating time series, with t as its only dimension.
fit	GrowthFit	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(), envelope()

Examples

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

divergencefunction#

View source

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(): spectral around periodic angles; a flat direction (2-D run) is left out. E.g. a div B check of an MHD run.

Parameters

Name	Type	Default	Description
vector	xarray.DataArray	required	The vector field, with a component dimension of size 3 and the dimensions eta1, eta2, eta3.
components	(cartesian, contravariant)	"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	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().

See Also

gradient(), curl()

Examples

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

driftfunction#

View source

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

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

Parameters

Name	Type	Default	Description
data	xarray.DataArray	required	The array, with a t dimension.
ref	float or array_like or xarray.DataArray	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() : The absolute relative deviation.

Examples

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

envelopefunction#

View source

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

Name	Type	Description
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.

errorfunction#

View source

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()), 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

Name	Type	Default	Description
data	xarray.DataArray	required	The numerical solution.
exact	(callable, array_like, number or xarray.DataArray)	required	The exact solution: an array (aligned with data, or broadcast to it) or a function of its coordinates, evaluated with evaluate_on(), e.g. lambda x, y, z, t: ....
norm	(rms, max, l1, l2, pointwise)	"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	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	None	The dimensions the norm is taken over. Default: every dimension but t. Weighted norms need exactly eta1, eta2, eta3.
weighted	bool	False	Integrate over the physical volume instead of averaging over the grid points (no effect on "max" and "pointwise"). Default: False.
domain	struphy domain	None	The mapping (out.domain), for the exact |√g| of weighted norms. Default: from the X, Y, Z coordinates.
args	sequence of str	None	The coordinates passed to a callable exact, as in 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() : The order of a series of errors.

Examples

>>> 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)

evaluate_onfunction#

View source

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

Name	Type	Default	Description
data	xarray.DataArray	required	The array whose coordinates and shape are used.
function	callable	required	Called with the coordinates named by args, positionally, as xarray.DataArrays; returns an array (or number) broadcastable to data.
args	sequence of str	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

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

field_energyfunction#

View source

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(), 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

Name	Type	Default	Description
data	xarray.DataArray	required	The field, with the dimensions eta1, eta2, eta3; vectors have a component dimension of size 3.
form	(None, 0, 1, 2, 3, v)	None	What data holds, see the table. Default: None, a function or Cartesian vector components.
weight	array_like	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	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()).
normalization	float	1.0	The prefactor α. Default: 1.
quadrature	dict	None	Explicit quadrature weights per logical dimension, as for 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()

Examples

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

fit_dispersion_branchesfunction#

View source

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(), 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

Name	Type	Default	Description
spectrum	xarray.DataArray	required	The power spectrum, with dimensions omega and k.
n_branches	int	required	The number of branches, at least 1.
k_range	(float, float)	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	0.5	Maxima count only above this fraction of the column’s peak power. Default: 0.5.
order	int	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()

Examples

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

flux_functionfunction#

View source

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

Name	Type	Description
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

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

gradientfunction#

View source

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()). 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

Name	Type	Default	Description
data	xarray.DataArray	required	The scalar field, with the dimensions eta1, eta2, eta3 (at least one with more than one point) and no component dimension.
domain	struphy domain	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(), curl()

Examples

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

growth_ratefunction#

View source

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

Name	Type	Default	Description
data	xarray.DataArray	required	The time series, with t as its only dimension (e.g. a field energy).
fit	GrowthFit	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() : The same fit to the envelope of an oscillating series.

Examples

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

loss_mapfunction#

View source

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() (absB needed for the first two).

Parameters

Name	Type	Default	Description
markers	xarray.Dataset	required	An orbits product over (t, marker).
x	str	'v_par'	The first quantity. Default: "v_par".
y	str	None	The second quantity. Default: "mu", or "v_perp" without mu.
t	int or float	0	The time of the plotted values: an integer position (default 0, the initial one) or a float nearest value.
absB	callable	None	|B|(x, y, z), for "energy" and "pitch" (see 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() : The losses over time. plasma_plots.plotting.plot_loss_map() : The plot.

Examples

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

lost_fractionfunction#

View source

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

Name	Type	Default	Description
markers	xarray.Dataset	required	A marker Dataset over (t, marker), e.g. an orbits product.
weight	str	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() : Which markers are lost, and when. classify_orbits() : Lost markers as a class.

Examples

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

marker_densityfunction#

View source

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

Name	Type	Default	Description
markers	xarray.Dataset	required	A marker Dataset with position variables over (t, marker) (or marker).
dims	str or sequence of str	('eta1')	The position variables to bin over, e.g. ("eta1",) or ("x", "y"). Default: ("eta1",).
bins	int or sequence of int	32	The number of bins, for all or for each. Default: 32.
weight	str	None	Weigh each marker by this variable, e.g. "weight". Default: count the markers.
ranges	dict	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() : The weights’ statistics. plasma_plots.plotting.plot_marker_density() : Sampling against physical density.

Examples

>>> 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()

normfunction#

View source

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() and volume_integral() for those).

Parameters

Name	Type	Default	Description
data	xarray.DataArray	required	The field.
dims	str or sequence of str	None	The dimensions summed over. Default: every dimension except t.
squared	bool	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

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

orbit_invariantsfunction#

View source

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() of the energy, measures the pusher’s accuracy.

Parameters

Name	Type	Default	Description
orbits	xarray.Dataset	required	The orbits product, with variables over (t, marker).
absB	callable	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(), bounce_period()

Examples

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

oscillation_frequencyfunction#

View source

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()) resolves several frequencies at once; this measures one, from a few periods, better than a frequency bin.

Parameters

Name	Type	Default	Description
data	xarray.DataArray	required	The time series, with t as its only dimension, e.g. a probe or a mode amplitude.
window	(float or None, float or None)	(None, None)	The time interval (t0, t1) used; None for an open end. Default: every sample.
method	(zero_crossings, peaks)	"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	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() : The decay of the same oscillation.

Examples

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

polar_coordinatesfunction#

View source

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

Name	Type	Default	Description
data	xarray.DataArray	required	The array, with the X and Y coordinates.
center	(float, float)	(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

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

power_spectrumfunction#

View source

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()), independent of Struphy. Built on 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

Name	Type	Default	Description
data	xarray.DataArray	required	The signal, with exactly the two dimensions t and dim, each on a uniform grid.
dim	str	None	The spatial dimension. Default: the sole dimension other than t.
detrend	bool	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()

Examples

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

project_modefunction#

View source

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

Name	Type	Default	Description
data	xarray.DataArray	required	The field. dim must sample one full period uniformly (a duplicate endpoint is dropped).
dim	str	required	The periodic dimension, e.g. "eta1".
number	float	required	The mode number, not a wavenumber: k = 2π number / L.
kind	(sin, cos, complex)	"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	1.0	The period of dim. Default: 1, the logical unit interval.
bin_correction	bool	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

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

quadrature_weightsfunction#

View source

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

Name	Type	Default	Description
coordinate	array_like of float	required	The sample points along one logical direction.
period	float	None	The direction’s period, if it is an angle (see plasma_plots.arrays.angle_period()). Default: none.

Returns

numpy.ndarray

One weight per sample point.

See Also

volume_integral(), field_energy()

quasisymmetry_errorfunction#

View source

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(), 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

Name	Type	Default	Description
data	xarray.DataArray	required	|B| on a Boozer grid, as for boozer_spectrum().
helicity	(QA, QP, QH)	"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)	"boozer"	As for 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(), or for an unknown helicity.

See Also

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

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

rational_surfacesfunction#

View source

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

Name	Type	Default	Description
profile	xarray.DataArray	required	A 1-D profile, e.g. iota over rho.
count	int	4	The number of rational values. Default: 4.
nfp	int	None	The numerators are multiples of it. Default: the profile’s nfp attribute (GVEC’s, see plasma_plots.gvec.from_gvec()), else 1.
max_denominator	int	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() : rationals= marks them on the profile.

Examples

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

reconnected_fluxfunction#

View source

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()); 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()); for the GEM challenge it is the usual reconnected-flux curve.

Parameters

Name	Type	Default	Description
data	xarray.DataArray	required	The flux function over t and two logical directions, e.g. flux_function() of B.
relative	bool	True	Subtract the value at the first time. Default: True.
o_point	(float, float)	None	The logical coordinates near which to look for the O-point. Default: the dominant island’s.
x_point	(float, float)	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() : The points themselves.

Examples

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

relative_errorfunction#

View source

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

Absolute relative deviation from an explicit reference or first sample.

Parameters

Name	Type	Default	Description
data	xarray.DataArray	required	The array, with a t dimension, e.g. a conserved energy.
ref	float or array_like or xarray.DataArray	None	The reference, broadcast against data; must be non-zero everywhere. Default: data at the first time sample.
skip_first	bool	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() : The signed deviation.

Examples

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

spatial_averagefunction#

View source

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() for that). Physical X, Y, Z coordinates that depend on the averaged dimensions are dropped.

Parameters

Name	Type	Default	Description
data	xarray.DataArray	required	The array, e.g. a binned distribution function.
dims	str or sequence of str	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.

surface_averagefunction#

View source

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()), and the average is a function of the radius and every other dimension. The angles are integrated as in 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

Name	Type	Default	Description
data	xarray.DataArray	required	The field, over the radial and both angular logical dimensions.
jacobian	xarray.DataArray	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	None	The mapping, for the exact √g of a Struphy run.
quadrature	dict	None	Explicit weights for the angles, as for 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()

Examples

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

toroidal_componentsfunction#

View source

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

Name	Type	Default	Description
vector	xarray.DataArray	required	The field in Cartesian components, with a component dimension of size 3 and the X, Y, Z coordinates.
R0	float	required	The major radius of the magnetic axis.
Z0	float	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()

Examples

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

velocity_momentsfunction#

View source

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

Name	Type	Default	Description
f	xarray.DataArray	required	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	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

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

volume_integralfunction#

View source

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

Name	Type	Default	Description
data	xarray.DataArray	required	The integrand, with the three logical dimensions (eta1, eta2, eta3, or GVEC’s rho, theta, zeta; see plasma_plots.arrays.logical_dims()).
form	(0, 3)	0	0 (default) for a function, 3 for a density.
weight	array_like	None	An extra weight w over the logical grid, broadcast against data.
domain	struphy domain	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()). Only needed for form=0.
quadrature	dict	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()), else 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	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(), quadrature_weights()

Examples

>>> 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

weight_statisticsfunction#

View source

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

Name	Type	Default	Description
markers	xarray.Dataset	required	A marker Dataset with the weights over (t, marker) (or marker), e.g. an orbits product.
weight	str	'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() : Where the markers are, against what they represent. plasma_plots.plotting.plot_weight_histogram() : The distribution of the weights.

Examples

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