# plasma_plots.fieldlines

*module*

Magnetic field lines on the mapped grid, traced, cut, classified and sampled along.

A field line follows the direction of a vector field, ``dx/ds = B/|B|`` with ``s`` the arc length.
Here it is traced in the logical coordinates of the grid (``eta1``, ``eta2``, ``eta3``, or GVEC's
``rho``, ``theta``, ``zeta``; see [`plasma_plots.arrays.logical_dims()`][plasma_plots.arrays.logical_dims]), where the grid is
rectangular and the angles wrap around: ``dη/ds = J⁻¹ B / |B|``, with the Jacobian ``J`` of the
mapping from the ``X``, ``Y``, ``Z`` coordinates. The contravariant unit field is interpolated
trilinearly between the grid points and integrated with a fourth-order Runge-Kutta scheme in fixed
steps of arc length. A line stops when it leaves the grid through a bounded direction (the radius,
or the ends of a torus sector), after a number of toroidal transits, or after a length.

The result is a Dataset over ``(s, line)``, like an orbits product is over ``(t, marker)``: the
logical coordinates, the positions ``x``, ``y``, ``z`` and ``|B|`` along each line, and per line
its rotational transform, its toroidal transits, where it ended and its connection length. The
punctures of a poloidal section are recorded while tracing, at every step, so a Poincaré plot does
not depend on how densely the line is saved. The other functions work on that Dataset:
[`poincare_section()`][poincare_section], [`classify_field_lines()`][classify_field_lines] and [`islands()`][islands], [`footprint()`][footprint] and
[`seed_grid()`][seed_grid], [`sample_along()`][sample_along] and [`parallel_wavenumber()`][parallel_wavenumber];
[`plasma_plots.fieldline_plots`][plasma_plots.fieldline_plots] draws them.

Trilinear interpolation of the direction field is second-order accurate in the grid spacing; a
Poincaré plot on a coarse grid shows that as a slow drift of the punctures, not as islands with a
rational transform. Refine the evaluation grid, not the step, when that matters.

**Examples**

```pycon
>>> lines = trace_field_lines(B.isel(t=-1), seeds=12, turns=50)
>>> lines.plasma.plot.poincare()  # R, Z punctures of the plane zeta = 0
>>> lines.iota  # the rotational transform of each line
>>> # the mode along one field line, and its parallel wavenumber
>>> along = sample_along(phi, lines)
>>> parallel_wavenumber(along)
```

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

## plasma_plots.fieldlines.LINE_CLASSES

*attribute* · *module attribute*

```python
LINE_CLASSES = {0: 'surface', 1: 'island', 2: 'chaotic'}
```

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

## plasma_plots.fieldlines.classify_field_lines

*function*

```python
def classify_field_lines(section: xr.Dataset, *, max_denominator: int = 12, tolerance: float | None = None, threshold: float = 0.1, min_spread: float | None = None) -> xr.DataArray
```

Classify each field line of a Poincaré section: on a flux surface (0), in an island (1) or chaotic (2).

A heuristic on the punctures of each line, in the order of their arc length. A line whose
punctures spread less than ``min_spread`` in the radial coordinate lies on a flux surface.
Otherwise its rotational transform picks the lowest-order rational ``n/m`` within
``tolerance`` (``n`` a multiple of nfp when the field has it), and the punctures are folded
into one island by the angle ``ψ = m θ mod period``. Successive punctures of a line on a
surface near the rational keep advancing in ``ψ`` (rotation: their accumulated advance grows
without bound), while those of a line inside an island turn back (libration about the
O-point): a line is in an island when its accumulated ``ψ`` advance never spans a full
period. Otherwise the punctures sorted by the poloidal angle should trace a smooth curve, the
flux surface; when the median radial jump between neighbours in the angle exceeds
``threshold`` times the line's radial spread, the line is chaotic. Lines with fewer than six punctures are surfaces. A line on a surface
whose ι is closer to the rational than one over the number of transits passes as an island;
trace longer to tell them apart.

**Parameters**

- `section` (`xarray.Dataset`) — A [`poincare_section()`][plasma_plots.fieldlines.poincare_section], or the lines of [`trace_field_lines()`][plasma_plots.fieldlines.trace_field_lines] (cut at their section).
- `max_denominator` (`int`) (default: `12`) — The largest ``m`` of the rationals considered. Default: 12.
- `tolerance` (`float`) (default: `None`) — How close ι must be to ``n/m``. Default: two over the number of toroidal transits of the line (the resolution of ι from the trace), at least 1e-3.
- `threshold` (`float`) (default: `0.1`) — The median radial jump between punctures that are neighbours in the poloidal angle, relative to the line's radial spread, above which a non-librating line is chaotic (a smooth curve through 100 punctures gives a few hundredths). Default: 0.1.
- `min_spread` (`float`) (default: `None`) — The radial spread of the punctures (in the radial logical coordinate) below which a line is on a surface. Default: half the radial grid spacing of the traced field.

