"""Classes for tracing magnetic field lines."""
import warnings
import xarray as xr
import numpy as np
from typing import Union, List, Callable, Sequence
from scipy.integrate import solve_ivp
from scipy.integrate._ivp.ivp import OdeResult
from scipy.interpolate import InterpolatedUnivariateSpline
from .vector_field_tracer_m import VectorFieldTracer
from torx.decorators import autodoc_class, autodoc_function
from torx.arrays import make_xarray
[docs]
@autodoc_class
class MagneticFieldTracer(VectorFieldTracer):
"""
Allows to trace along a magnetic fieldline.
Toroidal angle is denoted as phi. Poloidal angle as theta.
"""
[docs]
def __init__(self, equi, max_flux_surface_wander: float=1e-3):
"""Initialize the magnetic field tracer."""
# Only for NumericalEquilibrium the init of the parent class can be
# called with a the discrete magnetic vector field
if type(equi).__name__ == "NumericalEquilibrium":
r_norm, z_norm = equi.spline_basis_r, equi.spline_basis_z
vec_field = equi.magfield_vector(r_norm, z_norm)
super().__init__(vec_field)
self.equi = equi
self.max_flux_surface_wander = max_flux_surface_wander
[docs]
def toroidal_integration(
self,
r_initial: float,
z_initial: float,
phi_initial: float,
phi_to_max: float,
tolerance: float=1e-4,
method: str="DOP853",
direction: str="fwd",
check_out_of_domain: bool=False,
check_out_of_bounds: bool=False,
events: Union[List[Callable], Callable]=None,
**integrator_kwargs
) -> OdeResult:
"""
Perform an integration with phi as the independent variable.
Returns the (x,y) positions of the trace, as well as the fieldline
length.
"""
self.direction = direction
self._check_initial_state(r_initial, z_initial, phi_initial)
# Handle events and check format of the provided events
if check_out_of_domain:
events = self._add_event(events, self._out_of_domain)
if check_out_of_bounds:
events = self._add_event(events, self._out_of_bounds)
solution = solve_ivp(
fun=self.toroidal_integration_equation,
t_span=[phi_initial, phi_initial + self._dir_fac * phi_to_max],
y0=np.array([r_initial, z_initial, 0.0]),
events=events,
dense_output=True,
rtol=tolerance,
atol=tolerance,
method=method,
**integrator_kwargs,
)
if not np.isnan(self.max_flux_surface_wander):
self.check_flux_surface_wander(solution)
return solution
[docs]
def toroidal_integration_equation(self, phi, state):
"""
Equation that is required to integrate toroidally around the torus.
The fieldline length is returned as the third element of the state
vector.
"""
r_norm = state[0]
z_norm = state[1]
b_r, b_z, b_t, jacobian = self.equi._magfield_scalars(r_norm,
z_norm, phi)
d_state = np.zeros_like(state)
d_state[0] = b_r / b_t * jacobian
d_state[1] = b_z / b_t * jacobian
d_state[2] = (
np.sqrt(b_r * b_r / (b_t * b_t) + b_z * b_z / (b_t * b_t) + 1.0) \
* jacobian
)
return d_state
[docs]
def poloidal_integration(
self,
r_initial: float,
z_initial: float,
theta_max: float,
tolerance: float=1e-4,
method: str="DOP853",
check_out_of_domain: bool=False,
check_out_of_bounds: bool=False,
events: Union[List[Callable], Callable]=None,
**integrator_kwargs,
) -> OdeResult:
"""
Perform an integration with theta as the independent variable.
Valid only for the confined region.
Returns the (x,y) positions of the trace, as well as the fieldline
length.
"""
self._check_initial_state(r_initial, z_initial, 0.0)
# Handle events and check format of the provided events
if check_out_of_domain:
events = self._add_event(events, self._out_of_domain)
if check_out_of_bounds:
events = self._add_event(events, self._out_of_bounds)
# Make sure we are in the closed field line region
rho_initial = self.equi.normalized_flux_surface_label(r_initial,
z_initial).values
if rho_initial >= 1.0:
raise RuntimeError(
"Cannot perform poloidal integration for open fieldlines." \
+f"x={r_initial}, y= {z_initial}, rho={rho_initial}"
)
solution = solve_ivp(
fun=self.poloidal_integration_equation,
t_span=[0.0, theta_max],
y0=np.array([r_initial, z_initial, 0.0]),
events=events,
dense_output=True,
rtol=tolerance,
atol=tolerance,
method=method,
**integrator_kwargs,
)
self.check_flux_surface_wander(solution)
return solution
[docs]
def poloidal_integration_equation(self, _, state):
"""
Equation that is required to integrate poloidally around the torus.
The fieldline length is returned as the third element of the state
vector.
"""
r_norm = state[0]
z_norm = state[1]
b_r = self.equi.magfield_component_r(r_norm, z_norm)
b_z = self.equi.magfield_component_z(r_norm, z_norm)
theta = self.theta(r_norm, z_norm)
radius = self.minor_radius(r_norm, z_norm)
# Magnetic field component along the integration contour (theta)
b_along_fs = -b_r * np.sin(theta) + b_z * np.cos(theta)
# First check magnetic field obtained
b_pol = np.sqrt(b_r**2 + b_z**2)
# Magnetic field component across the integration contour (radial)
b_across_fs = b_r * np.cos(theta) + b_z * np.sin(theta)
b_pol_check = np.sqrt(b_across_fs**2 + b_along_fs**2)
assert np.isclose(b_pol_check / b_pol, 1.0
), f"Poloidal field on curve was {b_pol_check}, should be {b_pol}!"
# After check passed, calculate increments
d_state = np.zeros_like(state)
# The increments dR and dZ are given by simply integrating the
# equations for the magnetic field line (4.1.6) in [1] on
# the poloidal curve.
# [1] D'haeseleer, Flux coordinates and magnetic field structure
# NOTE: The poloidal line element is r*dtheta.
d_state[0] = radius * b_r / np.abs(b_along_fs)
d_state[1] = radius * b_z / np.abs(b_along_fs)
# The increment in symmetry angle is given by eq. (24) in [2].
# [2] Ribeiro, Conformal Tokamak Geometry for Turbulence Computations
# NOTE: Here we return the product q*dtheta_s. Since the contravariant
# component of the poloidal field is used, a factor 1/r enters the
# denominator.
# NOTE: The current I used in the formula is I = B_tor * R = B_tor_axis
# since we set B_tor = B_tor_axis / R per definition. We use the
# call to the toroidal field with axis coordinates since this
# returns the field with the same normalization as used above.
current = self.equi.magfield_component_toroidal(self.equi.axis_r_norm,
self.equi.axis_z_norm)
d_state[2] = current * np.abs(radius / b_along_fs / r_norm**2)
return d_state
[docs]
def check_flux_surface_wander(self, solution):
"""
Check the maximum deviation from the flux-surface during trace.
Throws an error if it exceeds the threshold.
"""
rho_initial = self.equi.normalized_flux_surface_label(
solution.y[0][0], solution.y[1][0]
).values
dense_output = solution.sol(
np.linspace(np.min(solution.t), np.max(solution.t), 1000)
)
rho_trace = self.equi.normalized_flux_surface_label(
dense_output[0], dense_output[1],
).values
flux_surface_wander = np.abs(rho_trace - rho_initial)
assert (
np.max(flux_surface_wander) < self.max_flux_surface_wander
), (
f"Flux surface wander exceeded threshold "
f"({np.max(flux_surface_wander):3.2e}"
f" > {self.max_flux_surface_wander:3.2e}) "
"Try decreasing the tolerance of the integrator, "
"or increasing the max_flux_surface_wander"
)
[docs]
def find_neighboring_points(
self, r_initial, z_initial, n_toroidal_planes: int=16
):
"""
Find the neighbors in both directions of an array of sample points.
Although this loop should be trivially parallel, dask parallelism
doesn't seem to work here.
"""
r_initial = np.atleast_1d(r_initial)
z_initial = np.atleast_1d(z_initial)
assert r_initial.shape == z_initial.shape
assert (
r_initial.ndim == 1 and z_initial.ndim == 1
), f"Should provide r and z as 1D arrays"
forward_trace = np.zeros((r_initial.size, 3))
reverse_trace = np.zeros((r_initial.size, 3))
print("Tracing", end=" ")
for i in range(r_initial.size):
print(f"{i}/{r_initial.size}", end=", ")
r_in, z_in = r_initial[i], z_initial[i]
forward_trace[i, :] = self.toroidal_integration(
r_in,
z_in,
+2.0 * np.pi / n_toroidal_planes,
).y[:, -1]
reverse_trace[i, :] = self.toroidal_integration(
r_in,
z_in,
-2.0 * np.pi / n_toroidal_planes,
).y[:, -1]
print("Done")
return forward_trace, reverse_trace
[docs]
def trace_reverse(self, r_initial, z_initial, n_phi_total, phi_initial = 0,
n_phi=1):
"""
Trace field lines back for a given amount of toroidal planes.
Starts at r,z initial and returns the new points r,z reverse.
"""
r_reverse = np.zeros(r_initial.size)
z_reverse = np.zeros(r_initial.size)
for i in range(r_initial.size):
r_in, z_in = r_initial[i], z_initial[i]
sol = self.toroidal_integration(
r_in,
z_in,
phi_initial,
(-2.0 * np.pi / n_phi_total) * n_phi,
)
r_reverse[i] = sol.y[0, -1]
z_reverse[i] = sol.y[1, -1]
return r_reverse, z_reverse
[docs]
def minor_radius(self, r_norm, z_norm):
"""Minor radius defined as Cartesian distance to magnetic axis."""
return np.sqrt((r_norm - self.equi.axis_r_norm.values)**2 +
(z_norm - self.equi.axis_z_norm.values)**2)
[docs]
def theta(self, r_norm, z_norm):
"""
Poloidal angle from the magnetic axis.
Note
----
The zero point is the outboard midplane, and theta ranges
from [-pi, pi].
"""
return np.arctan2(z_norm - self.equi.axis_z_norm.values,
r_norm - self.equi.axis_r_norm.values)
def _out_of_bounds(self, _, state):
"""
Check if integrator state is out of bounds.
Return 1 if the state of the integrator is inside the boundary
polygon, otherwise -1
"""
r, z = state[0], state[1]
boundary_polygon = self.equi.get_boundary_polygon()
return 1 if boundary_polygon.point_inside(r, z) else -1
_out_of_bounds.terminal = True
_out_of_bounds.direction = -1
@autodoc_class
class NormalTracer(MagneticFieldTracer):
"""Finds a line which is continuously normal to the poloidal field."""
def toroidal_integration_equation(self, _, state):
"""
Equation that is required to integrate toroidally around the torus.
The fieldline length is returned as the third element of the state
vector.
"""
r_norm, z_norm = state[0], state[1]
b_r, b_z, b_t, jacobian = (
self.equi.magfield_component_r(r_norm, z_norm),
self.equi.magfield_component_z(r_norm, z_norm),
self.equi.magfield_component_toroidal(r_norm, z_norm),
self.equi.jacobian(r_norm, z_norm),
)
d_state = np.zeros_like(state)
d_state[0] = -b_z / b_t * jacobian
d_state[1] = +b_r / b_t * jacobian
d_state[2] = 0
return d_state
def check_flux_surface_wander(self, solution):
"""
Check the maximum deviation from the flux-surface during trace.
Since we're perpendicularly tracing, we necessarily will go off the
flux surface.
"""
pass
[docs]
@autodoc_function
def polar_to_cart_chord(equi, rho, theta, max_radius=1.0, n_samples=1000):
"""Convert from polar to cartesian coordinates."""
xaxis = np.array(equi.axis_r_norm)
yaxis = np.array(equi.axis_z_norm)
radius = np.linspace(0, max_radius, n_samples)
r_norm = radius * np.cos(theta) + xaxis
z_norm = radius * np.sin(theta) + yaxis
rho_values = equi.normalized_flux_surface_label(r_norm, z_norm)
# NOTE: We manually set the first entry to zero due to possible very small
# inaccuracies that might lead to finite rho at axis
rho_values[0] = 0.0
xsol, ysol, rsol = [], [], []
for rhoval in np.atleast_1d(rho):
roots = InterpolatedUnivariateSpline(radius,
rho_values - rhoval).roots()
if len(roots) == 0:
print(f"rho value {rhoval} not found. Skipping")
continue
elif len(roots) > 1:
raise RuntimeError(
"Multiple roots found. Set max radius to ensure that rho is " \
"monotonic"
)
xsol.append(roots[0] * np.cos(theta) + xaxis)
ysol.append(roots[0] * np.sin(theta) + yaxis)
rsol.append(rhoval)
return np.array(xsol), np.array(ysol), np.array(rsol)
[docs]
def flare_fieldline_trace(
equi,
r_initial: float,
z_initial: float,
phi_initial: float=0.0,
direction: int=1,
step_size: float=0,
nsteps: int=16,
stop_at_boundary: bool=True,
coordinates: str="cylindrical",
angular_units: str="rad",
normalize: bool=True,
):
"""
Trace a magnetic field line using the FLARE field-line tracer.
This is a wrapper around ``flare.analysis.fieldline_trace`` with a more
convenient interface for the FLARE equilibrium. The initial point is
specified by ``(r_initial, z_initial, phi_initial)``.
If ``normalize=True`` (default), ``r_initial`` and ``z_initial`` are
assumed to be given as normalized coordinates and are converted to
the coordinates expected by FLARE (real space) before tracing.
The tracing direction is specified by ``direction``:
- ``+1`` for forward tracing
- ``-1`` for backward tracing
Parameters
----------
equi : FlareEquilibrium
Equilibrium providing the FLARE backend and the normalization used
to convert between normalized and real-space coordinates.
r_initial : float
Initial radial coordinate of the field line.
z_initial : float
Initial vertical coordinate of the field line.
phi_initial : float, default=0.0
Initial toroidal angle in radians.
direction : int, default=1
Tracing direction: ``+1`` forward, ``-1`` backward.
step_size : float, default=0
Integration step size forwarded to the FLARE backend.
nsteps : int, default=16
Number of integration steps.
stop_at_boundary : bool, default=True
Whether to stop the trace when it reaches the boundary.
coordinates : str, default="cylindrical"
Coordinate system forwarded to the FLARE backend.
angular_units : str, default="rad"
Angular units forwarded to the FLARE backend.
normalize : bool, default=True
If True, treat the initial point as normalized coordinates and
return the traced coordinates in normalized units.
Returns
-------
tuple[numpy.ndarray, numpy.ndarray, int]
The traced coordinates ``(R, Z, phi)``, the field-line length, and
the FLARE exit code.
Note
----
This function forwards most tracing options directly to the FLARE
backend. See the FLARE documentation for the detailed meaning of
``step_size``, ``coordinates``, ``angular_units``, and the returned
``exit_code``.
"""
from torx.specializations.flare import import_flare
analysis = import_flare()[2]
if normalize:
r_initial, z_initial = equi._denormalize_flare(r_initial, z_initial)
x0 = np.array([r_initial, z_initial, phi_initial], dtype=float)
coords, fieldline_length, exit_code = analysis.fieldline_trace(
x0,
idir=direction,
ds=step_size,
nsteps=nsteps,
stop_at_boundary=stop_at_boundary,
coordinates=coordinates,
angular_units=angular_units,
)
if normalize:
coords[0], coords[1] = equi._normalize_flare(coords[0], coords[1])
return coords, fieldline_length, exit_code
def _bidirectional_toroidal_tracer_flare(
equi,
r_initial: float,
z_initial: float,
phi_initial: float,
n_turns: int,
n_planes: int
):
"""
Trace a single field line using the FLARE backend.
Returns
-------
tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, int, int]
Arrays `(R, Z, phi, s)` of shape `(2 * n_turns * n_planes + 1,)`,
followed by backward/forward status codes.
"""
n_steps_per_direction = n_turns * n_planes
dphi = 2.0 * np.pi / n_planes
xyz_bwd, s_bwd, istat_bwd = flare_fieldline_trace(
equi,
r_initial=r_initial,
z_initial=z_initial,
phi_initial=phi_initial,
direction=-1,
step_size=dphi,
nsteps=n_steps_per_direction,
coordinates="cylindrical",
angular_units="rad",
)
xyz_fwd, s_fwd, istat_fwd = flare_fieldline_trace(
equi,
r_initial=r_initial,
z_initial=z_initial,
phi_initial=phi_initial,
direction=1,
step_size=dphi,
nsteps=n_steps_per_direction,
coordinates="cylindrical",
angular_units="rad",
)
R_bwd, Z_bwd, phi_bwd = xyz_bwd
R_fwd, Z_fwd, phi_fwd = xyz_fwd
# Reverse backward branch so the full trace runs from
# phi_min -> phi_initial -> phi_max.
R_bwd = R_bwd[::-1]
Z_bwd = Z_bwd[::-1]
phi_bwd = phi_bwd[::-1]
s_bwd = s_bwd[::-1]
# Re-center arc length so the initial point has s = 0.
s_bwd = s_bwd - s_bwd[-1]
s_fwd = s_fwd - s_fwd[0]
# Combine both branches, skipping duplicate initial point
R = np.concatenate([R_bwd, R_fwd[1:]])
Z = np.concatenate([Z_bwd, Z_fwd[1:]])
phi = np.concatenate([phi_bwd, phi_fwd[1:]])
s = np.concatenate([s_bwd, s_fwd[1:]])
return R, Z, phi, s, istat_bwd, istat_fwd
def _bidirectional_toroidal_tracer(
tracer,
r_initial: float,
z_initial: float,
phi_initial: float,
n_turns: int,
n_planes: int,
):
"""
Trace a single field line using the standard MagneticFieldTracer backend.
Returns
-------
tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, bool, str | None]
Arrays `(R, Z, phi, s)`, a success flag, and an error message.
"""
n_steps = n_turns * n_planes
dphi = 2.0 * np.pi / n_planes
t_eval_fwd = phi_initial + dphi * np.arange(1, n_steps + 1)
t_eval_bwd = phi_initial - dphi * np.arange(1, n_steps + 1)
t_eval_fwd[-1] = phi_initial + n_turns * 2.0 * np.pi
t_eval_bwd[-1] = phi_initial - n_turns * 2.0 * np.pi
y0 = np.array([r_initial, z_initial, 0.0])
needs_dense = not np.isnan(tracer.max_flux_surface_wander)
# Stop early if R leaves the physically valid range.
def r_nonpositive(_, state):
return state[0]
r_nonpositive.terminal = True
r_nonpositive.direction = -1
r_escape_threshold = 100.0 * max(abs(r_initial), 1.0)
def r_escaped(_, state):
return r_escape_threshold - state[0]
r_escaped.terminal = True
r_escaped.direction = -1
ivp_kw = dict(
rtol=1e-4, atol=1e-4, method="DOP853",
dense_output=needs_dense,
events=[r_nonpositive, r_escaped],
)
phi = np.concatenate([t_eval_bwd[::-1], [phi_initial], t_eval_fwd])
n_total = 2 * n_steps + 1
result_fwd = solve_ivp(
fun=tracer.toroidal_integration_equation,
t_span=[phi_initial, phi_initial + n_turns * 2.0 * np.pi],
y0=y0,
t_eval=t_eval_fwd,
**ivp_kw,
)
if result_fwd.status != 0:
message = (
"Terminated: R left valid range" if result_fwd.status == 1
else result_fwd.message
)
return (
np.full(n_total, np.nan),
np.full(n_total, np.nan),
phi,
np.full(n_total, np.nan),
False,
message,
)
result_bwd = solve_ivp(
fun=tracer.toroidal_integration_equation,
t_span=[phi_initial, phi_initial - n_turns * 2.0 * np.pi],
y0=y0,
t_eval=t_eval_bwd,
**ivp_kw,
)
if result_bwd.status != 0:
message = (
"Terminated: R left valid range" if result_bwd.status == 1
else result_bwd.message
)
return (
np.full(n_total, np.nan),
np.full(n_total, np.nan),
phi,
np.full(n_total, np.nan),
False,
message,
)
if needs_dense:
tracer.check_flux_surface_wander(result_fwd)
tracer.check_flux_surface_wander(result_bwd)
# Backward branch is reversed so index 0 is furthest back. Arc length is
# negative (phi decreasing), so negate to get positive distances from start.
R = np.concatenate([result_bwd.y[0, ::-1], [r_initial], result_fwd.y[0, :]])
Z = np.concatenate([result_bwd.y[1, ::-1], [z_initial], result_fwd.y[1, :]])
s = np.concatenate([-result_bwd.y[2, ::-1], [0.0], result_fwd.y[2, :]])
return R, Z, phi, s, True, None
[docs]
def bidirectional_toroidal_trace(
equi,
r_initial: float | Sequence[float],
z_initial: float | Sequence[float],
phi_initial: float,
n_turns: int = 2,
n_planes: int = 50,
max_flux_surface_wander: float = np.nan,
verbosity: int = 0,
) -> tuple[xr.DataArray, xr.DataArray, xr.DataArray, xr.DataArray]:
"""
Trace magnetic field lines in both toroidal directions.
For each initial point `(r_initial, z_initial)` at toroidal angle
`phi_initial`, the field line is traced backward and forward for
`n_turns` full toroidal turns. Sampling is performed at `n_planes`
equally spaced toroidal planes per turn.
If `equi` is a FlareEquilibrium instance, the FLARE backend is used.
Otherwise the standard `MagneticFieldTracer` backend is used.
Parameters
----------
equi
Equilibrium object.
r_initial : float or Sequence[float]
Initial radial coordinate(s).
z_initial : float or Sequence[float]
Initial vertical coordinate(s). Must have the same length as
`r_initial`.
phi_initial : float
Initial toroidal angle in radians.
n_turns : int, default=2
Number of full toroidal turns to trace in each direction.
n_planes : int, default=50
Number of sampling planes per turn.
max_flux_surface_wander : float, default is NaN
Passed to `MagneticFieldTracer` for non-FLARE equilibria.
verbosity : int, default=0
Verbosity level.
Returns
-------
tuple[
xarray.DataArray, xarray.DataArray,
xarray.DataArray, xarray.DataArray,
]
`R`, `Z`, `phi`, and `s`, each with dimensions `("line", "plane")`.
The `s` array is the field-line length centered such that the initial
point has value zero, with negative values on the backward branch and
positive values on the forward branch.
The returned `phi` array is in radians.
"""
from torx.equilibrium import FlareEquilibrium
n_turns = int(n_turns)
n_planes = int(n_planes)
if n_turns < 1:
raise ValueError("n_turns must be >= 1")
if n_planes < 1:
raise ValueError("n_planes must be >= 1")
r_initial = np.atleast_1d(np.asarray(r_initial, dtype=float))
z_initial = np.atleast_1d(np.asarray(z_initial, dtype=float))
if len(r_initial) != len(z_initial):
raise ValueError("r_initial and z_initial must have the same length")
n_field_lines = len(r_initial)
n_steps_total = 2 * n_turns * n_planes
shape = (n_field_lines, n_steps_total + 1)
use_flare = isinstance(equi, FlareEquilibrium)
backend = "FLARE" if use_flare else "standard"
R_data = np.full(shape, np.nan)
Z_data = np.full(shape, np.nan)
phi_data = np.full(shape, np.nan)
s_data = np.full(shape, np.nan)
tracer = None
if not use_flare:
tracer = MagneticFieldTracer(
equi,
max_flux_surface_wander=max_flux_surface_wander,
)
for l in range(n_field_lines):
if verbosity >= 1:
print(f"field line {l + 1} / {n_field_lines}")
if use_flare:
R, Z, phi, s, istat_bwd, istat_fwd = (
_bidirectional_toroidal_tracer_flare(
equi=equi,
r_initial=r_initial[l],
z_initial=z_initial[l],
phi_initial=phi_initial,
n_turns=n_turns,
n_planes=n_planes,
)
)
if istat_bwd != 0:
print(
f"Backward trace error for field line {l + 1}: "
f"istat = {istat_bwd}"
)
if istat_fwd != 0:
print(
f"Forward trace error for field line {l + 1}: "
f"istat = {istat_fwd}"
)
else:
R, Z, phi, s, success, message = _bidirectional_toroidal_tracer(
tracer=tracer,
r_initial=r_initial[l],
z_initial=z_initial[l],
phi_initial=phi_initial,
n_turns=n_turns,
n_planes=n_planes,
)
if not success:
if "R left valid range" in message:
warnings.warn(
f"Field line {l + 1} / {n_field_lines}: "
f"R left valid range before completing "
f"{n_turns} turns."
)
else:
warnings.warn(
f"Field line {l + 1} / {n_field_lines} trace "
f"failed: {message}"
)
R_data[l, :] = R
Z_data[l, :] = Z
phi_data[l, :] = phi
s_data[l, :] = s
attrs = {
"n_turns": n_turns,
"n_planes": n_planes,
"backend": backend,
}
dims = ("line", "plane")
coords = {
"line": np.arange(shape[0]),
"plane": np.arange(shape[1]),
}
arrays = (R_data, Z_data, phi_data, s_data)
names = ("R", "Z", "phi", "s")
R_da, Z_da, phi_da, s_da = tuple(
make_xarray(
array,
name=name,
dims=dims,
coords=coords,
attrs=attrs.copy(),
)
for array, name in zip(arrays, names, strict=True)
)
return R_da, Z_da, phi_da, s_da
