From 6f99f1ae18265f37924afbb84c7ded30946534c0 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Fri, 3 Aug 2018 09:38:39 -0500 Subject: [PATCH 01/24] Tally the particle only when it is within the range. --- src/tallies/tally_filter_zernike.F90 | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/tallies/tally_filter_zernike.F90 b/src/tallies/tally_filter_zernike.F90 index d005416138..ad4f592770 100644 --- a/src/tallies/tally_filter_zernike.F90 +++ b/src/tallies/tally_filter_zernike.F90 @@ -67,15 +67,17 @@ contains x = p % coord(1) % xyz(1) - this % x y = p % coord(1) % xyz(2) - this % y r = sqrt(x*x + y*y)/this % r - theta = atan2(y, x) + if (r <= 1) then + theta = atan2(y, x) + + ! Get moments for Zernike polynomial orders 0..n + call calc_zn(this % order, r, theta, zn) - ! Get moments for Zernike polynomial orders 0..n - call calc_zn(this % order, r, theta, zn) - - do i = 1, this % n_bins - call match % bins % push_back(i) - call match % weights % push_back(zn(i)) - end do + do i = 1, this % n_bins + call match % bins % push_back(i) + call match % weights % push_back(zn(i)) + end do + endif end subroutine get_all_bins subroutine to_statepoint(this, filter_group) From cdfa3565c9f95f793206a194563a1948a48e1572 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Fri, 3 Aug 2018 12:42:18 -0500 Subject: [PATCH 02/24] Add math function for radial Zernike --- src/math.F90 | 10 ++++++++++ src/math_functions.cpp | 29 +++++++++++++++++++++++++++++ src/math_functions.h | 23 +++++++++++++++++++++++ 3 files changed, 62 insertions(+) diff --git a/src/math.F90 b/src/math.F90 index 33f1ab4d30..87c747f579 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -11,6 +11,7 @@ module math public :: calc_pn public :: calc_rn public :: calc_zn + public :: calc_zn_rad public :: evaluate_legendre public :: rotate_angle public :: maxwell_spectrum @@ -70,6 +71,15 @@ module math real(C_DOUBLE), intent(out) :: zn(((n + 1) * (n + 2)) / 2) end subroutine calc_zn + pure subroutine calc_zn_rad(n, rho, phi, zn_rad) bind(C, name='calc_zn_rad_c') + use ISO_C_BINDING + implicit none + integer(C_INT), value, intent(in) :: n + real(C_DOUBLE), value, intent(in) :: rho + real(C_DOUBLE), value, intent(in) :: phi + real(C_DOUBLE), intent(out) :: zn_rad((n / 2) + 1) + end subroutine calc_zn + subroutine rotate_angle_c_intfc(uvw, mu, phi) bind(C, name='rotate_angle_c') use ISO_C_BINDING implicit none diff --git a/src/math_functions.cpp b/src/math_functions.cpp index 31f86ea5f6..2a79bb4b24 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -587,6 +587,35 @@ void calc_zn_c(int n, double rho, double phi, double zn[]) { } +void calc_zn_rad_c(int n, double rho, double phi, double zn_rad[]) { + // Calculate R_p0(rho) as Zn_p0(rho) + // Set up the array of the coefficients + int length = int(n/2) + 1; + double zn_rad[length]; + double q = 0; + + // R_00 is always 1 + zn_rad[0] = 1; + + // Fill in the rest of the array (Eq 3.8 and Eq 3.10 in Chong) + for (int p = 2; p <= n; p += 2) { + int index = int(p/2); + if (p == 2) { + // Setting up R_22 to calculate R_20 (Eq 3.10 in Chong) + R_22 = std::pow(rho, 2); + zn_rad[index] = 2 * R_22 - zn_rad[0] + } + 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); + double k4 = (-p * (p + q - 2) * (p - q - 2)) / 2.; + zn_rad[index] = + ((k2 * rho * rho + k3) * zn_rad[index-1] + k4 * zn_rad[index-2]) / k1; + } + } +} + void rotate_angle_c(double uvw[3], double mu, double* phi) { // Copy original directional cosines diff --git a/src/math_functions.h b/src/math_functions.h index 2ceea98d87..98882463df 100644 --- a/src/math_functions.h +++ b/src/math_functions.h @@ -91,6 +91,29 @@ extern "C" void calc_rn_c(int n, const double uvw[3], double rn[]); extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]); +//============================================================================== +//! Calculate only the even radial components of n-th order modified Zernike +//! polynomial moment with azimuthal dependency m = 0 for a given angle +//! (rho, theta) location on the unit disk. +//! +//! Since m = 0, n could only be even orders. Z_q0 = R_q0 +//! +//! This procedure uses the modified Kintner's method for calculating Zernike +//! polynomials as outlined in Chong, C. W., Raveendran, P., & Mukundan, +//! R. (2003). A comparative analysis of algorithms for fast computation of +//! Zernike moments. Pattern Recognition, 36(3), 731-742. +//! The normalization of the polynomials is such that the integral of Z_pq^2 +//! over the unit disk is exactly pi. +//! +//! @param n The maximum order requested +//! @param rho The radial parameter to specify location on the unit disk +//! @param phi The angle parameter to specify location on the unit disk +//! @param zn_rad The requested moments of order 0 to n (inclusive) +//! evaluated at rho and phi when m = 0. +//============================================================================== + +extern "C" void calc_zn_rad_c(int n, double rho, double phi, double zn_rad[]); + //============================================================================== //! Rotate the direction cosines through a polar angle whose cosine is mu and //! through an azimuthal angle sampled uniformly. From 1013782e125669bffcebb289707ab01a7dd8369a Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Fri, 3 Aug 2018 16:40:42 -0500 Subject: [PATCH 03/24] Add new filter file in CMakeLists --- CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 15c89142e0..878a41b935 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -377,6 +377,7 @@ add_library(libopenmc SHARED src/tallies/tally_filter_surface.F90 src/tallies/tally_filter_universe.F90 src/tallies/tally_filter_zernike.F90 + src/tallies/tally_filter_zernike_radial.F90 src/tallies/tally_header.F90 src/tallies/trigger.F90 src/tallies/trigger_header.F90 From dfa52096c51917741a1d1324f92e9aee326362b5 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Fri, 3 Aug 2018 16:42:14 -0500 Subject: [PATCH 04/24] Add new filter for radial zernike --- src/tallies/tally_filter_zernike_radial.F90 | 203 ++++++++++++++++++++ 1 file changed, 203 insertions(+) create mode 100644 src/tallies/tally_filter_zernike_radial.F90 diff --git a/src/tallies/tally_filter_zernike_radial.F90 b/src/tallies/tally_filter_zernike_radial.F90 new file mode 100644 index 0000000000..6eaf97f536 --- /dev/null +++ b/src/tallies/tally_filter_zernike_radial.F90 @@ -0,0 +1,203 @@ +module tally_filter_zernike_radial + + use, intrinsic :: ISO_C_BINDING + + use constants + use error + use hdf5_interface + use math, only: calc_zn_rad + use particle_header, only: Particle + use string, only: to_str + use tally_filter_header + use xml_interface + + implicit none + private + +!=============================================================================== +! ZERNIKERADIALFILTER gives even order radial Zernike polynomial moments of a +! particle's position +!=============================================================================== + + type, public, extends(TallyFilter) :: ZernikeRadialFilter + integer :: order + real(8) :: x + real(8) :: y + real(8) :: r + contains + procedure :: from_xml + procedure :: get_all_bins + procedure :: to_statepoint + procedure :: text_label + end type ZernikeRadialFilter + +contains + +!=============================================================================== +! ZernikeRadialFilter methods +!=============================================================================== + + subroutine from_xml(this, node) + class(ZernikeRadialFilter), intent(inout) :: this + type(XMLNode), intent(in) :: node + + integer :: n + + ! Get center of cylinder and radius + call get_node_value(node, "x", this % x) + call get_node_value(node, "y", this % y) + call get_node_value(node, "r", this % r) + + ! Get specified order + call get_node_value(node, "order", n) + this % order = n + this % n_bins = n/2 + 1 + end subroutine from_xml + + subroutine get_all_bins(this, p, estimator, match) + class(ZernikeRadialFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + type(TallyFilterMatch), intent(inout) :: match + + integer :: i + real(8) :: x, y, r + real(C_DOUBLE) :: zn_rad(this % n_bins) + + ! Determine normalized (r,theta) positions + x = p % coord(1) % xyz(1) - this % x + y = p % coord(1) % xyz(2) - this % y + r = sqrt(x*x + y*y)/this % r + if (r <= 1) then + + ! Get moments for even order Zernike polynomial orders 0..n + call calc_zn_rad(this % order, r, zn_rad) + + do i = 1, this % n_bins + call match % bins % push_back(i) + call match % weights % push_back(zn_rad(i)) + end do + endif + end subroutine get_all_bins + + subroutine to_statepoint(this, filter_group) + class(ZernikeRadialFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + call write_dataset(filter_group, "type", "zernikeradial") + call write_dataset(filter_group, "n_bins", this % n_bins) + call write_dataset(filter_group, "order", this % order) + call write_dataset(filter_group, "x", this % x) + call write_dataset(filter_group, "y", this % y) + call write_dataset(filter_group, "r", this % r) + end subroutine to_statepoint + + function text_label(this, bin) result(label) + class(ZernikeRadialFilter), intent(in) :: this + integer, intent(in) :: bin + 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 + +!=============================================================================== +! C API FUNCTIONS +!=============================================================================== + + function openmc_zernike_radial_filter_get_order(index, order) result(err) bind(C) + ! Get the order of an expansion filter + integer(C_INT32_T), value :: index + integer(C_INT), intent(out) :: order + integer(C_INT) :: err + + err = verify_filter(index) + if (err == 0) then + select type (f => filters(index) % obj) + type is (ZernikeRadialFilter) + order = f % order + class default + err = E_INVALID_TYPE + call set_errmsg("Not a Zernike Radial filter.") + end select + end if + end function openmc_zernike_radial_filter_get_order + + + function openmc_zernike_radial_filter_get_params(index, x, y, r) result(err) bind(C) + ! Get the Zernike filter parameters + integer(C_INT32_T), value :: index + real(C_DOUBLE), intent(out) :: x + real(C_DOUBLE), intent(out) :: y + real(C_DOUBLE), intent(out) :: r + integer(C_INT) :: err + + err = verify_filter(index) + if (err == 0) then + select type (f => filters(index) % obj) + type is (ZernikeRadialFilter) + x = f % x + y = f % y + r = f % r + class default + err = E_INVALID_TYPE + call set_errmsg("Not a Zernike Radial filter.") + end select + end if + end function openmc_zernike_radial_filter_get_params + + + function openmc_zernike_radial_filter_set_order(index, order) result(err) bind(C) + ! Set the order of an expansion filter + integer(C_INT32_T), value :: index + integer(C_INT), value :: order + integer(C_INT) :: err + + err = verify_filter(index) + if (err == 0) then + select type (f => filters(index) % obj) + type is (ZernikeRadialFilter) + f % order = order + f % n_bins = order/2 + 1 + class default + err = E_INVALID_TYPE + call set_errmsg("Not a Zernike Radial filter.") + end select + end if + end function openmc_zernike_radial_filter_set_order + + + function openmc_zernike_radial_filter_set_params(index, x, y, r) result(err) bind(C) + ! Set the Zernike filter parameters + integer(C_INT32_T), value :: index + real(C_DOUBLE), intent(in), optional :: x + real(C_DOUBLE), intent(in), optional :: y + real(C_DOUBLE), intent(in), optional :: r + integer(C_INT) :: err + + err = verify_filter(index) + if (err == 0) then + select type (f => filters(index) % obj) + type is (ZernikeRadialFilter) + if (present(x)) f % x = x + if (present(y)) f % y = y + if (present(r)) f % r = r + class default + err = E_INVALID_TYPE + call set_errmsg("Not a Zernike Radial filter.") + end select + end if + end function openmc_zernike_radial_filter_set_params + +end module tally_filter_zernike_radial From f8cac6004980586262732837997d8f8393e20ce9 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Fri, 3 Aug 2018 16:44:21 -0500 Subject: [PATCH 05/24] Add C API for the new radial zernike filters --- openmc/capi/filter.py | 30 ++++++++++++++++++++++++++++-- openmc/capi/math.py | 24 ++++++++++++++++++++++++ 2 files changed, 52 insertions(+), 2 deletions(-) diff --git a/openmc/capi/filter.py b/openmc/capi/filter.py index 391e63f073..621c2a4f6c 100644 --- a/openmc/capi/filter.py +++ b/openmc/capi/filter.py @@ -20,7 +20,7 @@ __all__ = ['Filter', 'AzimuthalFilter', 'CellFilter', 'EnergyFunctionFilter', 'LegendreFilter', 'MaterialFilter', 'MeshFilter', 'MeshSurfaceFilter', 'MuFilter', 'PolarFilter', 'SphericalHarmonicsFilter', 'SpatialLegendreFilter', 'SurfaceFilter', - 'UniverseFilter', 'ZernikeFilter', 'filters'] + 'UniverseFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters'] # Tally functions _dll.openmc_cell_filter_get_bins.argtypes = [ @@ -95,6 +95,12 @@ _dll.openmc_zernike_filter_get_order.errcheck = _error_handler _dll.openmc_zernike_filter_set_order.argtypes = [c_int32, c_int] _dll.openmc_zernike_filter_set_order.restype = c_int _dll.openmc_zernike_filter_set_order.errcheck = _error_handler +_dll.openmc_zernike_radial_filter_get_order.argtypes = [c_int32, POINTER(c_int)] +_dll.openmc_zernike_radial_filter_get_order.restype = c_int +_dll.openmc_zernike_radial_filter_get_order.errcheck = _error_handler +_dll.openmc_zernike_radial_filter_set_order.argtypes = [c_int32, c_int] +_dll.openmc_zernike_radial_filter_set_order.restype = c_int +_dll.openmc_zernike_radial_filter_set_order.errcheck = _error_handler class Filter(_FortranObjectWithID): __instances = WeakValueDictionary() @@ -360,6 +366,25 @@ class ZernikeFilter(Filter): _dll.openmc_zernike_filter_set_order(self._index, order) +class ZernikeRadialFilter(Filter): + filter_type = 'zernikeradial' + + def __init__(self, order=None, uid=None, new=True, index=None): + super().__init__(uid, new, index) + if order is not None: + self.order = order + + @property + def order(self): + temp_order = c_int() + _dll.openmc_zernike_radial_filter_get_order(self._index, temp_order) + return temp_order.value + + @order.setter + def order(self, order): + _dll.openmc_zernike_radial_filter_set_order(self._index, order) + + _FILTER_TYPE_MAP = { 'azimuthal': AzimuthalFilter, 'cell': CellFilter, @@ -380,7 +405,8 @@ _FILTER_TYPE_MAP = { 'spatiallegendre': SpatialLegendreFilter, 'surface': SurfaceFilter, 'universe': UniverseFilter, - 'zernike': ZernikeFilter + 'zernike': ZernikeFilter, + 'zernikeradial': ZernikeRadialFilter } diff --git a/openmc/capi/math.py b/openmc/capi/math.py index d1d7abdf64..4f8bb5cdbe 100644 --- a/openmc/capi/math.py +++ b/openmc/capi/math.py @@ -21,6 +21,9 @@ _dll.calc_rn_c.argtypes = [c_int, ndpointer(c_double), ndpointer(c_double)] _dll.calc_zn_c.restype = None _dll.calc_zn_c.argtypes = [c_int, c_double, c_double, ndpointer(c_double)] +_dll.calc_zn_rad_c.restype = None +_dll.calc_zn_rad_c.argtypes = [c_int, c_double, ndpointer(c_double)] + _dll.rotate_angle_c.restype = None _dll.rotate_angle_c.argtypes = [ndpointer(c_double), c_double, POINTER(c_double)] @@ -154,6 +157,27 @@ 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 + + Parameters + ---------- + n : int + Maximum order + rho : float + Radial location in the unit disk + + Returns + ------- + numpy.ndarray + Corresponding resulting list of coefficients + + """ + + num_bins = n // 2 + 1 + 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 From aa28f4a1452089e0e83ff2b82543d31f5f67ccb8 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Fri, 3 Aug 2018 16:46:00 -0500 Subject: [PATCH 06/24] Add new python class for radial zernike filter --- openmc/filter.py | 2 +- openmc/filter_expansion.py | 133 +++++++++++++++++++++++++++++++++++++ 2 files changed, 134 insertions(+), 1 deletion(-) diff --git a/openmc/filter.py b/openmc/filter.py index 3180989eb7..ae44024941 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -22,7 +22,7 @@ _FILTER_TYPES = ( 'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy', 'energyout', 'mu', 'polar', 'azimuthal', 'distribcell', 'delayedgroup', 'energyfunction', 'cellfrom', 'legendre', 'spatiallegendre', - 'sphericalharmonics', 'zernike' + 'sphericalharmonics', 'zernike', 'zernikeradial' ) _CURRENT_NAMES = ( diff --git a/openmc/filter_expansion.py b/openmc/filter_expansion.py index 1970c07710..842d03477f 100644 --- a/openmc/filter_expansion.py +++ b/openmc/filter_expansion.py @@ -463,3 +463,136 @@ class ZernikeFilter(ExpansionFilter): subelement.text = str(self.r) return element + +class ZernikeRadialFilter(ExpansionFilter): + r"""Score radial Zernike expansion moments in space up to specified order. + + The Zernike polynomials are defined the same as in ZernikeFilter. + + .. math:: + Z_n^m(\rho, \theta) = R_n^m(\rho) \cos (m\theta), \quad m > 0 + + Z_n^{m}(\rho, \theta) = R_n^{m}(\rho) \sin (m\theta), \quad m < 0 + + Z_n^{m}(\rho, \theta) = R_n^{m}(\rho), \quad m = 0 + + where the radial polynomials are + + .. math:: + R_n^m(\rho) = \sum\limits_{k=0}^{(n-m)/2} \frac{(-1)^k (n-k)!}{k! ( + \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. + + 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 + ---------- + order : int + Maximum radial Zernike polynomial order + x : float + x-coordinate of center of circle for normalization + y : float + y-coordinate of center of circle for normalization + r : int or None + Radius of circle for normalization + + Attributes + ---------- + order : int + Maximum radial Zernike polynomial order + x : float + x-coordinate of center of circle for normalization + y : float + y-coordinate of center of circle for normalization + r : int or None + Radius of circle for normalization + id : int + Unique identifier for the filter + num_bins : int + The number of filter bins + + """ + + def __init__(self, order, x=0.0, y=0.0, r=1.0, filter_id=None): + super().__init__(order, filter_id) + self.x = x + self.y = y + self.r = r + + def __hash__(self): + string = type(self).__name__ + '\n' + string += '{: <16}=\t{}\n'.format('\tOrder', self.order) + return hash(string) + + def __repr__(self): + string = type(self).__name__ + '\n' + string += '{: <16}=\t{}\n'.format('\tOrder', self.order) + string += '{: <16}=\t{}\n'.format('\tID', self.id) + return string + + @ExpansionFilter.order.setter + def order(self, order): + ExpansionFilter.order.__set__(self, order) + self.bins = ['Z{},0'.format(n) for n in range(0, order+1, 2)] + + @property + def x(self): + return self._x + + @x.setter + def x(self, x): + cv.check_type('x', x, Real) + self._x = x + + @property + def y(self): + return self._y + + @y.setter + def y(self, y): + cv.check_type('y', y, Real) + self._y = y + + @property + def r(self): + return self._r + + @r.setter + def r(self, r): + cv.check_type('r', r, Real) + self._r = r + + @classmethod + def from_hdf5(cls, group, **kwargs): + if group['type'].value.decode() != cls.short_name.lower(): + raise ValueError("Expected HDF5 data for filter type '" + + cls.short_name.lower() + "' but got '" + + group['type'].value.decode() + " instead") + + filter_id = int(group.name.split('/')[-1].lstrip('filter ')) + order = group['order'].value + x, y, r = group['x'].value, group['y'].value, group['r'].value + + return cls(order, x, y, r, filter_id) + + def to_xml_element(self): + """Return XML Element representing the filter. + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing radial Zernike filter data + + """ + element = super().to_xml_element() + subelement = ET.SubElement(element, 'x') + subelement.text = str(self.x) + subelement = ET.SubElement(element, 'y') + subelement.text = str(self.y) + subelement = ET.SubElement(element, 'r') + subelement.text = str(self.r) + + return element \ No newline at end of file From 76c4ba0c7823899e5c929c3caa7d15179de21d8b Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Fri, 3 Aug 2018 16:47:32 -0500 Subject: [PATCH 07/24] Math for radial components of zernike polynomials --- src/math.F90 | 5 ++--- src/math_functions.cpp | 9 ++++----- src/math_functions.h | 2 +- 3 files changed, 7 insertions(+), 9 deletions(-) diff --git a/src/math.F90 b/src/math.F90 index 87c747f579..ab5eff14c1 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -71,14 +71,13 @@ module math real(C_DOUBLE), intent(out) :: zn(((n + 1) * (n + 2)) / 2) end subroutine calc_zn - pure subroutine calc_zn_rad(n, rho, phi, zn_rad) bind(C, name='calc_zn_rad_c') + pure subroutine calc_zn_rad(n, rho, zn_rad) bind(C, name='calc_zn_rad_c') use ISO_C_BINDING implicit none integer(C_INT), value, intent(in) :: n real(C_DOUBLE), value, intent(in) :: rho - real(C_DOUBLE), value, intent(in) :: phi real(C_DOUBLE), intent(out) :: zn_rad((n / 2) + 1) - end subroutine calc_zn + end subroutine calc_zn_rad subroutine rotate_angle_c_intfc(uvw, mu, phi) bind(C, name='rotate_angle_c') use ISO_C_BINDING diff --git a/src/math_functions.cpp b/src/math_functions.cpp index 2a79bb4b24..53bf296dfe 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -587,11 +587,10 @@ void calc_zn_c(int n, double rho, double phi, double zn[]) { } -void calc_zn_rad_c(int n, double rho, double phi, double zn_rad[]) { +void calc_zn_rad_c(int n, double rho, double zn_rad[]) { // Calculate R_p0(rho) as Zn_p0(rho) // Set up the array of the coefficients - int length = int(n/2) + 1; - double zn_rad[length]; + double q = 0; // R_00 is always 1 @@ -602,8 +601,8 @@ void calc_zn_rad_c(int n, double rho, double phi, double zn_rad[]) { int index = int(p/2); if (p == 2) { // Setting up R_22 to calculate R_20 (Eq 3.10 in Chong) - R_22 = std::pow(rho, 2); - zn_rad[index] = 2 * R_22 - zn_rad[0] + double R_22 = std::pow(rho, 2); + zn_rad[index] = 2 * R_22 - zn_rad[0]; } else { double k1 = ((p + q) * (p - q) * (p - 2)) / 2.; diff --git a/src/math_functions.h b/src/math_functions.h index 98882463df..0993b1a01b 100644 --- a/src/math_functions.h +++ b/src/math_functions.h @@ -112,7 +112,7 @@ extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]); //! evaluated at rho and phi when m = 0. //============================================================================== -extern "C" void calc_zn_rad_c(int n, double rho, double phi, double zn_rad[]); +extern "C" void calc_zn_rad_c(int n, double rho, double zn_rad[]); //============================================================================== //! Rotate the direction cosines through a polar angle whose cosine is mu and From 0f9920719f1f1278ba340eb36839e8289c5c02bf Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Fri, 3 Aug 2018 16:48:10 -0500 Subject: [PATCH 08/24] Update the constant for radial zernike filter --- src/constants.F90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index ab6887610e..f377f55f33 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -343,7 +343,7 @@ module constants integer, parameter :: NO_BIN_FOUND = -1 ! Tally filter and map types - integer, parameter :: N_FILTER_TYPES = 20 + integer, parameter :: N_FILTER_TYPES = 21 integer, parameter :: & FILTER_UNIVERSE = 1, & FILTER_MATERIAL = 2, & @@ -364,7 +364,8 @@ module constants FILTER_LEGENDRE = 17, & FILTER_SPH_HARMONICS = 18, & FILTER_SPTL_LEGENDRE = 19, & - FILTER_ZERNIKE = 20 + FILTER_ZERNIKE = 20, & + FILTER_ZERNIKE_RADIAL = 21 ! Mesh types integer, parameter :: & From 8e9ed84e111fa01e827d8fc38fc8a96b0e4749a8 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Fri, 3 Aug 2018 16:49:51 -0500 Subject: [PATCH 09/24] Add the radial zernike filter in all the F90 files --- src/tallies/tally_filter.F90 | 5 +++++ src/tallies/tally_header.F90 | 3 +++ 2 files changed, 8 insertions(+) diff --git a/src/tallies/tally_filter.F90 b/src/tallies/tally_filter.F90 index 4b589cedbf..8d8bfadd57 100644 --- a/src/tallies/tally_filter.F90 +++ b/src/tallies/tally_filter.F90 @@ -26,6 +26,7 @@ module tally_filter use tally_filter_surface use tally_filter_universe use tally_filter_zernike + use tally_filter_zernike_radial implicit none @@ -88,6 +89,8 @@ contains type_ = 'universe' type is (ZernikeFilter) type_ = 'zernike' + type is (ZernikeRadialFilter) + type_ = 'zernikeradial' end select ! Convert Fortran string to null-terminated C string. We assume the @@ -159,6 +162,8 @@ contains allocate(UniverseFilter :: filters(index) % obj) case ('zernike') allocate(ZernikeFilter :: filters(index) % obj) + case ('zernikeradial') + allocate(ZernikeRadialFilter :: filters(index) % obj) case default err = E_UNASSIGNED call set_errmsg("Unknown filter type: " // trim(type_)) diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index ee44685924..224fafc3b8 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -339,6 +339,9 @@ contains type is (ZernikeFilter) j = FILTER_ZERNIKE this % estimator = ESTIMATOR_COLLISION + type is (ZernikeRadialFilter) + j = FILTER_ZERNIKE_RADIAL + this % estimator = ESTIMATOR_COLLISION end select this % find_filter(j) = i end do From 4c1c1c2cf00581e20a677df3a1e825d029d83690 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Mon, 6 Aug 2018 08:37:23 -0500 Subject: [PATCH 10/24] Add unit test for Radial Zernike filter --- tests/unit_tests/test_filters.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/tests/unit_tests/test_filters.py b/tests/unit_tests/test_filters.py index 488badbfaf..5817c4c3d2 100644 --- a/tests/unit_tests/test_filters.py +++ b/tests/unit_tests/test_filters.py @@ -100,6 +100,23 @@ def test_zernike(): assert elem.attrib['type'] == 'zernike' assert elem.find('order').text == str(n) +def test_zernike_radial(): + n = 4 + f = openmc.ZernikeRadialFilter(n, 0., 0., 1.) + assert f.order == n + assert f.bins[0] == 'Z0,0' + assert f.bins[-1] == 'Z{},0'.format(n) + assert len(f.bins) == n//2 + 1 + + # Make sure __repr__ works + repr(f) + + # to_xml_element() + elem = f.to_xml_element() + assert elem.tag == 'filter' + assert elem.attrib['type'] == 'zernikeradial' + assert elem.find('order').text == str(n) + def test_first_moment(run_in_tmpdir, box_model): plain_tally = openmc.Tally() From d65b7a11bd73426ef142056ada287841f8db0c4c Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Tue, 7 Aug 2018 10:19:29 -0500 Subject: [PATCH 11/24] Some travial modfication to meet the coding standard. --- openmc/capi/math.py | 7 ++++++- openmc/filter_expansion.py | 23 +++++++++++++-------- src/math_functions.cpp | 5 ++--- src/tallies/tally_filter_zernike_radial.F90 | 12 ----------- 4 files changed, 22 insertions(+), 25 deletions(-) diff --git a/openmc/capi/math.py b/openmc/capi/math.py index 4f8bb5cdbe..390ee466ff 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 842d03477f..a756f7ba4a 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 3d60ce4d0b..8a36f7b922 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 6eaf97f536..1a4ed5ed55 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 !=============================================================================== From 87209a8af217f0cf9035c35cac014428e5b2cf73 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Tue, 7 Aug 2018 11:22:11 -0500 Subject: [PATCH 12/24] ZernikeRadial as an extension of Zernike --- src/tallies/tally_filter_zernike.F90 | 119 ++++++++++-- src/tallies/tally_filter_zernike_radial.F90 | 191 -------------------- 2 files changed, 106 insertions(+), 204 deletions(-) delete mode 100644 src/tallies/tally_filter_zernike_radial.F90 diff --git a/src/tallies/tally_filter_zernike.F90 b/src/tallies/tally_filter_zernike.F90 index ad4f592770..d574c542d5 100644 --- a/src/tallies/tally_filter_zernike.F90 +++ b/src/tallies/tally_filter_zernike.F90 @@ -5,7 +5,7 @@ module tally_filter_zernike use constants use error use hdf5_interface - use math, only: calc_zn + use math, only: calc_zn, calc_zn_rad use particle_header, only: Particle use string, only: to_str use tally_filter_header @@ -24,19 +24,32 @@ module tally_filter_zernike real(8) :: y real(8) :: r contains - procedure :: from_xml - procedure :: get_all_bins - procedure :: to_statepoint - procedure :: text_label + procedure :: from_xml => from_xml_zn + procedure :: get_all_bins => get_all_bins_zn + procedure :: to_statepoint => to_statepoint_zn + procedure :: text_label => text_label_zn end type ZernikeFilter +!=============================================================================== +! ZERNIKERADIALFILTER gives even order radial Zernike polynomial moments of a +! particle's position +!=============================================================================== + + type, public, extends(ZernikeFilter) :: ZernikeRadialFilter + contains + procedure :: from_xml => from_xml_zn_rad + procedure :: get_all_bins => get_all_bins_zn_rad + procedure :: to_statepoint => to_statepoint_zn_rad + procedure :: text_label => text_label_zn_rad + end type ZernikeRadialFilter + contains !=============================================================================== ! ZernikeFilter methods !=============================================================================== - subroutine from_xml(this, node) + subroutine from_xml_zn(this, node) class(ZernikeFilter), intent(inout) :: this type(XMLNode), intent(in) :: node @@ -51,9 +64,9 @@ contains call get_node_value(node, "order", n) this % order = n this % n_bins = ((n + 1)*(n + 2))/2 - end subroutine from_xml + end subroutine from_xml_zn - subroutine get_all_bins(this, p, estimator, match) + subroutine get_all_bins_zn(this, p, estimator, match) class(ZernikeFilter), intent(in) :: this type(Particle), intent(in) :: p integer, intent(in) :: estimator @@ -78,9 +91,9 @@ contains call match % weights % push_back(zn(i)) end do endif - end subroutine get_all_bins + end subroutine get_all_bins_zn - subroutine to_statepoint(this, filter_group) + subroutine to_statepoint_zn(this, filter_group) class(ZernikeFilter), intent(in) :: this integer(HID_T), intent(in) :: filter_group @@ -90,9 +103,9 @@ contains call write_dataset(filter_group, "x", this % x) call write_dataset(filter_group, "y", this % y) call write_dataset(filter_group, "r", this % r) - end subroutine to_statepoint + end subroutine to_statepoint_zn - function text_label(this, bin) result(label) + function text_label_zn(this, bin) result(label) class(ZernikeFilter), intent(in) :: this integer, intent(in) :: bin character(MAX_LINE_LEN) :: label @@ -110,7 +123,74 @@ contains exit end if end do - end function text_label + end function text_label_zn + +!=============================================================================== +! ZernikeRadialFilter methods +!=============================================================================== + + subroutine from_xml_zn_rad(this, node) + class(ZernikeRadialFilter), intent(inout) :: this + type(XMLNode), intent(in) :: node + + integer :: n + + ! Get center of cylinder and radius + call get_node_value(node, "x", this % x) + call get_node_value(node, "y", this % y) + call get_node_value(node, "r", this % r) + + ! Get specified order + call get_node_value(node, "order", n) + this % order = n + this % n_bins = n/2 + 1 + end subroutine from_xml_zn_rad + + subroutine get_all_bins_zn_rad(this, p, estimator, match) + class(ZernikeRadialFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + type(TallyFilterMatch), intent(inout) :: match + + integer :: i + real(8) :: x, y, r + real(C_DOUBLE) :: zn_rad(this % n_bins) + + ! Determine normalized (r,theta) positions + x = p % coord(1) % xyz(1) - this % x + y = p % coord(1) % xyz(2) - this % y + r = sqrt(x*x + y*y)/this % r + if (r <= 1) then + + ! Get moments for even order Zernike polynomial orders 0..n + call calc_zn_rad(this % order, r, zn_rad) + + do i = 1, this % n_bins + call match % bins % push_back(i) + call match % weights % push_back(zn_rad(i)) + end do + endif + end subroutine get_all_bins_zn_rad + + subroutine to_statepoint_zn_rad(this, filter_group) + class(ZernikeRadialFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + call write_dataset(filter_group, "type", "zernikeradial") + call write_dataset(filter_group, "n_bins", this % n_bins) + call write_dataset(filter_group, "order", this % order) + call write_dataset(filter_group, "x", this % x) + call write_dataset(filter_group, "y", this % y) + call write_dataset(filter_group, "r", this % r) + end subroutine to_statepoint_zn_rad + + function text_label_zn_rad(this, bin) result(label) + class(ZernikeRadialFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + + label = "Zernike expansion, Z" // trim(to_str(2*bin)) // ",0" + end function text_label_zn_rad !=============================================================================== ! C API FUNCTIONS @@ -127,6 +207,8 @@ contains select type (f => filters(index) % obj) type is (ZernikeFilter) order = f % order + type is (ZernikeRadialFilter) + order = f % order class default err = E_INVALID_TYPE call set_errmsg("Not a Zernike filter.") @@ -150,6 +232,10 @@ contains x = f % x y = f % y r = f % r + type is (ZernikeRadialFilter) + x = f % x + y = f % y + r = f % r class default err = E_INVALID_TYPE call set_errmsg("Not a Zernike filter.") @@ -170,6 +256,9 @@ contains type is (ZernikeFilter) f % order = order f % n_bins = ((order + 1)*(order + 2))/2 + type is (ZernikeRadialFilter) + f % order = order + f % n_bins = order/2 + 1 class default err = E_INVALID_TYPE call set_errmsg("Not a Zernike filter.") @@ -193,6 +282,10 @@ contains if (present(x)) f % x = x if (present(y)) f % y = y if (present(r)) f % r = r + type is (ZernikeRadialFilter) + if (present(x)) f % x = x + if (present(y)) f % y = y + if (present(r)) f % r = r class default err = E_INVALID_TYPE call set_errmsg("Not a Zernike filter.") diff --git a/src/tallies/tally_filter_zernike_radial.F90 b/src/tallies/tally_filter_zernike_radial.F90 deleted file mode 100644 index 1a4ed5ed55..0000000000 --- a/src/tallies/tally_filter_zernike_radial.F90 +++ /dev/null @@ -1,191 +0,0 @@ -module tally_filter_zernike_radial - - use, intrinsic :: ISO_C_BINDING - - use constants - use error - use hdf5_interface - use math, only: calc_zn_rad - use particle_header, only: Particle - use string, only: to_str - use tally_filter_header - use xml_interface - - implicit none - private - -!=============================================================================== -! ZERNIKERADIALFILTER gives even order radial Zernike polynomial moments of a -! particle's position -!=============================================================================== - - type, public, extends(TallyFilter) :: ZernikeRadialFilter - integer :: order - real(8) :: x - real(8) :: y - real(8) :: r - contains - procedure :: from_xml - procedure :: get_all_bins - procedure :: to_statepoint - procedure :: text_label - end type ZernikeRadialFilter - -contains - -!=============================================================================== -! ZernikeRadialFilter methods -!=============================================================================== - - subroutine from_xml(this, node) - class(ZernikeRadialFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: n - - ! Get center of cylinder and radius - call get_node_value(node, "x", this % x) - call get_node_value(node, "y", this % y) - call get_node_value(node, "r", this % r) - - ! Get specified order - call get_node_value(node, "order", n) - this % order = n - this % n_bins = n/2 + 1 - end subroutine from_xml - - subroutine get_all_bins(this, p, estimator, match) - class(ZernikeRadialFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: i - real(8) :: x, y, r - real(C_DOUBLE) :: zn_rad(this % n_bins) - - ! Determine normalized (r,theta) positions - x = p % coord(1) % xyz(1) - this % x - y = p % coord(1) % xyz(2) - this % y - r = sqrt(x*x + y*y)/this % r - if (r <= 1) then - - ! Get moments for even order Zernike polynomial orders 0..n - call calc_zn_rad(this % order, r, zn_rad) - - do i = 1, this % n_bins - call match % bins % push_back(i) - call match % weights % push_back(zn_rad(i)) - end do - endif - end subroutine get_all_bins - - subroutine to_statepoint(this, filter_group) - class(ZernikeRadialFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "zernikeradial") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "order", this % order) - call write_dataset(filter_group, "x", this % x) - call write_dataset(filter_group, "y", this % y) - call write_dataset(filter_group, "r", this % r) - end subroutine to_statepoint - - function text_label(this, bin) result(label) - class(ZernikeRadialFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = "Zernike expansion, Z" // trim(to_str(2*bin)) // ",0" - end function text_label - -!=============================================================================== -! C API FUNCTIONS -!=============================================================================== - - function openmc_zernike_radial_filter_get_order(index, order) result(err) bind(C) - ! Get the order of an expansion filter - integer(C_INT32_T), value :: index - integer(C_INT), intent(out) :: order - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (ZernikeRadialFilter) - order = f % order - class default - err = E_INVALID_TYPE - call set_errmsg("Not a Zernike Radial filter.") - end select - end if - end function openmc_zernike_radial_filter_get_order - - - function openmc_zernike_radial_filter_get_params(index, x, y, r) result(err) bind(C) - ! Get the Zernike filter parameters - integer(C_INT32_T), value :: index - real(C_DOUBLE), intent(out) :: x - real(C_DOUBLE), intent(out) :: y - real(C_DOUBLE), intent(out) :: r - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (ZernikeRadialFilter) - x = f % x - y = f % y - r = f % r - class default - err = E_INVALID_TYPE - call set_errmsg("Not a Zernike Radial filter.") - end select - end if - end function openmc_zernike_radial_filter_get_params - - - function openmc_zernike_radial_filter_set_order(index, order) result(err) bind(C) - ! Set the order of an expansion filter - integer(C_INT32_T), value :: index - integer(C_INT), value :: order - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (ZernikeRadialFilter) - f % order = order - f % n_bins = order/2 + 1 - class default - err = E_INVALID_TYPE - call set_errmsg("Not a Zernike Radial filter.") - end select - end if - end function openmc_zernike_radial_filter_set_order - - - function openmc_zernike_radial_filter_set_params(index, x, y, r) result(err) bind(C) - ! Set the Zernike filter parameters - integer(C_INT32_T), value :: index - real(C_DOUBLE), intent(in), optional :: x - real(C_DOUBLE), intent(in), optional :: y - real(C_DOUBLE), intent(in), optional :: r - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (ZernikeRadialFilter) - if (present(x)) f % x = x - if (present(y)) f % y = y - if (present(r)) f % r = r - class default - err = E_INVALID_TYPE - call set_errmsg("Not a Zernike Radial filter.") - end select - end if - end function openmc_zernike_radial_filter_set_params - -end module tally_filter_zernike_radial From abdd753d5926c92ef1ea15226084a800bb56d1b7 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Tue, 7 Aug 2018 11:23:45 -0500 Subject: [PATCH 13/24] Delete unecessary names in files --- CMakeLists.txt | 1 - src/tallies/tally_filter.F90 | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 5f5706124f..73b3a3851b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -380,7 +380,6 @@ add_library(libopenmc SHARED src/tallies/tally_filter_surface.F90 src/tallies/tally_filter_universe.F90 src/tallies/tally_filter_zernike.F90 - src/tallies/tally_filter_zernike_radial.F90 src/tallies/tally_header.F90 src/tallies/trigger.F90 src/tallies/trigger_header.F90 diff --git a/src/tallies/tally_filter.F90 b/src/tallies/tally_filter.F90 index e2dcdd87d4..efeac180e2 100644 --- a/src/tallies/tally_filter.F90 +++ b/src/tallies/tally_filter.F90 @@ -27,7 +27,7 @@ module tally_filter use tally_filter_surface use tally_filter_universe use tally_filter_zernike - use tally_filter_zernike_radial + ! use tally_filter_zernike_radial implicit none From 0b45dc66c4874b3d55f2a19b756ecfec21d87f8f Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Tue, 7 Aug 2018 11:34:29 -0500 Subject: [PATCH 14/24] Change the ZernikeRadial filter on the Python API side as a subclass of Zernike filter. --- openmc/capi/filter.py | 21 ---------- openmc/filter_expansion.py | 78 +------------------------------------- 2 files changed, 1 insertion(+), 98 deletions(-) diff --git a/openmc/capi/filter.py b/openmc/capi/filter.py index 621c2a4f6c..5bd451b300 100644 --- a/openmc/capi/filter.py +++ b/openmc/capi/filter.py @@ -95,12 +95,6 @@ _dll.openmc_zernike_filter_get_order.errcheck = _error_handler _dll.openmc_zernike_filter_set_order.argtypes = [c_int32, c_int] _dll.openmc_zernike_filter_set_order.restype = c_int _dll.openmc_zernike_filter_set_order.errcheck = _error_handler -_dll.openmc_zernike_radial_filter_get_order.argtypes = [c_int32, POINTER(c_int)] -_dll.openmc_zernike_radial_filter_get_order.restype = c_int -_dll.openmc_zernike_radial_filter_get_order.errcheck = _error_handler -_dll.openmc_zernike_radial_filter_set_order.argtypes = [c_int32, c_int] -_dll.openmc_zernike_radial_filter_set_order.restype = c_int -_dll.openmc_zernike_radial_filter_set_order.errcheck = _error_handler class Filter(_FortranObjectWithID): __instances = WeakValueDictionary() @@ -369,21 +363,6 @@ class ZernikeFilter(Filter): class ZernikeRadialFilter(Filter): filter_type = 'zernikeradial' - def __init__(self, order=None, uid=None, new=True, index=None): - super().__init__(uid, new, index) - if order is not None: - self.order = order - - @property - def order(self): - temp_order = c_int() - _dll.openmc_zernike_radial_filter_get_order(self._index, temp_order) - return temp_order.value - - @order.setter - def order(self, order): - _dll.openmc_zernike_radial_filter_set_order(self._index, order) - _FILTER_TYPE_MAP = { 'azimuthal': AzimuthalFilter, diff --git a/openmc/filter_expansion.py b/openmc/filter_expansion.py index a756f7ba4a..ad761f2ff7 100644 --- a/openmc/filter_expansion.py +++ b/openmc/filter_expansion.py @@ -465,7 +465,7 @@ class ZernikeFilter(ExpansionFilter): return element -class ZernikeRadialFilter(ExpansionFilter): +class ZernikeRadialFilter(ZernikeFilter): r"""Score radial Zernike expansion moments in space up to specified order. The Zernike polynomials are defined the same as in ZernikeFilter. @@ -521,83 +521,7 @@ class ZernikeRadialFilter(ExpansionFilter): """ - def __init__(self, order, x=0.0, y=0.0, r=1.0, filter_id=None): - super().__init__(order, filter_id) - self.x = x - self.y = y - self.r = r - - def __hash__(self): - string = type(self).__name__ + '\n' - string += '{: <16}=\t{}\n'.format('\tOrder', self.order) - return hash(string) - - def __repr__(self): - string = type(self).__name__ + '\n' - string += '{: <16}=\t{}\n'.format('\tOrder', self.order) - string += '{: <16}=\t{}\n'.format('\tID', self.id) - return string - @ExpansionFilter.order.setter def order(self, order): ExpansionFilter.order.__set__(self, order) self.bins = ['Z{},0'.format(n) for n in range(0, order+1, 2)] - - @property - def x(self): - return self._x - - @x.setter - def x(self, x): - cv.check_type('x', x, Real) - self._x = x - - @property - def y(self): - return self._y - - @y.setter - def y(self, y): - cv.check_type('y', y, Real) - self._y = y - - @property - def r(self): - return self._r - - @r.setter - def r(self, r): - cv.check_type('r', r, Real) - self._r = r - - @classmethod - def from_hdf5(cls, group, **kwargs): - if group['type'].value.decode() != cls.short_name.lower(): - raise ValueError("Expected HDF5 data for filter type '" - + cls.short_name.lower() + "' but got '" - + group['type'].value.decode() + " instead") - - filter_id = int(group.name.split('/')[-1].lstrip('filter ')) - order = group['order'].value - x, y, r = group['x'].value, group['y'].value, group['r'].value - - return cls(order, x, y, r, filter_id) - - def to_xml_element(self): - """Return XML Element representing the filter. - - Returns - ------- - element : xml.etree.ElementTree.Element - XML element containing radial Zernike filter data - - """ - element = super().to_xml_element() - subelement = ET.SubElement(element, 'x') - subelement.text = str(self.x) - subelement = ET.SubElement(element, 'y') - subelement.text = str(self.y) - subelement = ET.SubElement(element, 'r') - subelement.text = str(self.r) - - return element \ No newline at end of file From a4cdef712d75ce2934b34bc120d99653c37dd32b Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Tue, 7 Aug 2018 11:39:44 -0500 Subject: [PATCH 15/24] Remove comment --- src/tallies/tally_filter.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/tallies/tally_filter.F90 b/src/tallies/tally_filter.F90 index efeac180e2..81b1958169 100644 --- a/src/tallies/tally_filter.F90 +++ b/src/tallies/tally_filter.F90 @@ -27,7 +27,6 @@ module tally_filter use tally_filter_surface use tally_filter_universe use tally_filter_zernike - ! use tally_filter_zernike_radial implicit none From 1074da87f8a6277276c0cd5abc34afed5b566cf5 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Wed, 8 Aug 2018 13:27:41 -0500 Subject: [PATCH 16/24] Made some change accordingly to the suggestions. --- openmc/filter.py | 2 +- src/constants.F90 | 5 +++-- src/math_functions.h | 13 ++++--------- 3 files changed, 8 insertions(+), 12 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index 36053ec46c..e63d72412f 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -22,7 +22,7 @@ _FILTER_TYPES = ( 'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy', 'energyout', 'mu', 'polar', 'azimuthal', 'distribcell', 'delayedgroup', 'energyfunction', 'cellfrom', 'legendre', 'spatiallegendre', - 'sphericalharmonics', 'zernike', 'particle', 'zernikeradial' + 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle' ) _CURRENT_NAMES = ( diff --git a/src/constants.F90 b/src/constants.F90 index 1bdcec74fb..4de52322e0 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -379,8 +379,9 @@ module constants FILTER_SPH_HARMONICS = 18, & FILTER_SPTL_LEGENDRE = 19, & FILTER_ZERNIKE = 20, & - FILTER_PARTICLE = 21, & - FILTER_ZERNIKE_RADIAL = 22 + FILTER_ZERNIKE_RADIAL = 21, & + FILTER_PARTICLE = 22 + ! Mesh types integer, parameter :: & diff --git a/src/math_functions.h b/src/math_functions.h index e0ecad60ec..7ef39c4461 100644 --- a/src/math_functions.h +++ b/src/math_functions.h @@ -92,18 +92,13 @@ extern "C" void calc_rn_c(int n, const double uvw[3], double rn[]); extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]); //============================================================================== -//! Calculate only the even radial components of n-th order modified Zernike -//! polynomial moment with azimuthal dependency m = 0 for a given angle -//! (rho, theta) location on the unit disk. +//! Calculate only the even order components of n-th order modified Zernike +//! polynomial moment with azimuthal dependency m = 0 for a given radial (rho) +//! location on the unit disk. //! //! Since m = 0, n could only be even orders. Z_q0 = R_q0 //! -//! This procedure uses the modified Kintner's method for calculating Zernike -//! polynomials as outlined in Chong, C. W., Raveendran, P., & Mukundan, -//! R. (2003). A comparative analysis of algorithms for fast computation of -//! Zernike moments. Pattern Recognition, 36(3), 731-742. -//! The normalization of the polynomials is such that the integral of Z_pq^2 -//! over the unit disk is exactly pi. +//! See calc_zn_c for methodology. //! //! @param n The maximum order requested //! @param rho The radial parameter to specify location on the unit disk From a7ce79c62f0cb9cd9f327126435b68c570547fb2 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Wed, 8 Aug 2018 13:32:08 -0500 Subject: [PATCH 17/24] ONE instead of 1 --- src/tallies/tally_filter_zernike.F90 | 16 ++++------------ 1 file changed, 4 insertions(+), 12 deletions(-) diff --git a/src/tallies/tally_filter_zernike.F90 b/src/tallies/tally_filter_zernike.F90 index d574c542d5..1ece62dbb2 100644 --- a/src/tallies/tally_filter_zernike.F90 +++ b/src/tallies/tally_filter_zernike.F90 @@ -80,7 +80,7 @@ contains x = p % coord(1) % xyz(1) - this % x y = p % coord(1) % xyz(2) - this % y r = sqrt(x*x + y*y)/this % r - if (r <= 1) then + if (r <= ONE) then theta = atan2(y, x) ! Get moments for Zernike polynomial orders 0..n @@ -160,7 +160,7 @@ contains x = p % coord(1) % xyz(1) - this % x y = p % coord(1) % xyz(2) - this % y r = sqrt(x*x + y*y)/this % r - if (r <= 1) then + if (r <= ONE) then ! Get moments for even order Zernike polynomial orders 0..n call calc_zn_rad(this % order, r, zn_rad) @@ -228,11 +228,7 @@ contains err = verify_filter(index) if (err == 0) then select type (f => filters(index) % obj) - type is (ZernikeFilter) - x = f % x - y = f % y - r = f % r - type is (ZernikeRadialFilter) + class is (ZernikeFilter) x = f % x y = f % y r = f % r @@ -278,11 +274,7 @@ contains err = verify_filter(index) if (err == 0) then select type (f => filters(index) % obj) - type is (ZernikeFilter) - if (present(x)) f % x = x - if (present(y)) f % y = y - if (present(r)) f % r = r - type is (ZernikeRadialFilter) + class is (ZernikeFilter) if (present(x)) f % x = x if (present(y)) f % y = y if (present(r)) f % r = r From 1bd1a1b855fccb26cecf7300032bcc4a5a90a092 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Wed, 8 Aug 2018 14:52:09 -0500 Subject: [PATCH 18/24] Remove some redundancy --- src/tallies/tally_filter_zernike.F90 | 39 ++++++++++++++++------------ 1 file changed, 23 insertions(+), 16 deletions(-) diff --git a/src/tallies/tally_filter_zernike.F90 b/src/tallies/tally_filter_zernike.F90 index 1ece62dbb2..96e016d2ee 100644 --- a/src/tallies/tally_filter_zernike.F90 +++ b/src/tallies/tally_filter_zernike.F90 @@ -13,6 +13,10 @@ module tally_filter_zernike implicit none private + public :: openmc_zernike_filter_get_order + public :: openmc_zernike_filter_get_params + public :: openmc_zernike_filter_set_order + public :: openmc_zernike_filter_set_params !=============================================================================== ! ZERNIKEFILTER gives Zernike polynomial moments of a particle's position @@ -24,6 +28,7 @@ module tally_filter_zernike real(8) :: y real(8) :: r contains + procedure :: calc_n_bins => calc_n_bins_zn procedure :: from_xml => from_xml_zn procedure :: get_all_bins => get_all_bins_zn procedure :: to_statepoint => to_statepoint_zn @@ -37,7 +42,8 @@ module tally_filter_zernike type, public, extends(ZernikeFilter) :: ZernikeRadialFilter contains - procedure :: from_xml => from_xml_zn_rad + procedure :: calc_n_bins => calc_n_bins_zn_rad + ! Inherit from_xml from ZernikeFilter procedure :: get_all_bins => get_all_bins_zn_rad procedure :: to_statepoint => to_statepoint_zn_rad procedure :: text_label => text_label_zn_rad @@ -49,6 +55,15 @@ contains ! ZernikeFilter methods !=============================================================================== + function calc_n_bins_zn(this) result(num_n_bins) + class(ZernikeFilter), intent(in) :: this + integer :: n + integer :: num_n_bins + + n = this % order + num_n_bins = ((n+1) * (n+2))/2 + end function calc_n_bins_zn + subroutine from_xml_zn(this, node) class(ZernikeFilter), intent(inout) :: this type(XMLNode), intent(in) :: node @@ -63,7 +78,7 @@ contains ! Get specified order call get_node_value(node, "order", n) this % order = n - this % n_bins = ((n + 1)*(n + 2))/2 + this % n_bins = this % calc_n_bins() end subroutine from_xml_zn subroutine get_all_bins_zn(this, p, estimator, match) @@ -129,22 +144,14 @@ contains ! ZernikeRadialFilter methods !=============================================================================== - subroutine from_xml_zn_rad(this, node) - class(ZernikeRadialFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - + function calc_n_bins_zn_rad(this) result(num_n_bins) + class(ZernikeRadialFilter), intent(in) :: this integer :: n + integer :: num_n_bins - ! Get center of cylinder and radius - call get_node_value(node, "x", this % x) - call get_node_value(node, "y", this % y) - call get_node_value(node, "r", this % r) - - ! Get specified order - call get_node_value(node, "order", n) - this % order = n - this % n_bins = n/2 + 1 - end subroutine from_xml_zn_rad + n = this % order + num_n_bins = n/2 + 1 + end function calc_n_bins_zn_rad subroutine get_all_bins_zn_rad(this, p, estimator, match) class(ZernikeRadialFilter), intent(in) :: this From fd6c8906aca1c3cc1600030e71277725757f20f1 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Thu, 9 Aug 2018 09:42:36 -0500 Subject: [PATCH 19/24] Change order --- src/tallies/tally_header.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index c6f3772e1f..7793739a49 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -339,11 +339,11 @@ contains type is (ZernikeFilter) j = FILTER_ZERNIKE this % estimator = ESTIMATOR_COLLISION - type is (ParticleFilter) - j = FILTER_PARTICLE type is (ZernikeRadialFilter) j = FILTER_ZERNIKE_RADIAL this % estimator = ESTIMATOR_COLLISION + type is (ParticleFilter) + j = FILTER_PARTICLE end select this % find_filter(j) = i end do From 960211de0cee2a04ce538e8dd87dbdc56f4e67a2 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Thu, 9 Aug 2018 11:05:16 -0500 Subject: [PATCH 20/24] Add ZernikeRadialFilter in Sphinx --- docs/source/pythonapi/base.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index b9eb7033f1..752cf417e8 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -122,6 +122,7 @@ Constructing Tallies openmc.SpatialLegendreFilter openmc.SphericalHarmonicsFilter openmc.ZernikeFilter + openmc.ZernikeRadialFilter openmc.ParticleFilter openmc.Mesh openmc.Trigger From 11723c8263e17161a11e3d88645d2710f4c6c36a Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Thu, 9 Aug 2018 11:15:22 -0500 Subject: [PATCH 21/24] Fortran index starts from 1. --- src/tallies/tally_filter_zernike.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tallies/tally_filter_zernike.F90 b/src/tallies/tally_filter_zernike.F90 index 96e016d2ee..b680fc06d0 100644 --- a/src/tallies/tally_filter_zernike.F90 +++ b/src/tallies/tally_filter_zernike.F90 @@ -133,7 +133,7 @@ contains if (bin <= last) then first = last - n m = -n + (bin - first)*2 - label = "Zernike expansion, Z" // trim(to_str(n)) // "," & + label = "Zernike expansion, Z" // trim(to_str((n-1))) // "," & // trim(to_str(m)) exit end if From 7f015db3434f909c0d0e62d06d0f543bfc7eb1a8 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Thu, 9 Aug 2018 12:41:30 -0500 Subject: [PATCH 22/24] Change the description of ZernikeRadialFilter --- openmc/filter_expansion.py | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/openmc/filter_expansion.py b/openmc/filter_expansion.py index ad761f2ff7..7debebd6b9 100644 --- a/openmc/filter_expansion.py +++ b/openmc/filter_expansion.py @@ -466,32 +466,30 @@ class ZernikeFilter(ExpansionFilter): class ZernikeRadialFilter(ZernikeFilter): - r"""Score radial Zernike expansion moments in space up to specified order. + r"""Score the :math:`m = 0` (radial variation only) Zernike moments up to + specified order. - The Zernike polynomials are defined the same as in ZernikeFilter. + The Zernike polynomials are defined the same as in :class:`ZernikeFilter`. .. math:: - Z_n^m(\rho, \theta) = R_n^m(\rho) \cos (m\theta), \quad m > 0 - Z_n^{m}(\rho, \theta) = R_n^{m}(\rho) \sin (m\theta), \quad m < 0 - - Z_n^{m}(\rho, \theta) = R_n^{m}(\rho), \quad m = 0 + Z_n^{0}(\rho, \theta) = R_n^{0}(\rho) where the radial polynomials are .. math:: - R_n^m(\rho) = \sum\limits_{k=0}^{(n-m)/2} \frac{(-1)^k (n-k)!}{k! ( - \frac{n+m}{2} - k)! (\frac{n-m}{2} - k)!} \rho^{n-2k}. + R_n^{0}(\rho) = \sum\limits_{k=0}^{n/2} \frac{(-1)^k (n-k)!}{k! (( + \frac{n}{2} - k)!)^{2}} \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. + With this definition, the integral of :math:`(Z_n^0)^2` over the unit disk + is :math:`\frac{\pi}{n+1}`. 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. + there are :math:`\frac{n}{2} + 1` terms in total. The indexing is from the + lowest even order (0) to highest even order. Parameters ---------- From e2774801a5f509b30e539cb7f257d583aec1480b Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Fri, 10 Aug 2018 09:22:08 -0500 Subject: [PATCH 23/24] Fixed the typo --- src/tallies/tally_filter_zernike.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/tallies/tally_filter_zernike.F90 b/src/tallies/tally_filter_zernike.F90 index b680fc06d0..73fda7fade 100644 --- a/src/tallies/tally_filter_zernike.F90 +++ b/src/tallies/tally_filter_zernike.F90 @@ -133,7 +133,7 @@ contains if (bin <= last) then first = last - n m = -n + (bin - first)*2 - label = "Zernike expansion, Z" // trim(to_str((n-1))) // "," & + label = "Zernike expansion, Z" // trim(to_str(n)) // "," & // trim(to_str(m)) exit end if @@ -196,7 +196,7 @@ contains integer, intent(in) :: bin character(MAX_LINE_LEN) :: label - label = "Zernike expansion, Z" // trim(to_str(2*bin)) // ",0" + label = "Zernike expansion, Z" // trim(to_str(2*(bin-1))) // ",0" end function text_label_zn_rad !=============================================================================== From daf57e663ac10773274cb95b193a058b7fa08548 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Mon, 13 Aug 2018 08:58:39 -0500 Subject: [PATCH 24/24] Fixed some issues --- openmc/capi/filter.py | 2 +- src/tallies/tally_filter_zernike.F90 | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/openmc/capi/filter.py b/openmc/capi/filter.py index 5bd451b300..8d9d794d16 100644 --- a/openmc/capi/filter.py +++ b/openmc/capi/filter.py @@ -360,7 +360,7 @@ class ZernikeFilter(Filter): _dll.openmc_zernike_filter_set_order(self._index, order) -class ZernikeRadialFilter(Filter): +class ZernikeRadialFilter(ZernikeFilter): filter_type = 'zernikeradial' diff --git a/src/tallies/tally_filter_zernike.F90 b/src/tallies/tally_filter_zernike.F90 index 73fda7fade..04515cf318 100644 --- a/src/tallies/tally_filter_zernike.F90 +++ b/src/tallies/tally_filter_zernike.F90 @@ -55,13 +55,13 @@ contains ! ZernikeFilter methods !=============================================================================== - function calc_n_bins_zn(this) result(num_n_bins) + function calc_n_bins_zn(this) result(n_bins) class(ZernikeFilter), intent(in) :: this integer :: n - integer :: num_n_bins + integer :: n_bins n = this % order - num_n_bins = ((n+1) * (n+2))/2 + n_bins = ((n+1) * (n+2))/2 end function calc_n_bins_zn subroutine from_xml_zn(this, node)