**Returns**

- (`xarray.DataArray`) — The class of each line over ``line`` (``flag_meanings`` in the attrs, see [`LINE_CLASSES`][plasma_plots.fieldlines.LINE_CLASSES]), with the coordinates ``n`` and ``m`` of the rational found (0 when none) and ``iota``.

> **See Also**
>
> [`islands()`][plasma_plots.fieldlines.islands] : The island chains and their widths.

**Examples**

```pycon
>>> codes = classify_field_lines(poincare_section(lines))
>>> lines.isel(line=codes.values == 1)  # the lines inside islands
```

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

## plasma_plots.fieldlines.footprint

*function*

```python
def footprint(lines: xr.Dataset) -> xr.Dataset
```

Where the traced field lines left the grid, with their connection lengths.

**Parameters**

- `lines` (`xarray.Dataset`) — The lines of [`trace_field_lines()`][plasma_plots.fieldlines.trace_field_lines].

**Returns**

- (`xarray.Dataset`) — Over ``line``: the logical coordinates of the exit point (NaN for a line that did not leave the grid), ``x``, ``y``, ``z``, ``R``, ``connection_length`` and ``exited``, with the seeds as coordinates. The attrs name the dimensions.

> **See Also**
>
> [`plasma_plots.fieldline_plots.plot_footprint()`][plasma_plots.fieldline_plots.plot_footprint] : The exit points over the angles.
> [`plasma_plots.fieldline_plots.plot_connection_length()`][plasma_plots.fieldline_plots.plot_connection_length] : The connection lengths over the seeds.

**Examples**

```pycon
>>> hits = footprint(edge)
>>> hits.where(hits.exited, drop=True).to_dataframe()
```

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

## plasma_plots.fieldlines.islands

*function*

```python
def islands(section: xr.Dataset, *, max_denominator: int = 12, tolerance: float | None = None, threshold: float = 0.1, min_spread: float | None = None) -> xr.Dataset
```

The island chains a Poincaré section shows, with an estimate of their widths.

The lines [`classify_field_lines()`][plasma_plots.fieldlines.classify_field_lines] puts inside an island are grouped by their rational
``n/m``. The width of a chain is the radial extent of its widest island orbit, the one traced
closest to the separatrix (the full width of an island, at its O-point), so a chain traced
only near its O-point is underestimated: seed densely across it. The O-point's poloidal
angle is where that orbit's punctures are furthest apart radially, folded by ``m``; the
other O-points are a period over ``m`` apart.

**Parameters**

- `section` (`xarray.Dataset`) — A [`poincare_section()`][plasma_plots.fieldlines.poincare_section], or the lines of [`trace_field_lines()`][plasma_plots.fieldlines.trace_field_lines].
- `max_denominator` (`int`) (default: `12`) — The largest ``m`` of the rationals considered. Default: 12.
- `tolerance` (`float`) (default: `None`) — How close ι must be to ``n/m``; see [`classify_field_lines()`][plasma_plots.fieldlines.classify_field_lines].
- `threshold` (`float`) (default: `0.1`) — The chaos criterion of [`classify_field_lines()`][plasma_plots.fieldlines.classify_field_lines]. Default: 0.1.
- `min_spread` (`float`) (default: `None`) — The surface criterion of [`classify_field_lines()`][plasma_plots.fieldlines.classify_field_lines]. Default: half the radial grid spacing.

**Returns**

- (`xarray.Dataset`) — Over ``chain``: ``n``, ``m``, ``iota`` (``n/m``), ``width`` (in the radial logical coordinate), ``width_physical`` (the distance between the orbit's radially outermost and innermost punctures in ``x``, ``y``, ``z``), ``center`` (the mean radial coordinate of the chain's punctures), ``lines`` (how many lines lie in it) and ``o_point_<poloidal>`` (the poloidal angle of one O-point). Empty when there is no island.

> **See Also**
>
> [`classify_field_lines()`][plasma_plots.fieldlines.classify_field_lines] : Which lines are in islands.
> [`plasma_plots.fieldline_plots.plot_poincare()`][plasma_plots.fieldline_plots.plot_poincare] : ``islands=True`` annotates them.

