import numpy as np
from shapely import geometry
from pleque.core import Coordinates
from pleque.utils.decorators import *
[docs]class Surface(Coordinates):
[docs] def __init__(self, equilibrium, *coordinates, coord_type=None, grid=False, **coords):
"""
Calculates geometrical properties of a specified surface. To make the contour closed, the first and last points in
the passed coordinates have to be the same.
Instance is obtained by calling method `surface` in instance of `Equilibrium`.
:param coords: Instance of coordinate class
"""
super().__init__(equilibrium, *coordinates, coord_type=None, grid=False, **coords)
points_RZ = self.as_array(('R', 'Z'))
# closed surface has to have identical first and last points and then the shape is polygon
# opened surface is linestring
if np.isclose(points_RZ[0, 0], points_RZ[-1, 0]) and np.isclose(points_RZ[0, 1], points_RZ[-1, 1]):
self._poly = geometry.polygon.Polygon(points_RZ)
self._string = geometry.linestring.LineString(points_RZ)
self._closed = True
else:
self._string = geometry.linestring.LineString(points_RZ)
self._closed = False
@property
def closed(self):
"""
True if the fluxsurface is closed.
:return:
"""
return self._closed
@property
def area(self):
"""
Area of the closed fluxsurface.
:return:
"""
if self._closed:
return self._poly.area
else:
raise Exception("Opened Flux Surface does not have area")
@property
def length(self):
"""
Length of the fluxsurface contour
:return:
"""
return self._string.length
@property
def surface(self):
"""
Surface of fluxsurface calculated from the contour length
using Pappus centroid theorem : https://en.wikipedia.org/wiki/Pappus%27s_centroid_theorem
:return: float
"""
return self._string.length * 2 * np.pi * self.centroid.R[0]
@property
def centroid(self):
return self._eq.coordinates(R=np.array(self._string.centroid.coords)[0][0],
Z=np.array(self._string.centroid.coords)[0][0], coord_type=["R", "Z"])
@property
def volume(self):
"""
Volume of the closed fluxsurface calculated from the area
using Pappus centroid theorem : https://en.wikipedia.org/wiki/Pappus%27s_centroid_theorem
:return: float
"""
if self._closed:
return self._poly.area * 2 * np.pi * self.centroid.R[0]
else:
raise Exception("Opened Flux Surface does not have area")
@property
def diff_volume(self):
"""
Diferential volume :math:`V' = dV/d\psi`
Jardin, S.: Computational Methods in Plasma Physics
:return:
"""
if not hasattr(self, '_diff_volume'):
self._diff_volume = self._eval_diff_vol()
return self._diff_volume
def _eval_diff_vol(self):
Rs = (self.R[1:] + self.R[:-1]) / 2
Zs = (self.Z[1:] + self.Z[:-1]) / 2
dpsi = self._eq.diff_psi(Rs, Zs)
dl = self.dists
return 2 * np.pi * np.sum(Rs * dl / dpsi)
[docs]class FluxSurface(Surface):
[docs] def __init__(self, equilibrium, *coordinates, coord_type=None, grid=False, **coords):
"""
Calculates geometrical properties of the flux surface. To make the contour closed, the first and last points in
the passed coordinates have to be the same.
Instance is obtained by calling method `flux_surface` in instance of `Equilibrium`.
:param coords: Instance of coordinate class
"""
# FluxSurface can use only equilibrium default inner cocos if exists:
if equilibrium is not None:
super().__init__(equilibrium, *coordinates, coord_type=None, grid=False, cocos=equilibrium.cocos, **coords)
else:
super().__init__(equilibrium, *coordinates, coord_type=None, grid=False, **coords)
@property
def eval_q(self):
if not hasattr(self, '_q'):
self._q = self.get_eval_q('sum')
return self._q
[docs] def get_eval_q(self, method):
"""
Evaluete q usiong formula (5.35) from [Jardin, 2010: Computational methods in Plasma Physics]
:param method: str, ['sum', 'trapz', 'simps']
:return:
"""
# psi_sign = self._eq._psi_sign
cc = self.cocos_dict['sigma_Bp'] * self.cocos_dict['sigma_pol']
return cc * self._eq.F(psi_n=np.mean(self.psi_n), grid=False) / \
(2 * np.pi) ** (1 - self.cocos_dict['exp_Bp']) * \
self.surface_average(1 / self.R ** 2, method=method)
@property
def straight_fieldline_theta(self):
r"""
Calculate straight field line :math:`\theta^*` coordinate.
:return:
"""
from scipy.interpolate import CubicSpline
from pleque.utils.tools import arglis
if not hasattr(self, '_straight_fieldline_theta'):
h = 1 / self.R ** 2
izero = np.asscalar(np.argmin(np.mod(self.theta, 2 * np.pi)))
avg = self.surface_average(h)
cum_avg = self.cumsum_surface_average(h, roll=izero)
# ret = np.roll(ret, amin)
theta_star = 2 * np.pi * cum_avg / avg
# generate splines:
theta4spl = np.mod(self.theta, 2 * np.pi)
th_st4spl = theta_star
amin = np.asscalar(np.argmin(theta4spl))
theta4spl = np.roll(theta4spl, -amin)
th_st4spl = np.roll(th_st4spl, -amin)
# todo: this part is for testing
try:
self._theta2thetastar_spl = CubicSpline(theta4spl, th_st4spl, extrapolate='periodic')
except Exception:
# the sequence is not always strictly increasing.
# so get the increasing subsequence
alis = arglis(theta4spl)
theta4spl = theta4spl[alis]
th_st4spl = th_st4spl[alis]
self._theta2thetastar_spl = CubicSpline(theta4spl, th_st4spl, extrapolate='periodic')
try:
self._thetastar2theta_spl = CubicSpline(th_st4spl, theta4spl, extrapolate='periodic')
except Exception:
alis = arglis(th_st4spl)
theta4spl = theta4spl[alis]
th_st4spl = th_st4spl[alis]
self._thetastar2theta_spl = CubicSpline(th_st4spl, theta4spl, extrapolate='periodic')
self._straight_fieldline_theta = theta_star
return self._straight_fieldline_theta
@property
def tor_current(self):
"""
Return toroidal current through the closed flux surface
:return:
"""
if not hasattr(self, '_tor_current'):
self._tor_current = self._eval_tor_current()
return self._tor_current
def _eval_tor_current(self):
"""
Evalute magnitude of toroidal current through
:return:
"""
from scipy.constants import mu_0
# TODO: Precision of this method is quite poor. rel_err ~ 10e-5
# * Try improve precision by different integration techniques.
# * Increase fluxsurface position
# * Is it numeric error?
diff_psi = self._eq.diff_psi(self.R, self.Z)
cc_coef = self.cocos_dict['sigma_Bp'] / (2 * np.pi) ** self.cocos_dict['exp_Bp']
return cc_coef * 1 / mu_0 * self.surface_average(diff_psi ** 2 / self.R ** 2)
@property
@deprecated('Useless, will be removed. Use `abc` instead of `abc.contour`.')
def contour(self):
"""
Depracated. Fluxsurface contour points.
:return: numpy ndarray
"""
return self
[docs] def cumsum_surface_average(self, func, roll=0):
r"""
Return the surface average (over single magnetic surface) value of `func`.
Return the value of integration
.. math::
<func>(\psi)_i = \oint_0^{\theta_i} \frac{\mathrm{d}l R}{|\nabla \psi|}a(R, Z)
:param func: func(X, Y), Union[ndarray, int, float]
:return: ndarray
"""
import inspect
Rs = (self.R[1:] + self.R[:-1]) / 2
Zs = (self.Z[1:] + self.Z[:-1]) / 2
diff_psi = self._eq.diff_psi(Rs, Zs)
if inspect.isclass(func) or inspect.isfunction(func):
func_val = func(Rs, Zs)
elif isinstance(func, float) or isinstance(func, int):
func_val = func
else:
func_val = (func[1:] + func[:-1]) / 2
val = self.dists * Rs / diff_psi * func_val
val = np.roll(val, -roll)
ret = np.cumsum(val)
ret = np.hstack((0, ret))
ret = np.roll(ret, +roll)
return ret
[docs] def surface_average(self, func, method='sum'):
r"""
Return the surface average (over single magnetic surface) value of `func`.
Return the value of integration
.. math::
<func>(\psi) = \oint \frac{\mathrm{d}l R}{|\nabla \psi|}a(R, Z)
:param func: func(X, Y), Union[ndarray, int, float]
:param method: str, ['sum', 'trapz', 'simps']
:return:
"""
import inspect
from scipy.integrate import trapz, simps
if method == 'sum':
Rs = (self.R[1:] + self.R[:-1]) / 2
Zs = (self.Z[1:] + self.Z[:-1]) / 2
else:
Rs = self.R
Zs = self.Z
diff_psi = self._eq.diff_psi(Rs, Zs)
if inspect.isclass(func) or inspect.isfunction(func):
func_val = func(Rs, Zs)
elif isinstance(func, float) or isinstance(func, int):
func_val = func
else:
if method == 'sum':
func_val = (func[1:] + func[:-1]) / 2
else:
func_val = func
l = np.hstack((0, np.cumsum(self.dists)))
if method == 'sum':
ret = np.sum(self.dists * Rs / diff_psi * func_val)
elif method == 'trapz':
ret = trapz(Rs / diff_psi * func_val, l)
elif method == 'simps':
ret = simps(Rs / diff_psi * func_val, l)
else:
ret = None
return ret
[docs] def contains(self, coords: Coordinates):
points_RZ = coords.as_array(('R', 'Z'))[0, :]
if self._closed:
pnt = geometry.point.Point(points_RZ)
return self._poly.contains(pnt)
else:
raise Exception("Opened Flux Surface does not have area")
[docs] def distance(self, coords: Coordinates):
point = geometry.Point(coords.as_array()[0])
distance = self._string.distance(point)
return distance
@property
def max_radius(self):
'''
maximum radius on the given flux surface
:return:
'''
return np.max(self.R)
@property
def min_radius(self):
'''
minimum radius on the given flux surface
:return:
'''
return np.min(self.R)
@property
def minor_radius(self):
'''
a= (R_min - R_max)./2
:return:
'''
return (np.max(self.R) - np.min(self.R)) * 0.5
@property
def geom_radius(self):
'''
Geometrical radius
a= (R_min + R_max)./2
:return:
'''
return (np.max(self.R) + np.min(self.R)) * 0.5
@property
def elongation(self):
'''
Elongation
:return:
'''
return (np.max(self.Z) - np.min(self.Z)) * 0.5 / self.minor_radius
@property
def triangul_up(self):
'''
Upper triangularity
:return:
'''
imax = np.argmax(self.Z)
zma1 = self.R[imax]
rep = (self.geom_radius - zma1) / self.minor_radius
return rep
@property
def triangul_low(self):
'''
Lower triangularity
:return:
'''
imin = np.argmin(self.Z)
zma1 = self.R[imin]
rep= (self.geom_radius - zma1) / self.minor_radius
return rep
@property
def triangularity(self):
'''
triangularity
:return:
'''
rep = (self.triangul_up + self.triangul_low) * 0.5
return rep