diff --git a/openmc/capi/math.py b/openmc/capi/math.py index cd7ab40520..4de46691a6 100644 --- a/openmc/capi/math.py +++ b/openmc/capi/math.py @@ -11,15 +11,15 @@ _dll.t_percentile.restype = c_double _dll.t_percentile.argtypes = [POINTER(c_double), POINTER(c_int)] _dll.calc_pn.restype = c_double _dll.calc_pn.argtypes = [POINTER(c_int), POINTER(c_double)] +_dll.evaluate_legendre.restype = c_double +_dll.evaluate_legendre.argtypes = [POINTER(c_int), POINTER(c_double), + POINTER(c_double)] _dll.calc_rn.restype = None _dll.calc_rn.argtypes = [POINTER(c_int), ndpointer(c_double), ndpointer(c_double)] _dll.calc_zn.restype = None _dll.calc_zn.argtypes = [POINTER(c_int), POINTER(c_double), POINTER(c_double), ndpointer(c_double)] -_dll.evaluate_legendre.restype = c_double -_dll.evaluate_legendre.argtypes = [POINTER(c_int), POINTER(c_double), - POINTER(c_double)] _dll.rotate_angle.restype = None _dll.rotate_angle.argtypes = [ndpointer(c_double), POINTER(c_double), ndpointer(c_double), POINTER(c_double)] @@ -169,7 +169,7 @@ def calc_zn(n, rho, phi): """ - num_bins = ((n + 1) * (n + 2)) / 2 + num_bins = ((n + 1) * (n + 2)) // 2 zn = np.zeros(num_bins, dtype=np.float64) _dll.calc_zn(c_int(n), c_double(rho), c_double(phi), zn) return zn diff --git a/src/faddeeva/Faddeeva.h b/src/faddeeva/Faddeeva.h index 4293861907..9e26bc1ed9 100644 --- a/src/faddeeva/Faddeeva.h +++ b/src/faddeeva/Faddeeva.h @@ -1,5 +1,5 @@ /* Copyright (c) 2012 Massachusetts Institute of Technology - * + * * Permission is hereby granted, free of charge, to any person obtaining * a copy of this software and associated documentation files (the * "Software"), to deal in the Software without restriction, including @@ -7,17 +7,17 @@ * distribute, sublicense, and/or sell copies of the Software, and to * permit persons to whom the Software is furnished to do so, subject to * the following conditions: - * + * * The above copyright notice and this permission notice shall be * included in all copies or substantial portions of the Software. - * + * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ /* Available at: http://ab-initio.mit.edu/Faddeeva diff --git a/src/math.F90 b/src/math.F90 index d7e335284c..a386747a46 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -43,14 +43,6 @@ module math real(C_DOUBLE) :: pnx end function calc_pn_cc - subroutine calc_rn_cc(n, uvw, rn) bind(C, name='calc_rn_c') - use ISO_C_BINDING - implicit none - integer(C_INT), value, intent(in) :: n - real(C_DOUBLE), intent(in) :: uvw(3) - real(C_DOUBLE), intent(in) :: rn(2 * n + 1) - end subroutine calc_rn_cc - pure function evaluate_legendre_cc(n, data, x) & bind(C, name='evaluate_legendre_c') result(val) use ISO_C_BINDING @@ -61,6 +53,23 @@ module math real(C_DOUBLE) :: val end function evaluate_legendre_cc + subroutine calc_rn_cc(n, uvw, rn) bind(C, name='calc_rn_c') + use ISO_C_BINDING + implicit none + integer(C_INT), value, intent(in) :: n + real(C_DOUBLE), intent(in) :: uvw(3) + real(C_DOUBLE), intent(out) :: rn(2 * n + 1) + end subroutine calc_rn_cc + + subroutine calc_zn_cc(n, rho, phi, zn) bind(C, name='calc_zn_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(((n + 1) * (n + 2)) / 2) + end subroutine calc_zn_cc + subroutine rotate_angle_cc(uvw, mu, phi) bind(C, name='rotate_angle_c') use ISO_C_BINDING implicit none @@ -200,102 +209,13 @@ contains !=============================================================================== subroutine calc_zn(n, rho, phi, zn) bind(C) - ! 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. - integer(C_INT), intent(in) :: n ! Maximum order real(C_DOUBLE), intent(in) :: rho ! Radial location in the unit disk real(C_DOUBLE), intent(in) :: phi ! Theta (radians) location in the unit disk - real(C_DOUBLE), intent(out) :: zn(:) ! The resulting list of coefficients + ! The resulting list of coefficients + real(C_DOUBLE), intent(out) :: zn(((n + 1) * (n + 2)) / 2) - real(C_DOUBLE) :: sin_phi, cos_phi ! Sine and Cosine of phi - real(C_DOUBLE) :: sin_phi_vec(n+1) ! Contains sin(n*phi) - real(C_DOUBLE) :: cos_phi_vec(n+1) ! Contains cos(n*phi) - real(C_DOUBLE) :: zn_mat(n+1, n+1) ! Matrix form of the coefficients which is - ! easier to work with - real(C_DOUBLE) :: k1, k2, k3, k4 ! Variables for R_m_n calculation - integer(C_INT) :: i,p,q ! Loop counters - - real(C_DOUBLE), parameter :: SQRT_N_1(0:10) = [& - sqrt(1.0_8), sqrt(2.0_8), sqrt(3.0_8), sqrt(4.0_8), & - sqrt(5.0_8), sqrt(6.0_8), sqrt(7.0_8), sqrt(8.0_8), & - sqrt(9.0_8), sqrt(10.0_8), sqrt(11.0_8)] - real(C_DOUBLE), parameter :: SQRT_2N_2(0:10) = SQRT_N_1*sqrt(2.0_8) - - ! n == radial degree - ! m == azimuthal frequency - - ! ========================================================================== - ! Determine vector of sin(n*phi) and cos(n*phi). This takes advantage of the - ! following recurrence relations so that only a single sin/cos have to be - ! evaluated (http://mathworld.wolfram.com/Multiple-AngleFormulas.html) - ! - ! sin(nx) = 2 cos(x) sin((n-1)x) - sin((n-2)x) - ! cos(nx) = 2 cos(x) cos((n-1)x) - cos((n-2)x) - - sin_phi = sin(phi) - cos_phi = cos(phi) - - sin_phi_vec(1) = 1.0_8 - cos_phi_vec(1) = 1.0_8 - - sin_phi_vec(2) = 2.0_8 * cos_phi - cos_phi_vec(2) = cos_phi - - do i = 3, n+1 - sin_phi_vec(i) = 2.0_8 * cos_phi * sin_phi_vec(i-1) - sin_phi_vec(i-2) - cos_phi_vec(i) = 2.0_8 * cos_phi * cos_phi_vec(i-1) - cos_phi_vec(i-2) - end do - - do i = 1, n+1 - sin_phi_vec(i) = sin_phi_vec(i) * sin_phi - end do - - ! ========================================================================== - ! Calculate R_pq(rho) - - ! Fill the main diagonal first (Eq. 3.9 in Chong) - do p = 0, n - zn_mat(p+1, p+1) = rho**p - end do - - ! Fill in the second diagonal (Eq. 3.10 in Chong) - do q = 0, n-2 - zn_mat(q+2+1, q+1) = (q+2) * zn_mat(q+2+1, q+2+1) - (q+1) * zn_mat(q+1, q+1) - end do - - ! Fill in the rest of the values using the original results (Eq. 3.8 in Chong) - do p = 4, n - k2 = 2 * p * (p - 1) * (p - 2) - do q = p-4, 0, -2 - k1 = (p + q) * (p - q) * (p - 2) / 2 - k3 = -q**2*(p - 1) - p * (p - 1) * (p - 2) - k4 = -p * (p + q - 2) * (p - q - 2) / 2 - zn_mat(p+1, q+1) = ((k2 * rho**2 + k3) * zn_mat(p-2+1, q+1) + k4 * zn_mat(p-4+1, q+1)) / k1 - end do - end do - - ! Roll into a single vector for easier computation later - ! The vector is ordered (0,0), (1,-1), (1,1), (2,-2), (2,0), - ! (2, 2), .... in (n,m) indices - ! Note that the cos and sin vectors are offset by one - ! sin_phi_vec = [sin(x), sin(2x), sin(3x) ...] - ! cos_phi_vec = [1.0, cos(x), cos(2x)... ] - i = 1 - do p = 0, n - do q = -p, p, 2 - if (q < 0) then - zn(i) = zn_mat(p+1, abs(q)+1) * sin_phi_vec(abs(q)) * SQRT_2N_2(p) - else if (q == 0) then - zn(i) = zn_mat(p+1, q+1) * SQRT_N_1(p) - else - zn(i) = zn_mat(p+1, q+1) * cos_phi_vec(abs(q)+1) * SQRT_2N_2(p) - end if - i = i + 1 - end do - end do + call calc_zn_cc(n, rho, phi, zn) end subroutine calc_zn !=============================================================================== @@ -310,44 +230,11 @@ contains real(C_DOUBLE), intent(out) :: uvw(3) ! rotated directional cosine real(C_DOUBLE), optional :: phi ! azimuthal angle - real(C_DOUBLE) :: phi_ ! azimuthal angle - real(C_DOUBLE) :: sinphi ! sine of azimuthal angle - real(C_DOUBLE) :: cosphi ! cosine of azimuthal angle - real(C_DOUBLE) :: a ! sqrt(1 - mu^2) - real(C_DOUBLE) :: b ! sqrt(1 - w^2) - real(C_DOUBLE) :: u0 ! original cosine in x direction - real(C_DOUBLE) :: v0 ! original cosine in y direction - real(C_DOUBLE) :: w0 ! original cosine in z direction - - ! Copy original directional cosines - u0 = uvw0(1) - v0 = uvw0(2) - w0 = uvw0(3) - - ! Sample azimuthal angle in [0,2pi) if none provided + uvw = uvw0 if (present(phi)) then - phi_ = phi + call rotate_angle_cc(uvw, mu, phi) else - phi_ = TWO * PI * prn() - end if - - ! Precompute factors to save flops - sinphi = sin(phi_) - cosphi = cos(phi_) - a = sqrt(max(ZERO, ONE - mu*mu)) - b = sqrt(max(ZERO, ONE - w0*w0)) - - ! Need to treat special case where sqrt(1 - w**2) is close to zero by - ! expanding about the v component rather than the w component - if (b > 1e-10) then - uvw(1) = mu*u0 + a*(u0*w0*cosphi - v0*sinphi)/b - uvw(2) = mu*v0 + a*(v0*w0*cosphi + u0*sinphi)/b - uvw(3) = mu*w0 - a*b*cosphi - else - b = sqrt(ONE - v0*v0) - uvw(1) = mu*u0 + a*(u0*v0*cosphi + w0*sinphi)/b - uvw(2) = mu*v0 - a*b*cosphi - uvw(3) = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b + call rotate_angle_cc(uvw, mu, -10._8) end if end subroutine rotate_angle @@ -492,6 +379,9 @@ contains factors(i+3) = factors(i+1)*(E + (ONE + TWO * i) * half_inv_dopp2) end if end do + + ! call broaden_wmp_polynomials_cc(E, dopp, n, factors) + end subroutine broaden_wmp_polynomials end module math diff --git a/src/math_functions.cpp b/src/math_functions.cpp index 1466bde547..9f698f42a9 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -169,9 +169,10 @@ double __attribute__ ((const)) calc_pn_c(int n, double x) { double __attribute__ ((const)) evaluate_legendre_c(int n, double data[], double x) { double val; + int l; val = 0.5 * data[0]; - for (int l = 1; l < n; l++) { + for (l = 1; l < n; l++) { val += (static_cast(l) + 0.5) * data[l] * calc_pn_c(l, x); } return val; @@ -577,6 +578,122 @@ void calc_rn_c(int n, double uvw[3], double rn[]){ } } +//============================================================================== +// CALC_ZN calculates the n-th order modified Zernike polynomial moment for a +// given angle (rho, theta) location in the unit disk. The normalization of the +// polynomials is such tha the integral of Z_pq*Z_pq over the unit disk is +// exactly pi +//============================================================================== + +void calc_zn_c(int n, double rho, double phi, double zn[]) { + // 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. + double sin_phi; // Cosine of phi + double cos_phi; // Sine of phi + double sin_phi_vec[n + 1]; // Sin[n * phi] + double cos_phi_vec[n + 1]; // Cos[n * phi] + double zn_mat[n + 1][n + 1]; // Matrix forms of the coefficients which are + // easier to work with + // Variables for R_m_n calculation + double k1; + double k2; + double k3; + double k4; + // Loop counters + int i; + int p; + int q; + + const double SQRT_N_1[11] = + {std::sqrt(1.), std::sqrt(2.), std::sqrt(3.), std::sqrt(4.), + std::sqrt(5.), std::sqrt(6.), std::sqrt(7.), std::sqrt(8.), + std::sqrt(9.), std::sqrt(10.), std::sqrt(11.)}; + + const double SQRT_2N_2[11] = + {std::sqrt(2.) * std::sqrt(1.), std::sqrt(2.) * std::sqrt(2.), + std::sqrt(2.) * std::sqrt(3.), std::sqrt(2.) * std::sqrt(4.), + std::sqrt(2.) * std::sqrt(5.), std::sqrt(2.) * std::sqrt(6.), + std::sqrt(2.) * std::sqrt(7.), std::sqrt(2.) * std::sqrt(8.), + std::sqrt(2.) * std::sqrt(9.), std::sqrt(2.) * std::sqrt(10.), + std::sqrt(2.) * std::sqrt(11.)}; + + // n == radial degree + // m == azimuthal frequency + + // =========================================================================== + // Determine vector of sin(n*phi) and cos(n*phi). This takes advantage of the + // following recurrence relations so that only a single sin/cos have to be + // evaluated (http://mathworld.wolfram.com/Multiple-AngleFormulas.html) + // + // sin(nx) = 2 cos(x) sin((n-1)x) - sin((n-2)x) + // cos(nx) = 2 cos(x) cos((n-1)x) - cos((n-2)x) + + sin_phi = std::sin(phi); + cos_phi = std::cos(phi); + + sin_phi_vec[0] = 1.0; + cos_phi_vec[0] = 1.0; + + sin_phi_vec[1] = 2.0 * cos_phi; + cos_phi_vec[1] = cos_phi; + + for (i = 2; i <= n; i++) { + sin_phi_vec[i] = 2. * cos_phi * sin_phi_vec[i - 1] - sin_phi_vec[i - 2]; + cos_phi_vec[i] = 2. * cos_phi * cos_phi_vec[i - 1] - cos_phi_vec[i - 2]; + } + + for (i = 0; i <= n; i++) { + sin_phi_vec[i] *= sin_phi; + } + + // =========================================================================== + // Calculate R_pq(rho) + + // Fill the main diagonal first (Eq 3.9 in Chong) + for (p = 0; p <= n; p++) { + zn_mat[p][p] = std::pow(rho, p); + } + + // Fill the 2nd diagonal (Eq 3.10 in Chong) + for (q = 0; q <= n - 2; q++) { + zn_mat[q][q+2] = (q + 2) * zn_mat[q+2][q+2] - (q + 1) * zn_mat[q][q]; + } + + // Fill in the rest of the values using the original results (Eq. 3.8 in Chong) + for (p = 4; p <= n; p++) { + k2 = static_cast (2 * p * (p - 1) * (p - 2)); + for (q = p - 4; q >= 0; q -= 2) { + k1 = static_cast((p + q) * (p - q) * (p - 2)) / 2.; + k3 = static_cast(-q * q * (p - 1) - p * (p - 1) * (p - 2)); + k4 = static_cast(-p * (p + q - 2) * (p - q - 2)) / 2.; + zn_mat[q][p] = + ((k2 * rho * rho + k3) * zn_mat[q][p-2] + k4 * zn_mat[q][p-4]) / k1; + } + } + + // Roll into a single vector for easier computation later + // The vector is ordered (0,0), (1,-1), (1,1), (2,-2), (2,0), + // (2, 2), .... in (n,m) indices + // Note that the cos and sin vectors are offset by one + // sin_phi_vec = [sin(x), sin(2x), sin(3x) ...] + // cos_phi_vec = [1.0, cos(x), cos(2x)... ] + i = 0; + for (p = 0; p <= n; p++) { + for (q = -p; q <= p; q += 2) { + if (q < 0) { + zn[i] = zn_mat[std::abs(q)][p] * sin_phi_vec[std::abs(q) - 1] * SQRT_2N_2[p]; + } else if (q == 0) { + zn[i] = zn_mat[q][p] * SQRT_N_1[p]; + } else { + zn[i] = zn_mat[q][p] * cos_phi_vec[q] * SQRT_2N_2[p]; + } + i++; + } + } + +} //============================================================================== // ROTATE_ANGLE rotates direction std::cosines through a polar angle whose @@ -678,8 +795,8 @@ double watt_spectrum_c(double a, double b) { // FADDEEVA the Faddeeva function, using Stephen Johnson's implementation //============================================================================== -// std::complex faddeeva_c(std::complex z) { -// std::complex wv; // The resultant w(z) value +// double complex __attribute__ ((const)) faddeeva_c(double complex z) { +// double complex wv; // The resultant w(z) value // double relerr; // Target relative error in the inner loop of MIT Faddeeva // // Technically, the value we want is given by the equation: @@ -698,19 +815,19 @@ double watt_spectrum_c(double a, double b) { // // Note that faddeeva_w will interpret zero as machine epsilon // relerr = 0.; -// if (z.imag() > 0.) { -// wv = Faddeeva::w(z, relerr); +// if (cimag(z) > 0.) { +// wv = Faddeeva_w(z, relerr); // } else { -// wv = -std::conj(Faddeeva::w(std::conj(z), relerr)); +// wv = -conj(Faddeeva_w(conj(z), relerr)); // } // return wv; // } -// std::complex w_derivative_c(std::complex z, int order){ -// std::complex wv; // The resultant w(z) value +// double complex w_derivative_c(double complex z, int order){ +// double complex wv; // The resultant w(z) value -// const std::complex twoi_sqrtpi(0.0, 2.0 / std::sqrt(PI)); +// const double complex twoi_sqrtpi = 2.0 / std::sqrt(PI) * I; // switch(order) { // case 0: diff --git a/src/math_functions.h b/src/math_functions.h index aaba082615..4e7b2e8c64 100644 --- a/src/math_functions.h +++ b/src/math_functions.h @@ -2,11 +2,11 @@ #define MATH_FUNCTIONS_H #include -#include +#include #include #include "random_lcg.h" -// #include "faddeeva/Faddeeva.hh" +// #include "faddeeva/Faddeeva.h" namespace openmc { @@ -51,6 +51,15 @@ extern "C" double evaluate_legendre_c(int n, double data[], double x) extern "C" void calc_rn_c(int n, double uvw[3], double rn[]); +//============================================================================== +// CALC_ZN calculates the n-th order modified Zernike polynomial moment for a +// given angle (rho, theta) location in the unit disk. The normalization of the +// polynomials is such tha the integral of Z_pq*Z_pq over the unit disk is +// exactly pi +//============================================================================== + +extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]); + //============================================================================== // ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is // mu and through an azimuthal angle sampled uniformly. Note that this is done @@ -83,16 +92,17 @@ extern "C" double watt_spectrum_c(double a, double b); // FADDEEVA the Faddeeva function, using Stephen Johnson's implementation //============================================================================== -// extern "C" std::complex faddeeva_c(std::complex z); +// extern "C" double complex faddeeva_c(double complex z) __attribute__ ((const)); -// extern "C" std::complex w_derivative_c(std::complex z, int order); +// extern "C" double complex w_derivative_c(double complex z, int order); //============================================================================== // BROADEN_WMP_POLYNOMIALS Doppler broadens the windowed multipole curvefit. // The curvefit is a polynomial of the form a/E + b/sqrt(E) + c + d sqrt(E) ... //============================================================================== -extern "C" void broaden_wmp_polynomials_c(double E, double dopp, int n, double factors[]); +extern "C" void broaden_wmp_polynomials_c(double E, double dopp, int n, + double factors[]); } // namespace openmc #endif // MATH_FUNCTIONS_H \ No newline at end of file diff --git a/tests/unit_tests/test_math.py b/tests/unit_tests/test_math.py index bfa021cebd..03309bafe8 100644 --- a/tests/unit_tests/test_math.py +++ b/tests/unit_tests/test_math.py @@ -115,7 +115,38 @@ def test_calc_rn(): def test_calc_zn(): - pass + n = 10 + rho = 0.5 + phi = 0.5 + + # Reference solution from running the Fortran implementation + ref_vals = np.array( + [1.00000000e+00, 4.79425539e-01, 8.77582562e-01, + 5.15293637e-01, -8.66025404e-01, 3.30866239e-01, + 3.52667735e-01, -8.47512624e-01, -1.55136145e+00, + 2.50093775e-02, 1.79715684e-01, -1.33048245e+00, + -2.79508497e-01, -8.54292956e-01, -8.22482403e-02, + 6.47865100e-02, -1.18780200e+00, 5.18993370e-01, + 9.50010991e-01, -8.42327938e-02, -8.67263404e-02, + 8.25035501e-03, -7.44248626e-01, 1.52505281e+00, + 1.15751620e+00, 9.79225149e-01, 3.40611006e-01, + -5.78783240e-02, -1.09619759e-02, -3.17938327e-01, + 1.90147482e+00, 2.84658914e-01, 5.21064646e-01, + 1.34842791e-01, 4.25607546e-01, -2.92642715e-02, + -1.25423479e-02, -4.67751162e-02, 1.50696182e+00, + -6.13603897e-01, -8.67187500e-01, -3.93990531e-01, + -6.89672461e-01, 3.28139254e-01, -1.08327149e-02, + -8.53837419e-03, 7.04712042e-02, 7.73660979e-01, + -1.63799891e+00, -8.12396290e-01, -1.48708143e+00, + -1.16158437e-01, -1.03565982e+00, 1.88131088e-01, + -1.84122553e-03, -4.39233743e-03, 9.01295675e-02, + 1.32511582e-01, -1.83260987e+00, -7.93994967e-01, + -2.97978009e-01, -5.09818305e-01, 8.38707753e-01, + -9.29602211e-01, 7.78441102e-02, 1.29931014e-03]) + + test_vals = openmc.capi.math.calc_zn(n, rho, phi) + + assert np.allclose(ref_vals, test_vals) def test_rotate_angle(): @@ -128,7 +159,7 @@ def test_rotate_angle(): test_uvw = openmc.capi.math.rotate_angle(uvw0, mu, phi) - assert np.allclose(ref_uvw, test_uvw) + assert np.array_equal(ref_uvw, test_uvw) # Repeat for mu = 1 (no change) mu = 1. @@ -136,7 +167,7 @@ def test_rotate_angle(): test_uvw = openmc.capi.math.rotate_angle(uvw0, mu, phi) - assert np.allclose(ref_uvw, test_uvw) + assert np.array_equal(ref_uvw, test_uvw) # Need to test phi=None somehow... @@ -180,3 +211,5 @@ def test_broaden_wmp_polynomials(): test_val = openmc.capi.math.broaden_wmp_polynomials(test_E, test_dopp, n) assert np.allclose(ref_val, test_val) + +test_calc_zn() \ No newline at end of file