**Examples**

```pycon
>>> chains = islands(poincare_section(lines))
>>> chains.to_dataframe()
```

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

## plasma_plots.fieldlines.parallel_wavenumber

*function*

```python
def parallel_wavenumber(samples: xr.DataArray, *, method: str = 'fft', detrend: bool = True) -> xr.DataArray
```

The dominant wavenumber ``k∥`` of a field along each traced field line.

From the samples of [`sample_along()`][plasma_plots.fieldlines.sample_along], over the arc length ``s`` of each line (its valid,
uniformly spaced part). ``"fft"`` takes the peak of the power spectrum (Hann window), refined
by a parabola through the peak and its neighbours; ``"crossings"`` counts the zero crossings,
``k∥ = π × crossings / length``, which suits a few wavelengths. Compare with
[`plasma_plots.theory.waves.parallel_wavenumber()`][plasma_plots.theory.waves.parallel_wavenumber], ``(n + m/q)/R₀`` for a mode
``(m, n)`` in a cylinder.

**Parameters**

- `samples` (`xarray.DataArray`) — The field along the lines, over ``(s, line)`` and any other dimensions.
- `method` (`('fft', 'crossings')`) (default: `"fft"`) — How to estimate it. Default: ``"fft"``.
- `detrend` (`bool`) (default: `True`) — Remove the mean along each line first. Default: True.

**Returns**

- (`xarray.DataArray`) — ``k_parallel`` over ``line`` and the other dimensions, NaN for a line with fewer than four samples.

**Raises**

- `ValueError` — If ``samples`` has no ``s`` dimension, or ``method`` is unknown.

**Examples**

```pycon
>>> k_par = parallel_wavenumber(sample_along(phi.isel(t=-1), lines))
```

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

## plasma_plots.fieldlines.poincare_section

*function*

```python
def poincare_section(lines: xr.Dataset, *, angle: float | None = None) -> xr.Dataset
```

The punctures of a poloidal plane by traced field lines, as a Dataset over ``(puncture, line)``.

With ``angle`` left out (or equal to the ``section`` the lines were traced with) the punctures
recorded while tracing are returned, which count every step. Another ``angle`` finds the
crossings in the saved samples, interpolated linearly between them, which is as fine as the
lines were saved (see ``stride`` of [`trace_field_lines()`][plasma_plots.fieldlines.trace_field_lines]). A section passed in comes back
as it is.

**Parameters**

- `lines` (`xarray.Dataset`) — The lines of [`trace_field_lines()`][plasma_plots.fieldlines.trace_field_lines].
- `angle` (`float`) (default: `None`) — The toroidal logical coordinate of the plane. Default: the traced section.

**Returns**

- (`xarray.Dataset`) — Over ``(puncture, line)``: the radial and poloidal logical coordinates of each puncture (named as in the field), ``x``, ``y``, ``z``, ``R`` (``√(x² + y²)``) and ``s`` (the arc length there), padded with NaN; per ``line`` the coordinates ``iota``, ``transits``, ``connection_length`` and the seeds. The attrs hold the ``section`` and the dimension names.

**Raises**

- `ValueError` — If the toroidal direction does not wrap around, so there is no plane to cut, or a section
is asked to be cut again at another angle.

> **See Also**
>
> [`classify_field_lines()`][plasma_plots.fieldlines.classify_field_lines], [`islands()`][plasma_plots.fieldlines.islands]
> [`plasma_plots.fieldline_plots.plot_poincare()`][plasma_plots.fieldline_plots.plot_poincare] : The plot.

**Examples**

```pycon
>>> section = poincare_section(lines)
>>> section.plasma.plot.poincare()
>>> poincare_section(lines, angle=np.pi / 5).plasma.plot.poincare()
```

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

## plasma_plots.fieldlines.rotational_transform

*function*

```python
def rotational_transform(lines: xr.Dataset) -> xr.DataArray
```

The rotational transform of each traced field line, poloidal turns per toroidal turn.

