"""Provides functionality of chord/line segment lineouts."""
import numpy as np
import warnings
from numpy.polynomial import Polynomial
from torx.grid.grid_2d_m import Grid2D
from torx.equilibrium.equilibrium_m import EquilibriumBaseClass
from torx.analysis.lineouts.lineout_m import Lineout
from torx.analysis.lineouts.flux_surface_lineouts_m import FluxSurfaceLineout
from torx.geometry.line_geometry_m import intersect_lines
from torx.autodoc_decorators_m import autodoc_function
[docs]
@autodoc_function
def outboard_zaxis_chord(
grid: Grid2D,
equi: EquilibriumBaseClass,
n_samples: int = 1000
):
"""
Return equally spaced points along the outboard zaxis.
It is a horizontal line from the magnetic axis outwards to large radius at
geometric poloidal angle theta equals zero.
"""
r_start = equi.axis_r_norm
z_start = equi.axis_z_norm
r_end = grid.r_s.max()
z_end = equi.axis_z_norm
lineout = Lineout([r_start, r_end], [z_start, z_end])
lineout.find_points_on_grid(grid, n_samples=n_samples)
return lineout
[docs]
@autodoc_function
def outboard_midplane_chord(
grid: Grid2D,
equi: EquilibriumBaseClass,
n_samples: int = 1000,
separatrix_index: int = None,
):
"""
Return equally spaced points at the outboard midplane.
It is defined as a line going from the magnetic axis to the largest radial
position of the separatrix.
Shift very slightly inside to get just the confined-region separatrix.
"""
# TODO: Refactor the following block to be simpler
matched_contours = equi.get_separatrix(grid.r_s, grid.z_s, level=0.999)
if len(matched_contours) == 0:
raise ValueError("No separatrix found")
elif len(matched_contours) > 1:
if separatrix_index is None:
# Try auto find the correct separatrix
closed = np.zeros(len(matched_contours), dtype=int)
for i, _ in enumerate(matched_contours):
# Checks if start and end point of the contour are the same
if(np.abs(matched_contours[i][0,0] \
- matched_contours[i][-1,0]) < 1e-14 and
np.abs(matched_contours[i][0,1] \
- matched_contours[i][-1,1]) < 1e-14):
closed[i] = 1
# If single closed contour found close to the separatrix we take
# that. Otherwise the choice is ambiguous and separatrix_index has
# to be used.
# NOTE: Even if multiple open contours and no closed one are found
# it is not clear if get_separatrix failed or not.
if np.sum(closed) == 1:
idx_sep = np.where(closed == 1)[0][0]
matched_contours = matched_contours[idx_sep]
else:
raise RuntimeError(
f"Found {len(matched_contours)} contours for \
get_separatrix and could not auto-identify the correct \
one. Please identify which contour to use via the \
separatrix_index argument."
)
else:
matched_contours = matched_contours[separatrix_index]
elif separatrix_index not in [0, None]:
warnings.warn(
"Only a single contour found. Ignoring separatrix_index argument",
UserWarning,
)
matched_contours = matched_contours[0]
else:
matched_contours = matched_contours[0]
r_sep, z_sep = matched_contours.T
i_omp = np.argmax(r_sep)
r_start = equi.axis_r_norm
z_start = equi.axis_z_norm
# Linearly extrapolate to the end of the grid
# NOTE: this will fail if you try to do it for a vertical line
linear_extrap = Polynomial.fit(
(r_start, r_sep[i_omp]), (z_start, z_sep[i_omp]), deg=1
)
lineout = Lineout(
[r_start, grid.r_s.max()], [z_start, linear_extrap(grid.r_s.max())]
)
lineout.find_points_on_grid(grid, n_samples=n_samples)
return lineout
[docs]
@autodoc_function
def radial_line_chord(grid: Grid2D,
equi: EquilibriumBaseClass,
bounding_flux_surface: FluxSurfaceLineout,
n_samples: int=100,
theta=0.0):
"""
Return a radial lineout from the axis to a bounding flux surface.
Creates a straight radial line along a poloidal angle starting from the
magnetic axis. The end of the line is specified by a bounding flux surface
lineout which can be either a closed or open flux surface.
"""
# Idea: Create a line that goes far from the starting point along the
# specified poloidal angle. Then intersect with the bounding flux surface
# to obtain a line from start to the flux surface.
r_start = equi.axis_r_norm.values
z_start = equi.axis_z_norm.values
bounding_r = bounding_flux_surface.r_points
bounding_z = bounding_flux_surface.z_points
# Use twice the maximum distance of the bounding flux surface from the
# magnetic axis for the far away point
dist_from_axis = np.sqrt((bounding_r - r_start)**2 \
+ (bounding_z - z_start)**2)
radius = 2.0 * np.max(dist_from_axis)
r_end = r_start + radius * np.cos(theta)
z_end = z_start + radius * np.sin(theta)
# Close the bounding line by adding the first element to the end
# (required for correct intersection). Do only for closed flux surface.
if(bounding_flux_surface.rho(equi)[0] <= 1.0):
bounding_r = np.append(bounding_r, bounding_r[0])
bounding_z = np.append(bounding_z, bounding_z[0])
pint = intersect_lines(np.array([r_start, r_end]),
np.array([z_start, z_end]),
bounding_r,
bounding_z,
intersec_mode="closest")
r_end = pint[0]
z_end = pint[1]
lineout = Lineout([r_start, r_end], [z_start, z_end])
lineout.find_points_on_grid(grid, n_samples=n_samples)
return lineout
