diff --git a/openmc/capi/math.py b/openmc/capi/math.py index 4f8bb5cdb..390ee466f 100644 --- a/openmc/capi/math.py +++ b/openmc/capi/math.py @@ -157,8 +157,12 @@ def calc_zn(n, rho, phi): _dll.calc_zn_c(n, rho, phi, zn) return zn + def calc_zn_rad(n, rho): - """ Calculate the even orders in n-th order modified Zernike polynomial moment with no azimuthal dependency (m=0) for a given radial location in the unit disk. The normalization of the polynomials is such that the integral of Z_pq*Z_pq over the unit disk is exactly pi + """ Calculate the even orders in n-th order modified Zernike polynomial + moment with no azimuthal dependency (m=0) for a given radial location in + the unit disk. The normalization of the polynomials is such that the + integral of Z_pq*Z_pq over the unit disk is exactly pi. Parameters ---------- @@ -178,6 +182,7 @@ def calc_zn_rad(n, rho): zn_rad = np.zeros(num_bins, dtype=np.float64) _dll.calc_zn_rad_c(n, rho, zn_rad) return zn_rad + def rotate_angle(uvw0, mu, phi=None): """ Rotates direction cosines through a polar angle whose cosine is diff --git a/openmc/filter_expansion.py b/openmc/filter_expansion.py index 842d03477..a756f7ba4 100644 --- a/openmc/filter_expansion.py +++ b/openmc/filter_expansion.py @@ -234,8 +234,8 @@ class SphericalHarmonicsFilter(ExpansionFilter): r"""Score spherical harmonic expansion moments up to specified order. This filter allows you to obtain real spherical harmonic moments of either - the particle's direction or the cosine of the scattering angle. Specifying a - filter with order :math:`\ell` tallies moments for all orders from 0 to + the particle's direction or the cosine of the scattering angle. Specifying + a filter with order :math:`\ell` tallies moments for all orders from 0 to :math:`\ell`. Parameters @@ -342,11 +342,11 @@ class ZernikeFilter(ExpansionFilter): \frac{n+m}{2} - k)! (\frac{n-m}{2} - k)!} \rho^{n-2k}. With this definition, the integral of :math:`(Z_n^m)^2` over the unit disk - is :math:`\frac{\epsilon_m\pi}{2n+2}` for each polynomial where :math:`\epsilon_m` is - 2 if :math:`m` equals 0 and 1 otherwise. + is :math:`\frac{\epsilon_m\pi}{2n+2}` for each polynomial where + :math:`\epsilon_m` is 2 if :math:`m` equals 0 and 1 otherwise. - Specifying a filter with order N tallies moments for all :math:`n` from 0 to - N and each value of :math:`m`. The ordering of the Zernike polynomial + Specifying a filter with order N tallies moments for all :math:`n` from 0 + to N and each value of :math:`m`. The ordering of the Zernike polynomial moments follows the ANSI Z80.28 standard, where the one-dimensional index :math:`j` corresponds to the :math:`n` and :math:`m` by @@ -464,6 +464,7 @@ class ZernikeFilter(ExpansionFilter): return element + class ZernikeRadialFilter(ExpansionFilter): r"""Score radial Zernike expansion moments in space up to specified order. @@ -483,10 +484,14 @@ class ZernikeRadialFilter(ExpansionFilter): \frac{n+m}{2} - k)! (\frac{n-m}{2} - k)!} \rho^{n-2k}. With this definition, the integral of :math:`(Z_n^m)^2` over the unit disk - is :math:`\frac{\epsilon_m\pi}{2n+2}` for each polynomial where :math:`\epsilon_m` is - 2 if :math:`m` equals 0 and 1 otherwise. + is :math:`\frac{\epsilon_m\pi}{2n+2}` for each polynomial where + :math:`\epsilon_m` is 2 if :math:`m` equals 0 and 1 otherwise. - If there is only radial dependency, the polynomials are integrated over the azimuthal angles. The only terms left are :math:`Z_n^{0}(\rho, \theta) = R_n^{0}(\rho)`. Note that :math:`n` could only be even orders. Therefore, for a radial Zernike polynomials up to order of :math:`n`, there are :math:`\frac{n}{2} + 1` terms in total. + If there is only radial dependency, the polynomials are integrated over + the azimuthal angles. The only terms left are :math:`Z_n^{0}(\rho, \theta) + = R_n^{0}(\rho)`. Note that :math:`n` could only be even orders. + Therefore, for a radial Zernike polynomials up to order of :math:`n`, + there are :math:`\frac{n}{2} + 1` terms in total. Parameters ---------- diff --git a/src/math_functions.cpp b/src/math_functions.cpp index 3d60ce4d0..8a36f7b92 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -601,10 +601,9 @@ void calc_zn_rad_c(int n, double rho, double zn_rad[]) { int index = int(p/2); if (p == 2) { // Setting up R_22 to calculate R_20 (Eq 3.10 in Chong) - double R_22 = std::pow(rho, 2); + double R_22 = rho * rho; zn_rad[index] = 2 * R_22 - zn_rad[0]; - } - else { + } else { double k1 = ((p + q) * (p - q) * (p - 2)) / 2.; double k2 = 2 * p * (p - 1) * (p - 2); double k3 = -q * q * (p - 1) - p * (p - 1) * (p - 2); diff --git a/src/tallies/tally_filter_zernike_radial.F90 b/src/tallies/tally_filter_zernike_radial.F90 index 6eaf97f53..1a4ed5ed5 100644 --- a/src/tallies/tally_filter_zernike_radial.F90 +++ b/src/tallies/tally_filter_zernike_radial.F90 @@ -98,18 +98,6 @@ contains character(MAX_LINE_LEN) :: label label = "Zernike expansion, Z" // trim(to_str(2*bin)) // ",0" - ! integer :: n - ! integer :: first, last - - ! do n = 0, this % order, -2 - ! last = (n + 1)*(n + 2)/2 - ! if (bin <= last) then - ! first = last - n - ! m = -n + (bin - first)*2 - ! label = "Zernike expansion, Z" // trim(to_str(n)) // ",0" - ! exit - ! end if - ! end do end function text_label !===============================================================================