Computed while tracing from the whole line (its unwrapped poloidal and toroidal advance; a
toroidal turn is nfp periods of the toroidal coordinate: 2π for GVEC's angles, one period of
``eta3`` for Struphy's). NaN for a line that made no toroidal transit.

**Parameters**

- `lines` (`xarray.Dataset`) — The lines of [`trace_field_lines()`][plasma_plots.fieldlines.trace_field_lines], or a [`poincare_section()`][plasma_plots.fieldlines.poincare_section] of them.

**Returns**

- (`xarray.DataArray`) — ``iota`` over ``line``, with the seeds as coordinates.

**Raises**

- `ValueError` — If ``lines`` carries no ``iota``.

> **See Also**
>
> [`plasma_plots.analysis.rational_surfaces()`][plasma_plots.analysis.rational_surfaces] : Where a profile of ι is rational.

**Examples**

```pycon
>>> iota = rotational_transform(lines)
>>> iota.swap_dims(line="eta1_start").plasma.plot.lineout()
```

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

## plasma_plots.fieldlines.sample_along

*function*

```python
def sample_along(field: xr.DataArray, lines: xr.Dataset) -> xr.DataArray
```

Interpolate a scalar field along traced field lines, trilinearly on its logical grid.

**Parameters**

- `field` (`xarray.DataArray`) — A scalar field over the three logical dimensions the lines were traced in (its other dimensions, e.g. ``t``, are kept), with the ``X``, ``Y``, ``Z`` coordinates.
- `lines` (`xarray.Dataset`) — The lines of [`trace_field_lines()`][plasma_plots.fieldlines.trace_field_lines].

**Returns**

- (`xarray.DataArray`) — The field over ``(s, line)`` after the field's other dimensions, NaN where a line has stopped, named after the field and carrying its label and units; the ``s`` coordinate is the arc length, and the seeds, ``iota``, ``transits`` and ``connection_length`` are coordinates on ``line``.

**Raises**

- `ValueError` — If the field's logical dimensions differ from the lines', or it has a ``component``
dimension.

> **See Also**
>
> [`parallel_wavenumber()`][plasma_plots.fieldlines.parallel_wavenumber] : The dominant wavenumber along each line.
> [`plasma_plots.fieldline_plots.plot_along_field_lines()`][plasma_plots.fieldline_plots.plot_along_field_lines] : The profiles.

**Examples**

```pycon
>>> along = sample_along(phi, lines)  # (t, s, line)
>>> along.isel(line=0).plasma.plot.slice(x="s", y="t")
```

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

## plasma_plots.fieldlines.seed_grid

*function*

```python
def seed_grid(lines: xr.Dataset, name: str = 'connection_length') -> xr.DataArray
```

A per-line quantity over the two seed coordinates that vary, when the seeds form a grid.

**Parameters**

- `lines` (`xarray.Dataset`) — The lines of [`trace_field_lines()`][plasma_plots.fieldlines.trace_field_lines] (or their [`footprint()`][plasma_plots.fieldlines.footprint]), seeded with a dict of two 1-D arrays, so that the seeds form a grid.
- `name` (`str`) (default: `'connection_length'`) — The per-line variable. Default: ``"connection_length"``.

**Returns**

- (`xarray.DataArray`) — ``name`` over the two varying seed coordinates (named after the logical dimensions, e.g. ``(theta, zeta)``). With ``direction="both"`` the forward lines are taken (both carry the same connection length).

**Raises**

- `ValueError` — If the seeds do not vary along exactly two coordinates, or do not fill a grid.

**Examples**

```pycon
>>> seed_grid(edge, "connection_length").plasma.plot.slice()
```

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

## plasma_plots.fieldlines.trace_field_lines

*function*

```python
def trace_field_lines(data: xr.DataArray, *, seeds=8, turns: float | None = None, length: float | None = None, step: float | None = None, direction: str = 'forward', section: float | None = None, stride: int | None = None, components: str = 'cartesian', max_steps: int = 2000000) -> xr.Dataset
```

Trace field lines of a vector field through its mapped grid.

See [`plasma_plots.fieldlines`][plasma_plots.fieldlines] for the method. Each line starts at a seed and stops when it
has made ``turns`` toroidal transits, has reached the arc ``length``, or has left the grid
through a bounded direction (the radius, or the ends of a torus sector); afterwards its samples
are NaN. The punctures of the poloidal plane ``section`` (a value of the toroidal logical
coordinate) are recorded at every step while tracing, for [`poincare_section()`][plasma_plots.fieldlines.poincare_section].

**Parameters**

- `data` (`xarray.DataArray`) — The vector field, with a ``component`` dimension of size 3 and the three logical dimensions (every other dimension selected, e.g. ``B.isel(t=-1)``), with the ``X``, ``Y``, ``Z`` coordinates of the grid.
- `seeds` (`(int, dict, array_like or xarray.Dataset)`) (default: `8`) — Where the lines start, in logical coordinates: a number of seeds spread along the radial coordinate at the first poloidal and toroidal grid values (the outboard midplane of a torus whose angles start at 0); a dict of logical coordinates to values or 1-D arrays, which are combined into a grid of seeds (e.g. ``{"rho": 0.95, "theta": thetas, "zeta": zetas}`` for a connection-length map; a missing coordinate takes its first grid value); an ``(n, 3)`` array of points in the order of the logical dimensions; or a Dataset with one variable per logical coordinate. Default: 8.
- `turns` (`float`) (default: `None`) — Stop a line after this many toroidal transits (turns of the torus: nfp periods of the toroidal logical coordinate, i.e. 2π of GVEC's toroidal angle, one period of Struphy's ``eta3``). Default: 20 when the toroidal direction wraps around, else none.
- `length` (`float`) (default: `None`) — Stop a line after this arc length, in the units of ``X``, ``Y``, ``Z``. Default: with ``turns``, 1.5 times the length ``turns`` circles through the seeds' mean major radius would take (a safety net); without, four times the grid's extent.
- `step` (`float`) (default: `None`) — The step of arc length of the integrator. Default: the median spacing of the grid points (trilinear interpolation limits the accuracy before the step does).
- `direction` (`('forward', 'backward', 'both')`) (default: `"forward"`) — Along the field, against it, or each seed in both directions (two lines per seed, with the coordinates ``seed`` and ``direction`` telling them apart; the connection length is then the sum of both). Default: ``"forward"``.
- `section` (`float`) (default: `None`) — The toroidal logical coordinate of the poloidal plane whose punctures are recorded. Default: the first toroidal grid value (e.g. ``zeta = 0``).
- `stride` (`int`) (default: `None`) — Save every ``stride``-th step of each line. Default: as many as keep each line at about 20 000 samples at most; the punctures, transits and lengths count every step regardless.
- `components` (`('cartesian', 'contravariant')`) (default: `"cartesian"`) — What the components are: ``"cartesian"`` (default) ``(x, y, z)``, or contravariant logical components, as for [`plasma_plots.analysis.divergence()`][plasma_plots.analysis.divergence].
- `max_steps` (`int`) (default: `2000000`) — A cap on the number of steps of each line. Default: 2 000 000.

**Returns**

- (`xarray.Dataset`) — Over ``(s, line)``: the logical coordinates of each line (angles folded into one period), its positions ``x``, ``y``, ``z`` and ``absB`` (``|B|``), NaN after a line has stopped; over ``line``: ``iota`` (poloidal per toroidal turns, from the whole traced line; NaN without a toroidal transit), ``transits``, ``length`` (the arc length traced), ``exited`` (whether it left the grid), ``connection_length`` (the arc length to the exit, summed over both directions with ``direction="both"``; NaN while a line has not left the grid), and where each line ended, ``<dim>_end`` and ``x_end``, ``y_end``, ``z_end``; over ``(puncture, line)``: ``puncture_<radial>``, ``puncture_<poloidal>``, ``puncture_x``, ``puncture_y``, ``puncture_z`` and ``puncture_s``, where each line crosses ``section`` (padded with NaN). The coordinates ``<dim>_start``, ``seed`` and ``direction`` describe the seeds; the attrs hold the section, the dimension names, the periods and the toroidal turn.

**Raises**

- `ValueError` — If ``data`` is not a three-component field over the three logical dimensions with ``X``,
``Y``, ``Z``, has other dimensions left, or the seeds or ``direction`` are malformed.

> **See Also**
>
> [`poincare_section()`][plasma_plots.fieldlines.poincare_section] : The punctures as a Dataset of their own.
> [`footprint()`][plasma_plots.fieldlines.footprint] : Where the lines that left the grid did so.
> [`sample_along()`][plasma_plots.fieldlines.sample_along] : A field along the lines.
> [`plasma_plots.fieldline_plots.plot_poincare()`][plasma_plots.fieldline_plots.plot_poincare] : The Poincaré plot.

**Examples**

```pycon
>>> lines = trace_field_lines(B.isel(t=-1), seeds=12, turns=100)
>>> lines.iota.values  # the rotational transform of each line
>>> # connection lengths from a surface near the edge, both ways
>>> edge = trace_field_lines(
...     B.isel(t=-1),
...     seeds={"eta1": 0.98, "eta2": np.linspace(0, 1, 32), "eta3": 0.0},
...     direction="both",
...     turns=50,
... )
```

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