from collections.abc import Iterable
import math
import numpy as np
[docs]def legendre_from_expcoef(coef, domain=(-1, 1)):
"""Return a Legendre series object based on expansion coefficients.
Given a list of coefficients from FET tally and a array of down, return
the numpy Legendre object.
Parameters
----------
coef : Iterable of float
A list of coefficients of each term in Legendre polynomials
domain : (2,) List of float
Domain of the Legendre polynomial
Returns
-------
numpy.polynomial.Legendre
A numpy Legendre series class
"""
n = np.arange(len(coef))
c = (2*n + 1) * np.asarray(coef) / (domain[1] - domain[0])
return np.polynomial.Legendre(c, domain)
class Polynomial:
"""Abstract Polynomial Class for creating polynomials.
"""
def __init__(self, coef):
self.coef = np.asarray(coef)
[docs]class ZernikeRadial(Polynomial):
"""Create radial only Zernike polynomials given coefficients and domain.
The radial only Zernike polynomials are defined as in
:class:`ZernikeRadialFilter`.
.. versionadded:: 0.12
Parameters
----------
coef : Iterable of float
A list of coefficients of each term in radial only Zernike polynomials
radius : float
Domain of Zernike polynomials to be applied on. Default is 1.
Attributes
----------
order : int
The maximum (even) order of Zernike polynomials.
radius : float
Domain of Zernike polynomials to be applied on. Default is 1.
norm_coef : iterable of float
The list of coefficients of each term in the polynomials after
normalization.
"""
def __init__(self, coef, radius=1):
super().__init__(coef)
self._order = 2 * (len(self.coef) - 1)
self.radius = radius
norm_vec = (2 * np.arange(len(self.coef)) + 1) / (math.pi * radius**2)
self._norm_coef = norm_vec * self.coef
@property
def order(self):
return self._order
def __call__(self, r):
import openmc.lib as lib
if isinstance(r, Iterable):
return [np.sum(self._norm_coef * lib.calc_zn_rad(self.order, r_i / self.radius))
for r_i in r]
else:
return np.sum(self._norm_coef * lib.calc_zn_rad(self.order, r / self.radius))
class Zernike(Polynomial):
r"""Create Zernike polynomials given coefficients and domain.
The azimuthal Zernike polynomials are defined as in :class:`ZernikeFilter`.
.. versionadded:: 0.12
Parameters
----------
coef : Iterable of float
A list of coefficients of each term in radial only Zernike polynomials
radius : float
Domain of Zernike polynomials to be applied on. Default is 1.
Attributes
----------
order : int
The maximum (even) order of Zernike polynomials.
radius : float
Domain of Zernike polynomials to be applied on. Default is 1.
theta : float
Azimuthal of Zernike polynomial to be applied on. Default is 0.
norm_coef : iterable of float
The list of coefficients of each term in the polynomials after
normalization.
"""
def __init__(self, coef, radius=1):
super().__init__(coef)
# Solve order from number of coefficients
# N = (order + 1)(order + 2) / 2
self._order = int((math.sqrt(8 * len(self.coef) + 1) - 3) / 2)
self.radius = radius
norm_vec = np.ones(len(self.coef))
for n in range(self._order + 1):
for m in range(-n, n + 1, 2):
j = int((n*(n + 2) + m)/2)
if m == 0:
norm_vec[j] = n + 1
else:
norm_vec[j] = 2*n + 2
norm_vec /= (math.pi * radius**2)
self._norm_coef = norm_vec * self.coef
@property
def order(self):
return self._order
def __call__(self, r, theta=0.0):
import openmc.lib as lib
if isinstance(r, Iterable) and isinstance(theta, Iterable):
return [[np.sum(self._norm_coef * lib.calc_zn(self.order, r_i / self.radius, theta_i))
for r_i in r] for theta_i in theta]
elif isinstance(r, Iterable) and not isinstance(theta, Iterable):
return [np.sum(self._norm_coef * lib.calc_zn(self.order, r_i / self.radius, theta))
for r_i in r]
elif not isinstance(r, Iterable) and isinstance(theta, Iterable):
return [np.sum(self._norm_coef * lib.calc_zn(self.order, r / self.radius, theta_i))
for theta_i in theta]
else:
return np.sum(self._norm_coef * lib.calc_zn(self.order, r / self.radius, theta))