diff --git a/CMakeLists.txt b/CMakeLists.txt index 7ecd6e321..dc5558de1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -435,6 +435,7 @@ set(LIBOPENMC_CXX_SRC src/initialize.cpp src/finalize.cpp src/hdf5_interface.cpp + src/math_functions.cpp src/message_passing.cpp src/plot.cpp src/random_lcg.cpp diff --git a/openmc/capi/__init__.py b/openmc/capi/__init__.py index 217e782a8..672a6e811 100644 --- a/openmc/capi/__init__.py +++ b/openmc/capi/__init__.py @@ -47,3 +47,4 @@ from .mesh import * from .filter import * from .tally import * from .settings import settings +from .math import * diff --git a/openmc/capi/math.py b/openmc/capi/math.py new file mode 100644 index 000000000..d1d7abdf6 --- /dev/null +++ b/openmc/capi/math.py @@ -0,0 +1,250 @@ +from ctypes import (c_int, c_double, POINTER) + +import numpy as np +from numpy.ctypeslib import ndpointer + +from . import _dll + + +_dll.t_percentile_c.restype = c_double +_dll.t_percentile_c.argtypes = [c_double, c_int] + +_dll.calc_pn_c.restype = None +_dll.calc_pn_c.argtypes = [c_int, c_double, ndpointer(c_double)] + +_dll.evaluate_legendre_c.restype = c_double +_dll.evaluate_legendre_c.argtypes = [c_int, POINTER(c_double), c_double] + +_dll.calc_rn_c.restype = None +_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.rotate_angle_c.restype = None +_dll.rotate_angle_c.argtypes = [ndpointer(c_double), c_double, + POINTER(c_double)] +_dll.maxwell_spectrum_c.restype = c_double +_dll.maxwell_spectrum_c.argtypes = [c_double] + +_dll.watt_spectrum_c.restype = c_double +_dll.watt_spectrum_c.argtypes = [c_double, c_double] + +_dll.broaden_wmp_polynomials_c.restype = None +_dll.broaden_wmp_polynomials_c.argtypes = [c_double, c_double, c_int, + ndpointer(c_double)] + + +def t_percentile(p, df): + """ Calculate the percentile of the Student's t distribution with a + specified probability level and number of degrees of freedom + + Parameters + ---------- + p : float + Probability level + df : int + Degrees of freedom + + Returns + ------- + float + Corresponding t-value + + """ + + return _dll.t_percentile_c(p, df) + + +def calc_pn(n, x): + """ Calculate the n-th order Legendre polynomial at the value of x. + + Parameters + ---------- + n : int + Legendre order + x : float + Independent variable to evaluate the Legendre at + + Returns + ------- + float + Corresponding Legendre polynomial result + + """ + + pnx = np.empty(n + 1, dtype=np.float64) + _dll.calc_pn_c(n, x, pnx) + return pnx + + +def evaluate_legendre(data, x): + """ Finds the value of f(x) given a set of Legendre coefficients + and the value of x. + + Parameters + ---------- + data : iterable of float + Legendre coefficients + x : float + Independent variable to evaluate the Legendre at + + Returns + ------- + float + Corresponding Legendre expansion result + + """ + + data_arr = np.array(data, dtype=np.float64) + return _dll.evaluate_legendre_c(len(data), + data_arr.ctypes.data_as(POINTER(c_double)), + x) + + +def calc_rn(n, uvw): + """ Calculate the n-th order real Spherical Harmonics for a given angle; + all Rn,m values are provided for all n (where -n <= m <= n). + + Parameters + ---------- + n : int + Harmonics order + uvw : iterable of float + Independent variable to evaluate the Legendre at + + Returns + ------- + numpy.ndarray + Corresponding real harmonics value + + """ + + num_nm = (n + 1) * (n + 1) + rn = np.empty(num_nm, dtype=np.float64) + uvw_arr = np.array(uvw, dtype=np.float64) + _dll.calc_rn_c(n, uvw_arr, rn) + return rn + + +def calc_zn(n, rho, phi): + """ Calculate 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 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 + phi : float + Theta (radians) location in the unit disk + + Returns + ------- + numpy.ndarray + Corresponding resulting list of coefficients + + """ + + num_bins = ((n + 1) * (n + 2)) // 2 + zn = np.zeros(num_bins, dtype=np.float64) + _dll.calc_zn_c(n, rho, phi, zn) + return zn + + +def rotate_angle(uvw0, mu, phi=None): + """ Rotates direction cosines through a polar angle whose cosine is + mu and through an azimuthal angle sampled uniformly. + + Parameters + ---------- + uvw0 : iterable of float + Original direction cosine + mu : float + Polar angle cosine to rotate + phi : float, optional + Azimuthal angle; if None, one will be sampled uniformly + + Returns + ------- + numpy.ndarray + Rotated direction cosine + + """ + + uvw0_arr = np.array(uvw0, dtype=np.float64) + + if phi is None: + _dll.rotate_angle_c(uvw0_arr, mu, None) + else: + _dll.rotate_angle_c(uvw0_arr, mu, c_double(phi)) + uvw = uvw0_arr + + return uvw + + +def maxwell_spectrum(T): + """ Samples an energy from the Maxwell fission distribution based + on a direct sampling scheme. + + Parameters + ---------- + T : float + Spectrum parameter + + Returns + ------- + float + Sampled outgoing energy + + """ + + return _dll.maxwell_spectrum_c(T) + + +def watt_spectrum(a, b): + """ Samples an energy from the Watt energy-dependent fission spectrum. + + Parameters + ---------- + a : float + Spectrum parameter a + b : float + Spectrum parameter b + + Returns + ------- + float + Sampled outgoing energy + + """ + + return _dll.watt_spectrum_c(a, b) + + +def broaden_wmp_polynomials(E, dopp, n): + """ Doppler broadens the windowed multipole curvefit. The curvefit is a + polynomial of the form a/E + b/sqrt(E) + c + d sqrt(E) ... + + Parameters + ---------- + E : float + Energy to evaluate at + dopp : float + sqrt(atomic weight ratio / kT), with kT given in eV + n : int + Number of components to the polynomial + + Returns + ------- + numpy.ndarray + Resultant leading coefficients + + """ + + factors = np.zeros(n, dtype=np.float64) + _dll.broaden_wmp_polynomials_c(E, dopp, n, factors) + return factors diff --git a/src/api.F90 b/src/api.F90 index bc0757be5..da50007d6 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -10,6 +10,7 @@ module openmc_api use geometry_header use hdf5_interface use material_header + use math use mesh_header use message_passing use nuclide_header diff --git a/src/distribution_multivariate.F90 b/src/distribution_multivariate.F90 index 650ab26fa..19db1c297 100644 --- a/src/distribution_multivariate.F90 +++ b/src/distribution_multivariate.F90 @@ -119,7 +119,7 @@ contains else ! Sample azimuthal angle phi = this % phi % sample() - uvw(:) = rotate_angle(this % reference_uvw, mu, phi) + uvw = rotate_angle(this % reference_uvw, mu, phi) end if end function polar_azimuthal_sample diff --git a/src/faddeeva/Faddeeva.h b/src/faddeeva/Faddeeva.h index 429386190..9e26bc1ed 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 c63313164..33f1ab4d3 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -6,6 +6,18 @@ module math use random_lcg, only: prn implicit none + private + public :: t_percentile + public :: calc_pn + public :: calc_rn + public :: calc_zn + public :: evaluate_legendre + public :: rotate_angle + public :: maxwell_spectrum + public :: watt_spectrum + public :: faddeeva + public :: w_derivative + public :: broaden_wmp_polynomials !=============================================================================== ! FADDEEVA_W evaluates the scaled complementary error function. This @@ -13,6 +25,86 @@ module math !=============================================================================== interface + + pure function t_percentile(p, df) bind(C, name='t_percentile_c') & + result(t) + use ISO_C_BINDING + implicit none + real(C_DOUBLE), value, intent(in) :: p + integer(C_INT), value, intent(in) :: df + real(C_DOUBLE) :: t + end function t_percentile + + pure subroutine calc_pn(n, x, pnx) bind(C, name='calc_pn_c') + use ISO_C_BINDING + implicit none + integer(C_INT), value, intent(in) :: n + real(C_DOUBLE), value, intent(in) :: x + real(C_DOUBLE), intent(out) :: pnx(n + 1) + end subroutine calc_pn + + pure function evaluate_legendre_c_intfc(n, data, x) & + bind(C, name='evaluate_legendre_c') result(val) + use ISO_C_BINDING + implicit none + integer(C_INT), value, intent(in) :: n + real(C_DOUBLE), intent(in) :: data(n) + real(C_DOUBLE), value, intent(in) :: x + real(C_DOUBLE) :: val + end function evaluate_legendre_c_intfc + + pure subroutine calc_rn(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 + + pure subroutine calc_zn(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 + + subroutine rotate_angle_c_intfc(uvw, mu, phi) bind(C, name='rotate_angle_c') + use ISO_C_BINDING + implicit none + real(C_DOUBLE), intent(inout) :: uvw(3) + real(C_DOUBLE), value, intent(in) :: mu + real(C_DOUBLE), optional, intent(in) :: phi + end subroutine rotate_angle_c_intfc + + function maxwell_spectrum(T) bind(C, name='maxwell_spectrum_c') & + result(E_out) + use ISO_C_BINDING + implicit none + real(C_DOUBLE), value, intent(in) :: T + real(C_DOUBLE) :: E_out + end function maxwell_spectrum + + function watt_spectrum(a, b) bind(C, name='watt_spectrum_c') & + result(E_out) + use ISO_C_BINDING + implicit none + real(C_DOUBLE), value, intent(in) :: a + real(C_DOUBLE), value, intent(in) :: b + real(C_DOUBLE) :: E_out + end function watt_spectrum + + subroutine broaden_wmp_polynomials(E, dopp, n, factors) & + bind(C, name='broaden_wmp_polynomials_c') + use ISO_C_BINDING + implicit none + real(C_DOUBLE), value, intent(in) :: E + real(C_DOUBLE), value, intent(in) :: dopp + integer(C_INT), value, intent(in) :: n + real(C_DOUBLE), intent(inout) :: factors(n) + end subroutine broaden_wmp_polynomials + function faddeeva_w(z, relerr) bind(C, name='Faddeeva_w') result(w) use ISO_C_BINDING implicit none @@ -24,696 +116,17 @@ module math contains -!=============================================================================== -! NORMAL_PERCENTILE calculates the percentile of the standard normal -! distribution with a specified probability level -!=============================================================================== - - elemental function normal_percentile(p) result(z) - - real(8), intent(in) :: p ! probability level - real(8) :: z ! corresponding z-value - - real(8) :: q - real(8) :: r - real(8), parameter :: p_low = 0.02425_8 - real(8), parameter :: a(6) = (/ & - -3.969683028665376e1_8, 2.209460984245205e2_8, -2.759285104469687e2_8, & - 1.383577518672690e2_8, -3.066479806614716e1_8, 2.506628277459239e0_8 /) - real(8), parameter :: b(5) = (/ & - -5.447609879822406e1_8, 1.615858368580409e2_8, -1.556989798598866e2_8, & - 6.680131188771972e1_8, -1.328068155288572e1_8 /) - real(8), parameter :: c(6) = (/ & - -7.784894002430293e-3_8, -3.223964580411365e-1_8, -2.400758277161838_8, & - -2.549732539343734_8, 4.374664141464968_8, 2.938163982698783_8 /) - real(8), parameter :: d(4) = (/ & - 7.784695709041462e-3_8, 3.224671290700398e-1_8, & - 2.445134137142996_8, 3.754408661907416_8 /) - - ! The rational approximation used here is from an unpublished work at - ! http://home.online.no/~pjacklam/notes/invnorm/ - - if (p < p_low) then - ! Rational approximation for lower region. - - q = sqrt(-TWO*log(p)) - z = (((((c(1)*q + c(2))*q + c(3))*q + c(4))*q + c(5))*q + c(6)) / & - ((((d(1)*q + d(2))*q + d(3))*q + d(4))*q + ONE) - - elseif (p <= ONE - p_low) then - ! Rational approximation for central region - - q = p - HALF - r = q*q - z = (((((a(1)*r + a(2))*r + a(3))*r + a(4))*r + a(5))*r + a(6))*q / & - (((((b(1)*r + b(2))*r + b(3))*r + b(4))*r + b(5))*r + ONE) - - else - ! Rational approximation for upper region - - q = sqrt(-TWO*log(ONE - p)) - z = -(((((c(1)*q + c(2))*q + c(3))*q + c(4))*q + c(5))*q + c(6)) / & - ((((d(1)*q + d(2))*q + d(3))*q + d(4))*q + ONE) - endif - - ! Refinement based on Newton's method -#ifndef NO_F2008 - z = z - (HALF * erfc(-z/sqrt(TWO)) - p) * sqrt(TWO*PI) * exp(HALF*z*z) -#endif - - end function normal_percentile - -!=============================================================================== -! T_PERCENTILE calculates the percentile of the Student's t distribution with a -! specified probability level and number of degrees of freedom -!=============================================================================== - - elemental function t_percentile(p, df) result(t) - - real(8), intent(in) :: p ! probability level - integer, intent(in) :: df ! degrees of freedom - real(8) :: t ! corresponding t-value - - real(8) :: n ! degrees of freedom as a real(8) - real(8) :: k ! n - 2 - real(8) :: z ! percentile of normal distribution - real(8) :: z2 ! z * z - - if (df == 1) then - ! For one degree of freedom, the t-distribution becomes a Cauchy - ! distribution whose cdf we can invert directly - - t = tan(PI*(p - HALF)) - - elseif (df == 2) then - ! For two degrees of freedom, the cdf is given by 1/2 + x/(2*sqrt(x^2 + - ! 2)). This can be directly inverted to yield the solution below - - t = TWO*sqrt(TWO)*(p - HALF)/sqrt(ONE - FOUR*(p - HALF)**2) - - else - - ! This approximation is from E. Olusegun George and Meenakshi Sivaram, "A - ! modification of the Fisher-Cornish approximation for the student t - ! percentiles," Communication in Statistics - Simulation and Computation, - ! 16 (4), pp. 1123-1132 (1987). - - n = real(df,8) - k = ONE/(n - TWO) - z = normal_percentile(p) - z2 = z * z - t = sqrt(n*k) * (z + (z2 - THREE)*z*k/FOUR + ((5._8*z2 - 56._8)*z2 + & - 75._8)*z*k*k/96._8 + (((z2 - 27._8)*THREE*z2 + 417._8)*z2 - 315._8) & - *z*k*k*k/384._8) - - end if - - end function t_percentile - -!=============================================================================== -! CALC_PN calculates the n-th order Legendre polynomial at the value of x. -! Since this function is called repeatedly during the neutron transport process, -! neither n or x is checked to see if they are in the applicable range. -! This is left to the client developer to use where applicable. x is to be in -! the domain of [-1,1], and 0<=n<=5. If x is outside of the range, the return -! value will be outside the expected range; if n is outside the stated range, -! the return value will be 1.0. -!=============================================================================== - - elemental function calc_pn(n,x) result(pnx) - - integer, intent(in) :: n ! Legendre order requested - real(8), intent(in) :: x ! Independent variable the Legendre is to be - ! evaluated at; x must be in the domain [-1,1] - real(8) :: pnx ! The Legendre poly of order n evaluated at x - - select case(n) - case(1) - pnx = x - case(2) - pnx = 1.5_8 * x * x - HALF - case(3) - pnx = 2.5_8 * x * x * x - 1.5_8 * x - case(4) - pnx = 4.375_8 * (x ** 4) - 3.75_8 * x * x + 0.375_8 - case(5) - pnx = 7.875_8 * (x ** 5) - 8.75_8 * x * x * x + 1.875 * x - case(6) - pnx = 14.4375_8 * (x ** 6) - 19.6875_8 * (x ** 4) + & - 6.5625_8 * x * x - 0.3125_8 - case(7) - pnx = 26.8125_8 * (x ** 7) - 43.3125_8 * (x ** 5) + & - 19.6875_8 * x * x * x - 2.1875_8 * x - case(8) - pnx = 50.2734375_8 * (x ** 8) - 93.84375_8 * (x ** 6) + & - 54.140625 * (x ** 4) - 9.84375_8 * x * x + 0.2734375_8 - case(9) - pnx = 94.9609375_8 * (x ** 9) - 201.09375_8 * (x ** 7) + & - 140.765625_8 * (x ** 5) - 36.09375_8 * x * x * x + 2.4609375_8 * x - case(10) - pnx = 180.42578125_8 * (x ** 10) - 427.32421875_8 * (x ** 8) + & - 351.9140625_8 * (x ** 6) - 117.3046875_8 * (x ** 4) + & - 13.53515625_8 * x * x - 0.24609375_8 - case default - pnx = ONE ! correct for case(0), incorrect for the rest - end select - - end function calc_pn - -!=============================================================================== -! CALC_RN calculates the n-th order real spherical harmonics for a given angle -! (in terms of (u,v,w)). All Rn,m values are provided (where -n<=m<=n) -!=============================================================================== - - pure function calc_rn(n,uvw) result(rn) - - integer, intent(in) :: n ! Order requested - real(8), intent(in) :: uvw(3) ! Direction of travel, assumed to be on unit sphere - real(8) :: rn(2*n + 1) ! The resultant R_n(uvw) - - real(8) :: phi, w ! Azimuthal and Cosine of Polar angles (from uvw) - real(8) :: w2m1 ! (w^2 - 1), frequently used in these - - w = uvw(3) ! z = cos(polar) - if (uvw(1) == ZERO) then - phi = ZERO - else - phi = atan2(uvw(2), uvw(1)) - end if - - w2m1 = (ONE - w**2) - select case(n) - case (0) - ! l = 0, m = 0 - rn(1) = ONE - case (1) - ! l = 1, m = -1 - rn(1) = -(ONE*sqrt(w2m1) * sin(phi)) - ! l = 1, m = 0 - rn(2) = ONE * w - ! l = 1, m = 1 - rn(3) = -(ONE*sqrt(w2m1) * cos(phi)) - case (2) - ! l = 2, m = -2 - rn(1) = 0.288675134594813_8 * (-THREE * w**2 + THREE) * sin(TWO*phi) - ! l = 2, m = -1 - rn(2) = -(1.73205080756888_8 * w*sqrt(w2m1) * sin(phi)) - ! l = 2, m = 0 - rn(3) = 1.5_8 * w**2 - HALF - ! l = 2, m = 1 - rn(4) = -(1.73205080756888_8 * w*sqrt(w2m1) * cos(phi)) - ! l = 2, m = 2 - rn(5) = 0.288675134594813_8 * (-THREE * w**2 + THREE) * cos(TWO*phi) - case (3) - ! l = 3, m = -3 - rn(1) = -(0.790569415042095_8 * (w2m1)**(THREE/TWO) * sin(THREE * phi)) - ! l = 3, m = -2 - rn(2) = 1.93649167310371_8 * w*(w2m1) * sin(TWO*phi) - ! l = 3, m = -1 - rn(3) = -(0.408248290463863_8*sqrt(w2m1)*((15.0_8/TWO)*w**2 - THREE/TWO) * & - sin(phi)) - ! l = 3, m = 0 - rn(4) = 2.5_8 * w**3 - 1.5_8 * w - ! l = 3, m = 1 - rn(5) = -(0.408248290463863_8*sqrt(w2m1)*((15.0_8/TWO)*w**2 - THREE/TWO) * & - cos(phi)) - ! l = 3, m = 2 - rn(6) = 1.93649167310371_8 * w*(w2m1) * cos(TWO*phi) - ! l = 3, m = 3 - rn(7) = -(0.790569415042095_8 * (w2m1)**(THREE/TWO) * cos(THREE* phi)) - case (4) - ! l = 4, m = -4 - rn(1) = 0.739509972887452_8 * (w2m1)**2 * sin(4.0_8*phi) - ! l = 4, m = -3 - rn(2) = -(2.09165006633519_8 * w*(w2m1)**(THREE/TWO) * sin(THREE* phi)) - ! l = 4, m = -2 - rn(3) = 0.074535599249993_8 * (w2m1)*((105.0_8/TWO)*w**2 - 15.0_8/TWO) * & - sin(TWO*phi) - ! l = 4, m = -1 - rn(4) = -(0.316227766016838_8*sqrt(w2m1)*((35.0_8/TWO)*w**3 - 15.0_8/TWO*w)& - * sin(phi)) - ! l = 4, m = 0 - rn(5) = 4.375_8 * w**4 - 3.75_8 * w**2 + 0.375_8 - ! l = 4, m = 1 - rn(6) = -(0.316227766016838_8*sqrt(w2m1)*((35.0_8/TWO)*w**3 - 15.0_8/TWO*w)& - * cos(phi)) - ! l = 4, m = 2 - rn(7) = 0.074535599249993_8 * (w2m1)*((105.0_8/TWO)*w**2 - 15.0_8/TWO) * & - cos(TWO*phi) - ! l = 4, m = 3 - rn(8) = -(2.09165006633519_8 * w*(w2m1)**(THREE/TWO) * cos(THREE* phi)) - ! l = 4, m = 4 - rn(9) = 0.739509972887452_8 * (w2m1)**2 * cos(4.0_8*phi) - case (5) - ! l = 5, m = -5 - rn(1) = -(0.701560760020114_8 * (w2m1)**(5.0_8/TWO) * sin(5.0_8*phi)) - ! l = 5, m = -4 - rn(2) = 2.21852991866236_8 * w*(w2m1)**2 * sin(4.0_8*phi) - ! l = 5, m = -3 - rn(3) = -(0.00996023841111995_8 * (w2m1)**(THREE/TWO)* & - ((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * sin(THREE*phi)) - ! l = 5, m = -2 - rn(4) = 0.0487950036474267_8 * (w2m1) & - * ((315.0_8/TWO)*w**3 - 105.0_8/TWO*w) * sin(TWO*phi) - ! l = 5, m = -1 - rn(5) = -(0.258198889747161_8*sqrt(w2m1)* & - ((315.0_8/8.0_8)*w**4 - 105.0_8/4.0_8 * w**2 + 15.0_8/8.0_8) & - * sin(phi)) - ! l = 5, m = 0 - rn(6) = 7.875_8 * w**5 - 8.75_8 * w**3 + 1.875_8 * w - ! l = 5, m = 1 - rn(7) = -(0.258198889747161_8*sqrt(w2m1)* & - ((315.0_8/8.0_8)*w**4 - 105.0_8/4.0_8 * w**2 + 15.0_8/8.0_8) & - * cos(phi)) - ! l = 5, m = 2 - rn(8) = 0.0487950036474267_8 * (w2m1)* & - ((315.0_8/TWO)*w**3 - 105.0_8/TWO*w) * cos(TWO*phi) - ! l = 5, m = 3 - rn(9) = -(0.00996023841111995_8 * (w2m1)**(THREE/TWO)* & - ((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * cos(THREE*phi)) - ! l = 5, m = 4 - rn(10) = 2.21852991866236_8 * w*(w2m1)**2 * cos(4.0_8*phi) - ! l = 5, m = 5 - rn(11) = -(0.701560760020114_8 * (w2m1)**(5.0_8/TWO) * cos(5.0_8* phi)) - case (6) - ! l = 6, m = -6 - rn(1) = 0.671693289381396_8 * (w2m1)**3 * sin(6.0_8*phi) - ! l = 6, m = -5 - rn(2) = -(2.32681380862329_8 * w*(w2m1)**(5.0_8/TWO) * sin(5.0_8*phi)) - ! l = 6, m = -4 - rn(3) = 0.00104990131391452_8 * (w2m1)**2 * & - ((10395.0_8/TWO)*w**2 - 945.0_8/TWO) * sin(4.0_8*phi) - ! l = 6, m = -3 - rn(4) = -(0.00575054632785295_8 * (w2m1)**(THREE/TWO) * & - ((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * sin(THREE*phi)) - ! l = 6, m = -2 - rn(5) = 0.0345032779671177_8 * (w2m1) * & - ((3465.0_8/8.0_8)*w**4 - 945.0_8/4.0_8 * w**2 + 105.0_8/8.0_8) & - * sin(TWO*phi) - ! l = 6, m = -1 - rn(6) = -(0.218217890235992_8*sqrt(w2m1) * & - ((693.0_8/8.0_8)*w**5- 315.0_8/4.0_8 * w**3 + (105.0_8/8.0_8)*w) & - * sin(phi)) - ! l = 6, m = 0 - rn(7) = 14.4375_8 * w**6 - 19.6875_8 * w**4 + 6.5625_8 * w**2 - 0.3125_8 - ! l = 6, m = 1 - rn(8) = -(0.218217890235992_8*sqrt(w2m1) * & - ((693.0_8/8.0_8)*w**5- 315.0_8/4.0_8 * w**3 + (105.0_8/8.0_8)*w) & - * cos(phi)) - ! l = 6, m = 2 - rn(9) = 0.0345032779671177_8 * (w2m1) * & - ((3465.0_8/8.0_8)*w**4 -945.0_8/4.0_8 * w**2 + 105.0_8/8.0_8) & - * cos(TWO*phi) - ! l = 6, m = 3 - rn(10) = -(0.00575054632785295_8 * (w2m1)**(THREE/TWO) * & - ((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * cos(THREE*phi)) - ! l = 6, m = 4 - rn(11) = 0.00104990131391452_8 * (w2m1)**2 * & - ((10395.0_8/TWO)*w**2 - 945.0_8/TWO) * cos(4.0_8*phi) - ! l = 6, m = 5 - rn(12) = -(2.32681380862329_8 * w*(w2m1)**(5.0_8/TWO) * cos(5.0_8*phi)) - ! l = 6, m = 6 - rn(13) = 0.671693289381396_8 * (w2m1)**3 * cos(6.0_8*phi) - case (7) - ! l = 7, m = -7 - rn(1) = -(0.647259849287749_8 * (w2m1)**(7.0_8/TWO) * sin(7.0_8*phi)) - ! l = 7, m = -6 - rn(2) = 2.42182459624969_8 * w*(w2m1)**3 * sin(6.0_8*phi) - ! l = 7, m = -5 - rn(3) = -(9.13821798555235d-5*(w2m1)**(5.0_8/TWO)* & - ((135135.0_8/TWO)*w**2 - 10395.0_8/TWO) * sin(5.0_8*phi)) - ! l = 7, m = -4 - rn(4) = 0.000548293079133141_8 * (w2m1)**2* & - ((45045.0_8/TWO)*w**3 - 10395.0_8/TWO*w) * sin(4.0_8*phi) - ! l = 7, m = -3 - rn(5) = -(0.00363696483726654_8 * (w2m1)**(THREE/TWO)* & - ((45045.0_8/8.0_8)*w**4 - 10395.0_8/4.0_8 * w**2 + 945.0_8/8.0_8)* & - sin(THREE*phi)) - ! l = 7, m = -2 - rn(6) = 0.025717224993682_8 * (w2m1)* & - ((9009.0_8/8.0_8)*w**5 -3465.0_8/4.0_8 * w**3 + (945.0_8/8.0_8)*w)* & - sin(TWO*phi) - ! l = 7, m = -1 - rn(7) = -(0.188982236504614_8*sqrt(w2m1)* & - ((3003.0_8/16.0_8)*w**6 - 3465.0_8/16.0_8 * w**4 + & - (945.0_8/16.0_8)*w**2 - 35.0_8/16.0_8) * sin(phi)) - ! l = 7, m = 0 - rn(8) = 26.8125_8 * w**7 - 43.3125_8 * w**5 + 19.6875_8 * w**3 -2.1875_8 & - * w - ! l = 7, m = 1 - rn(9) = -(0.188982236504614_8*sqrt(w2m1)* & - ((3003.0_8/16.0_8)*w**6 - 3465.0_8/16.0_8 * w**4 + & - (945.0_8/16.0_8)*w**2 - 35.0_8/16.0_8) * cos(phi)) - ! l = 7, m = 2 - rn(10) = 0.025717224993682_8 * (w2m1)* & - ((9009.0_8/8.0_8)*w**5 -3465.0_8/4.0_8 * w**3 + (945.0_8/8.0_8)*w)* & - cos(TWO*phi) - ! l = 7, m = 3 - rn(11) = -(0.00363696483726654_8 * (w2m1)**(THREE/TWO)* & - ((45045.0_8/8.0_8)*w**4 - 10395.0_8/4.0_8 * w**2 + 945.0_8/8.0_8)* & - cos(THREE*phi)) - ! l = 7, m = 4 - rn(12) = 0.000548293079133141_8 * (w2m1)**2 * & - ((45045.0_8/TWO)*w**3 - 10395.0_8/TWO*w) * cos(4.0_8*phi) - ! l = 7, m = 5 - rn(13) = -(9.13821798555235d-5*(w2m1)**(5.0_8/TWO)* & - ((135135.0_8/TWO)*w**2 - 10395.0_8/TWO) * cos(5.0_8*phi)) - ! l = 7, m = 6 - rn(14) = 2.42182459624969_8 * w*(w2m1)**3 * cos(6.0_8*phi) - ! l = 7, m = 7 - rn(15) = -(0.647259849287749_8 * (w2m1)**(7.0_8/TWO) * cos(7.0_8*phi)) - case (8) - ! l = 8, m = -8 - rn(1) = 0.626706654240044_8 * (w2m1)**4 * sin(8.0_8*phi) - ! l = 8, m = -7 - rn(2) = -(2.50682661696018_8 * w*(w2m1)**(7.0_8/TWO) * sin(7.0_8*phi)) - ! l = 8, m = -6 - rn(3) = 6.77369783729086d-6*(w2m1)**3* & - ((2027025.0_8/TWO)*w**2 - 135135.0_8/TWO) * sin(6.0_8*phi) - ! l = 8, m = -5 - rn(4) = -(4.38985792528482d-5*(w2m1)**(5.0_8/TWO)* & - ((675675.0_8/TWO)*w**3 - 135135.0_8/TWO*w) * sin(5.0_8*phi)) - ! l = 8, m = -4 - rn(5) = 0.000316557156832328_8 * (w2m1)**2* & - ((675675.0_8/8.0_8)*w**4 - 135135.0_8/4.0_8 * w**2 & - + 10395.0_8/8.0_8) * sin(4.0_8*phi) - ! l = 8, m = -3 - rn(6) = -(0.00245204119306875_8 * (w2m1)**(THREE/TWO)* & - ((135135.0_8/8.0_8)*w**5 - 45045.0_8/4.0_8 * w**3 & - + (10395.0_8/8.0_8)*w) * sin(THREE*phi)) - ! l = 8, m = -2 - rn(7) = 0.0199204768222399_8 * (w2m1)* & - ((45045.0_8/16.0_8)*w**6- 45045.0_8/16.0_8 * w**4 + & - (10395.0_8/16.0_8)*w**2 - 315.0_8/16.0_8) * sin(TWO*phi) - ! l = 8, m = -1 - rn(8) = -(0.166666666666667_8*sqrt(w2m1)* & - ((6435.0_8/16.0_8)*w**7 - 9009.0_8/16.0_8 * w**5 + & - (3465.0_8/16.0_8)*w**3 - 315.0_8/16.0_8 * w) * sin(phi)) - ! l = 8, m = 0 - rn(9) = 50.2734375_8 * w**8 - 93.84375_8 * w**6 + 54.140625_8 * w**4 -& - 9.84375_8 * w**2 + 0.2734375_8 - ! l = 8, m = 1 - rn(10) = -(0.166666666666667_8*sqrt(w2m1)* & - ((6435.0_8/16.0_8)*w**7 - 9009.0_8/16.0_8 * w**5 + & - (3465.0_8/16.0_8)*w**3 - 315.0_8/16.0_8 * w) * cos(phi)) - ! l = 8, m = 2 - rn(11) = 0.0199204768222399_8 * (w2m1)*((45045.0_8/16.0_8)*w**6- & - 45045.0_8/16.0_8 * w**4 + (10395.0_8/16.0_8)*w**2 - & - 315.0_8/16.0_8) * cos(TWO*phi) - ! l = 8, m = 3 - rn(12) = -(0.00245204119306875_8 * (w2m1)**(THREE/TWO)* & - ((135135.0_8/8.0_8)*w**5 - 45045.0_8/4.0_8 * w**3 + & - (10395.0_8/8.0_8)*w) * cos(THREE*phi)) - ! l = 8, m = 4 - rn(13) = 0.000316557156832328_8 * (w2m1)**2*((675675.0_8/8.0_8)*w**4 - & - 135135.0_8/4.0_8 * w**2 + 10395.0_8/8.0_8) * cos(4.0_8*phi) - ! l = 8, m = 5 - rn(14) = -(4.38985792528482d-5*(w2m1)**(5.0_8/TWO)*((675675.0_8/TWO)*w**3 -& - 135135.0_8/TWO*w) * cos(5.0_8*phi)) - ! l = 8, m = 6 - rn(15) = 6.77369783729086d-6*(w2m1)**3*((2027025.0_8/TWO)*w**2 - & - 135135.0_8/TWO) * cos(6.0_8*phi) - ! l = 8, m = 7 - rn(16) = -(2.50682661696018_8 * w*(w2m1)**(7.0_8/TWO) * cos(7.0_8*phi)) - ! l = 8, m = 8 - rn(17) = 0.626706654240044_8 * (w2m1)**4 * cos(8.0_8*phi) - case (9) - ! l = 9, m = -9 - rn(1) = -(0.609049392175524_8 * (w2m1)**(9.0_8/TWO) * sin(9.0_8*phi)) - ! l = 9, m = -8 - rn(2) = 2.58397773170915_8 * w*(w2m1)**4 * sin(8.0_8*phi) - ! l = 9, m = -7 - rn(3) = -(4.37240315267812d-7*(w2m1)**(7.0_8/TWO)* & - ((34459425.0_8/TWO)*w**2 - 2027025.0_8/TWO) * sin(7.0_8*phi)) - ! l = 9, m = -6 - rn(4) = 3.02928976464514d-6*(w2m1)**3* & - ((11486475.0_8/TWO)*w**3 - 2027025.0_8/TWO*w) * sin(6.0_8*phi) - ! l = 9, m = -5 - rn(5) = -(2.34647776186144d-5*(w2m1)**(5.0_8/TWO)* & - ((11486475.0_8/8.0_8)*w**4 - 2027025.0_8/4.0_8 * w**2 + & - 135135.0_8/8.0_8) * sin(5.0_8*phi)) - ! l = 9, m = -4 - rn(6) = 0.000196320414650061_8 * (w2m1)**2*((2297295.0_8/8.0_8)*w**5 - & - 675675.0_8/4.0_8 * w**3 + (135135.0_8/8.0_8)*w) * sin(4.0_8*phi) - ! l = 9, m = -3 - rn(7) = -(0.00173385495536766_8 * (w2m1)**(THREE/TWO)* & - ((765765.0_8/16.0_8)*w**6 - 675675.0_8/16.0_8 * w**4 + & - (135135.0_8/16.0_8)*w**2 - 3465.0_8/16.0_8) * sin(THREE*phi)) - ! l = 9, m = -2 - rn(8) = 0.0158910431540932_8 * (w2m1)*((109395.0_8/16.0_8)*w**7- & - 135135.0_8/16.0_8 * w**5 + (45045.0_8/16.0_8)*w**3 & - - 3465.0_8/16.0_8 * w) * sin(TWO*phi) - ! l = 9, m = -1 - rn(9) = -(0.149071198499986_8*sqrt(w2m1)*((109395.0_8/128.0_8)*w**8 - & - 45045.0_8/32.0_8 * w**6 + (45045.0_8/64.0_8)*w**4 - 3465.0_8/32.0_8 & - * w**2 + 315.0_8/128.0_8) * sin(phi)) - ! l = 9, m = 0 - rn(10) = 94.9609375_8 * w**9 - 201.09375_8 * w**7 + 140.765625_8 * w**5- & - 36.09375_8 * w**3 + 2.4609375_8 * w - ! l = 9, m = 1 - rn(11) = -(0.149071198499986_8*sqrt(w2m1)*((109395.0_8/128.0_8)*w**8 - & - 45045.0_8/32.0_8 * w**6 + (45045.0_8/64.0_8)*w**4 -3465.0_8/32.0_8 & - * w**2 + 315.0_8/128.0_8) * cos(phi)) - ! l = 9, m = 2 - rn(12) = 0.0158910431540932_8 * (w2m1)*((109395.0_8/16.0_8)*w**7 - & - 135135.0_8/16.0_8 * w**5 + (45045.0_8/16.0_8)*w**3 & - - 3465.0_8/ 16.0_8 * w) * cos(TWO*phi) - ! l = 9, m = 3 - rn(13) = -(0.00173385495536766_8 * (w2m1)**(THREE/TWO)*((765765.0_8/16.0_8)& - *w**6 - 675675.0_8/16.0_8 * w**4 + (135135.0_8/16.0_8)*w**2 & - - 3465.0_8/16.0_8)* cos(THREE*phi)) - ! l = 9, m = 4 - rn(14) = 0.000196320414650061_8 * (w2m1)**2*((2297295.0_8/8.0_8)*w**5 - & - 675675.0_8/4.0_8 * w**3 + (135135.0_8/8.0_8)*w) * cos(4.0_8*phi) - ! l = 9, m = 5 - rn(15) = -(2.34647776186144d-5*(w2m1)**(5.0_8/TWO)*((11486475.0_8/8.0_8)* & - w**4 - 2027025.0_8/4.0_8 * w**2 + 135135.0_8/8.0_8) * cos(5.0_8*phi)) - ! l = 9, m = 6 - rn(16) = 3.02928976464514d-6*(w2m1)**3*((11486475.0_8/TWO)*w**3 - & - 2027025.0_8/TWO*w) * cos(6.0_8*phi) - ! l = 9, m = 7 - rn(17) = -(4.37240315267812d-7*(w2m1)**(7.0_8/TWO)* & - ((34459425.0_8/TWO)*w**2 - 2027025.0_8/TWO) * cos(7.0_8*phi)) - ! l = 9, m = 8 - rn(18) = 2.58397773170915_8 * w*(w2m1)**4 * cos(8.0_8*phi) - ! l = 9, m = 9 - rn(19) = -(0.609049392175524_8 * (w2m1)**(9.0_8/TWO) * cos(9.0_8*phi)) - case (10) - ! l = 10, m = -10 - rn(1) = 0.593627917136573_8 * (w2m1)**5 * sin(10.0_8*phi) - ! l = 10, m = -9 - rn(2) = -(2.65478475211798_8 * w*(w2m1)**(9.0_8/TWO) * sin(9.0_8*phi)) - ! l = 10, m = -8 - rn(3) = 2.49953651452314d-8*(w2m1)**4*((654729075.0_8/TWO)*w**2 - & - 34459425.0_8/TWO) * sin(8.0_8*phi) - ! l = 10, m = -7 - rn(4) = -(1.83677671621093d-7*(w2m1)**(7.0_8/TWO)* & - ((218243025.0_8/TWO)*w**3 - 34459425.0_8/TWO*w) * sin(7.0_8*phi)) - ! l = 10, m = -6 - rn(5) = 1.51464488232257d-6*(w2m1)**3*((218243025.0_8/8.0_8)*w**4 - & - 34459425.0_8/4.0_8 * w**2 + 2027025.0_8/8.0_8) * sin(6.0_8*phi) - ! l = 10, m = -5 - rn(6) = -(1.35473956745817d-5*(w2m1)**(5.0_8/TWO)* & - ((43648605.0_8/8.0_8)*w**5 - 11486475.0_8/4.0_8 * w**3 + & - (2027025.0_8/8.0_8)*w) * sin(5.0_8*phi)) - ! l = 10, m = -4 - rn(7) = 0.000128521880085575_8 * (w2m1)**2*((14549535.0_8/16.0_8)*w**6 - & - 11486475.0_8/16.0_8 * w**4 + (2027025.0_8/16.0_8)*w**2 - & - 45045.0_8/16.0_8) * sin(4.0_8*phi) - ! l = 10, m = -3 - rn(8) = -(0.00127230170115096_8 * (w2m1)**(THREE/TWO)* & - ((2078505.0_8/16.0_8)*w**7 - 2297295.0_8/16.0_8 * w**5 + & - (675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * sin(THREE*phi)) - ! l = 10, m = -2 - rn(9) = 0.012974982402692_8 * (w2m1)*((2078505.0_8/128.0_8)*w**8 - & - 765765.0_8/32.0_8 * w**6 + (675675.0_8/64.0_8)*w**4 - & - 45045.0_8/32.0_8 * w**2 + 3465.0_8/128.0_8) * sin(TWO*phi) - ! l = 10, m = -1 - rn(10) = -(0.134839972492648_8*sqrt(w2m1)*((230945.0_8/128.0_8)*w**9 - & - 109395.0_8/32.0_8 * w**7 + (135135.0_8/64.0_8)*w**5 - & - 15015.0_8/32.0_8 * w**3 + (3465.0_8/128.0_8)*w) * sin(phi)) - ! l = 10, m = 0 - rn(11) = 180.42578125_8 * w**10 - 427.32421875_8 * w**8 +351.9140625_8 & - * w**6 - 117.3046875_8 * w**4 + 13.53515625_8 * w**2 -0.24609375_8 - ! l = 10, m = 1 - rn(12) = -(0.134839972492648_8*sqrt(w2m1)*((230945.0_8/128.0_8)*w**9 - & - 109395.0_8/32.0_8 * w**7 + (135135.0_8/64.0_8)*w**5 -15015.0_8/ & - 32.0_8 * w**3 + (3465.0_8/128.0_8)*w) * cos(phi)) - ! l = 10, m = 2 - rn(13) = 0.012974982402692_8 * (w2m1)*((2078505.0_8/128.0_8)*w**8 - & - 765765.0_8/32.0_8 * w**6 + (675675.0_8/64.0_8)*w**4 -& - 45045.0_8/32.0_8 * w**2 + 3465.0_8/128.0_8) * cos(TWO*phi) - ! l = 10, m = 3 - rn(14) = -(0.00127230170115096_8 * (w2m1)**(THREE/TWO)* & - ((2078505.0_8/16.0_8)*w**7 - 2297295.0_8/16.0_8 * w**5 + & - (675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * cos(THREE*phi)) - ! l = 10, m = 4 - rn(15) = 0.000128521880085575_8 * (w2m1)**2*((14549535.0_8/16.0_8)*w**6 -& - 11486475.0_8/16.0_8 * w**4 + (2027025.0_8/16.0_8)*w**2 - & - 45045.0_8/16.0_8) * cos(4.0_8*phi) - ! l = 10, m = 5 - rn(16) = -(1.35473956745817d-5*(w2m1)**(5.0_8/TWO)* & - ((43648605.0_8/8.0_8)*w**5 - 11486475.0_8/4.0_8 * w**3 + & - (2027025.0_8/8.0_8)*w) * cos(5.0_8*phi)) - ! l = 10, m = 6 - rn(17) = 1.51464488232257d-6*(w2m1)**3*((218243025.0_8/8.0_8)*w**4 - & - 34459425.0_8/4.0_8 * w**2 + 2027025.0_8/8.0_8) * cos(6.0_8*phi) - ! l = 10, m = 7 - rn(18) = -(1.83677671621093d-7*(w2m1)**(7.0_8/TWO)* & - ((218243025.0_8/TWO)*w**3 - 34459425.0_8/TWO*w) * cos(7.0_8*phi)) - ! l = 10, m = 8 - rn(19) = 2.49953651452314d-8*(w2m1)**4* & - ((654729075.0_8/TWO)*w**2 - 34459425.0_8/TWO) * cos(8.0_8*phi) - ! l = 10, m = 9 - rn(20) = -(2.65478475211798_8 * w*(w2m1)**(9.0_8/TWO) * cos(9.0_8*phi)) - ! l = 10, m = 10 - rn(21) = 0.593627917136573_8 * (w2m1)**5 * cos(10.0_8*phi) - case default - rn = ONE - end select - - end function calc_rn - -!=============================================================================== -! CALC_ZN calculates the n-th order modified Zernike polynomial moment for a -! given angle (rho, theta) location in the unit disk. The normlization of the -! polynomials is such that the integral of Z_pq*Z_pq over the unit disk is -! exactly pi -!=============================================================================== - - subroutine calc_zn(n, rho, phi, 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. - - integer, intent(in) :: n ! Maximum order - real(8), intent(in) :: rho ! Radial location in the unit disk - real(8), intent(in) :: phi ! Theta (radians) location in the unit disk - real(8), intent(out) :: zn(:) ! The resulting list of coefficients - - real(8) :: sin_phi, cos_phi ! Sine and Cosine of phi - real(8) :: sin_phi_vec(n+1) ! Contains sin(n*phi) - real(8) :: cos_phi_vec(n+1) ! Contains cos(n*phi) - real(8) :: zn_mat(n+1, n+1) ! Matrix form of the coefficients which is - ! easier to work with - real(8) :: k1, k2, k3, k4 ! Variables for R_m_n calculation - integer :: i,p,q ! Loop counters - - ! 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)) - else if (q == 0) then - zn(i) = zn_mat(p+1, q+1) - else - zn(i) = zn_mat(p+1, q+1) * cos_phi_vec(abs(q)+1) - end if - i = i + 1 - end do - end do - end subroutine calc_zn - -!=============================================================================== -! EXPAND_HARMONIC expands a given series of real spherical harmonics -!=============================================================================== - - pure function expand_harmonic(data, order, uvw) result(val) - real(8), intent(in) :: data(:) - integer, intent(in) :: order - real(8), intent(in) :: uvw(3) - real(8) :: val - - integer :: l, lm_lo, lm_hi - - val = data(1) - lm_lo = 2 - lm_hi = 4 - do l = 1, order - 1 - val = val + sqrt(TWO * real(l,8) + ONE) * & - dot_product(calc_rn(l,uvw), data(lm_lo:lm_hi)) - lm_lo = lm_hi + 1 - lm_hi = lm_lo + 2 * (l + 1) - end do - - end function expand_harmonic - !=============================================================================== ! EVALUATE_LEGENDRE Find the value of f(x) given a set of Legendre coefficients ! and the value of x !=============================================================================== - pure function evaluate_legendre(data, x) result(val) - real(8), intent(in) :: data(:) - real(8), intent(in) :: x - real(8) :: val + pure function evaluate_legendre(data, x) result(val) bind(C) + real(C_DOUBLE), intent(in) :: data(:) + real(C_DOUBLE), intent(in) :: x + real(C_DOUBLE) :: val - integer :: l - - val = HALF * data(1) - do l = 1, size(data) - 1 - val = val + (real(l,8) + HALF) * data(l + 1) * calc_pn(l,x) - end do + val = evaluate_legendre_c_intfc(size(data) - 1, data, x) end function evaluate_legendre @@ -724,108 +137,24 @@ contains !=============================================================================== function rotate_angle(uvw0, mu, phi) result(uvw) - real(8), intent(in) :: uvw0(3) ! directional cosine - real(8), intent(in) :: mu ! cosine of angle in lab or CM - real(8), optional :: phi ! azimuthal angle - real(8) :: uvw(3) ! rotated directional cosine + real(C_DOUBLE), intent(in) :: uvw0(3) ! directional cosine + real(C_DOUBLE), intent(in) :: mu ! cosine of angle in lab or CM + real(C_DOUBLE), intent(in), optional :: phi ! azimuthal angle - real(8) :: phi_ ! azimuthal angle - real(8) :: sinphi ! sine of azimuthal angle - real(8) :: cosphi ! cosine of azimuthal angle - real(8) :: a ! sqrt(1 - mu^2) - real(8) :: b ! sqrt(1 - w^2) - real(8) :: u0 ! original cosine in x direction - real(8) :: v0 ! original cosine in y direction - real(8) :: w0 ! original cosine in z direction + real(C_DOUBLE) :: uvw(3) ! rotated directional cosine - ! Copy original directional cosines - u0 = uvw0(1) - v0 = uvw0(2) - w0 = uvw0(3) - - ! Sample azimuthal angle in [0,2pi) if none provided - if (present(phi)) then - phi_ = 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 - end if + uvw = uvw0 + call rotate_angle_c_intfc(uvw, mu, phi) end function rotate_angle -!=============================================================================== -! MAXWELL_SPECTRUM samples an energy from the Maxwell fission distribution based -! on a direct sampling scheme. The probability distribution function for a -! Maxwellian is given as p(x) = 2/(T*sqrt(pi))*sqrt(x/T)*exp(-x/T). This PDF can -! be sampled using rule C64 in the Monte Carlo Sampler LA-9721-MS. -!=============================================================================== - - function maxwell_spectrum(T) result(E_out) - - real(8), intent(in) :: T ! tabulated function of incoming E - real(8) :: E_out ! sampled energy - - real(8) :: r1, r2, r3 ! random numbers - real(8) :: c ! cosine of pi/2*r3 - - r1 = prn() - r2 = prn() - r3 = prn() - - ! determine cosine of pi/2*r - c = cos(PI/TWO*r3) - - ! determine outgoing energy - E_out = -T*(log(r1) + log(r2)*c*c) - - end function maxwell_spectrum - -!=============================================================================== -! WATT_SPECTRUM samples the outgoing energy from a Watt energy-dependent fission -! spectrum. Although fitted parameters exist for many nuclides, generally the -! continuous tabular distributions (LAW 4) should be used in lieu of the Watt -! spectrum. This direct sampling scheme is an unpublished scheme based on the -! original Watt spectrum derivation (See F. Brown's MC lectures). -!=============================================================================== - - function watt_spectrum(a, b) result(E_out) - - real(8), intent(in) :: a ! Watt parameter a - real(8), intent(in) :: b ! Watt parameter b - real(8) :: E_out ! energy of emitted neutron - - real(8) :: w ! sampled from Maxwellian - - w = maxwell_spectrum(a) - E_out = w + a*a*b/4. + (TWO*prn() - ONE)*sqrt(a*a*b*w) - - end function watt_spectrum - !=============================================================================== ! FADDEEVA the Faddeeva function, using Stephen Johnson's implementation !=============================================================================== - function faddeeva(z) result(wv) - complex(C_DOUBLE_COMPLEX), intent(in) :: z ! The point to evaluate Z at - complex(8) :: wv ! The resulting w(z) value + function faddeeva(z) result(wv) bind(C) + complex(C_DOUBLE_COMPLEX), intent(in) :: z ! The point to evaluate Z at + complex(C_DOUBLE_COMPLEX) :: wv ! The resulting w(z) value real(C_DOUBLE) :: relerr ! Target relative error in inner loop of MIT ! Faddeeva @@ -853,10 +182,10 @@ contains end function faddeeva - recursive function w_derivative(z, order) result(wv) + recursive function w_derivative(z, order) result(wv) bind(C) complex(C_DOUBLE_COMPLEX), intent(in) :: z ! The point to evaluate Z at - integer, intent(in) :: order - complex(8) :: wv ! The resulting w(z) value + integer(C_INT), intent(in) :: order + complex(C_DOUBLE_COMPLEX) :: wv ! The resulting w(z) value select case(order) case (0) @@ -869,62 +198,4 @@ contains end select end function w_derivative -!=============================================================================== -! 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) ... -!=============================================================================== - - subroutine broaden_wmp_polynomials(E, dopp, n, factors) - real(8), intent(in) :: E ! Energy to evaluate at - real(8), intent(in) :: dopp ! sqrt(atomic weight ratio / kT), - ! kT given in eV. - integer, intent(in) :: n ! number of components to polynomial - real(8), intent(out):: factors(n) ! output leading coefficient - - integer :: i - - real(8) :: sqrtE ! sqrt(energy) - real(8) :: beta ! sqrt(atomic weight ratio * E / kT) - real(8) :: half_inv_dopp2 ! 0.5 / dopp**2 - real(8) :: quarter_inv_dopp4 ! 0.25 / dopp**4 - real(8) :: erf_beta ! error function of beta - real(8) :: exp_m_beta2 ! exp(-beta**2) - - sqrtE = sqrt(E) - beta = sqrtE * dopp - half_inv_dopp2 = HALF / dopp**2 - quarter_inv_dopp4 = half_inv_dopp2**2 - - if (beta > 6.0_8) then - ! Save time, ERF(6) is 1 to machine precision. - ! beta/sqrtpi*exp(-beta**2) is also approximately 1 machine epsilon. - erf_beta = ONE - exp_m_beta2 = ZERO - else - erf_beta = erf(beta) - exp_m_beta2 = exp(-beta**2) - end if - - ! Assume that, for sure, we'll use a second order (1/E, 1/V, const) - ! fit, and no less. - - factors(1) = erf_beta / E - factors(2) = ONE / sqrtE - factors(3) = factors(1) * (half_inv_dopp2 + E) & - + exp_m_beta2 / (beta * SQRT_PI) - - ! Perform recursive broadening of high order components - do i = 1, n-3 - if (i /= 1) then - factors(i+3) = -factors(i-1) * (i - ONE) * i * quarter_inv_dopp4 & - + factors(i+1) * (E + (ONE + TWO * i) * half_inv_dopp2) - else - ! Although it's mathematically identical, factors(0) will contain - ! nothing, and we don't want to have to worry about memory. - factors(i+3) = factors(i+1)*(E + (ONE + TWO * i) * half_inv_dopp2) - end if - end do - end subroutine broaden_wmp_polynomials - end module math diff --git a/src/math_functions.cpp b/src/math_functions.cpp new file mode 100644 index 000000000..159e0a4e2 --- /dev/null +++ b/src/math_functions.cpp @@ -0,0 +1,693 @@ +#include "math_functions.h" + +namespace openmc { + +//============================================================================== +// Mathematical methods +//============================================================================== + +double normal_percentile_c(double p) { + constexpr double p_low = 0.02425; + constexpr double a[6] = {-3.969683028665376e1, 2.209460984245205e2, + -2.759285104469687e2, 1.383577518672690e2, + -3.066479806614716e1, 2.506628277459239e0}; + constexpr double b[5] = {-5.447609879822406e1, 1.615858368580409e2, + -1.556989798598866e2, 6.680131188771972e1, + -1.328068155288572e1}; + constexpr double c[6] = {-7.784894002430293e-3, -3.223964580411365e-1, + -2.400758277161838, -2.549732539343734, + 4.374664141464968, 2.938163982698783}; + constexpr double d[4] = {7.784695709041462e-3, 3.224671290700398e-1, + 2.445134137142996, 3.754408661907416}; + + // The rational approximation used here is from an unpublished work at + // http://home.online.no/~pjacklam/notes/invnorm/ + + double z; + double q; + + if (p < p_low) { + // Rational approximation for lower region. + + q = std::sqrt(-2.0 * std::log(p)); + z = (((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5]) / + ((((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0); + + } else if (p <= 1.0 - p_low) { + // Rational approximation for central region + q = p - 0.5; + double r = q * q; + z = (((((a[0]*r + a[1])*r + a[2])*r + a[3])*r + a[4])*r + a[5])*q / + (((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1.0); + + } else { + // Rational approximation for upper region + + q = std::sqrt(-2.0*std::log(1.0 - p)); + z = -(((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5]) / + ((((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0); + } + + // Refinement based on Newton's method + + z = z - (0.5 * std::erfc(-z / std::sqrt(2.0)) - p) * std::sqrt(2.0 * PI) * + std::exp(0.5 * z * z); + + return z; + +} + + +double t_percentile_c(double p, int df){ + double t; + + if (df == 1) { + // For one degree of freedom, the t-distribution becomes a Cauchy + // distribution whose cdf we can invert directly + + t = std::tan(PI*(p - 0.5)); + } else if (df == 2) { + // For two degrees of freedom, the cdf is given by 1/2 + x/(2*sqrt(x^2 + + // 2)). This can be directly inverted to yield the solution below + + t = 2.0 * std::sqrt(2.0)*(p - 0.5) / + std::sqrt(1. - 4. * std::pow(p - 0.5, 2.)); + } else { + // This approximation is from E. Olusegun George and Meenakshi Sivaram, "A + // modification of the Fisher-Cornish approximation for the student t + // percentiles," Communication in Statistics - Simulation and Computation, + // 16 (4), pp. 1123-1132 (1987). + double n = df; + double k = 1. / (n - 2.); + double z = normal_percentile_c(p); + double z2 = z * z; + t = std::sqrt(n * k) * (z + (z2 - 3.) * z * k / 4. + ((5. * z2 - 56.) * z2 + + 75.) * z * k * k / 96. + (((z2 - 27.) * 3. * z2 + 417.) * z2 - 315.) * + z * k * k * k / 384.); + } + + return t; +} + + +void calc_pn_c(int n, double x, double pnx[]) { + pnx[0] = 1.; + if (n >= 1) { + pnx[1] = x; + } + + // Use recursion relation to build the higher orders + for (int l = 1; l < n; l ++) { + pnx[l + 1] = ((2 * l + 1) * x * pnx[l] - l * pnx[l - 1]) / (l + 1); + } +} + + +double evaluate_legendre_c(int n, const double data[], double x) { + double pnx[n + 1]; + double val = 0.0; + calc_pn_c(n, x, pnx); + for (int l = 0; l <= n; l++) { + val += (l + 0.5) * data[l] * pnx[l]; + } + return val; +} + + +void calc_rn_c(int n, const double uvw[3], double rn[]){ + // rn[] is assumed to have already been allocated to the correct size + + // Store the cosine of the polar angle and the azimuthal angle + double w = uvw[2]; + double phi; + if (uvw[0] == 0.) { + phi = 0.; + } else { + phi = std::atan2(uvw[1], uvw[0]); + } + + // Store the shorthand of 1-w * w + double w2m1 = 1. - w * w; + + // Now evaluate the spherical harmonics function + rn[0] = 1.; + int i = 0; + for (int l = 1; l <= n; l++) { + // Set the index to the start of this order + i += 2 * (l - 1) + 1; + + // Now evaluate each + switch (l) { + case 1: + // l = 1, m = -1 + rn[i] = -(std::sqrt(w2m1) * std::sin(phi)); + // l = 1, m = 0 + rn[i + 1] = w; + // l = 1, m = 1 + rn[i + 2] = -(std::sqrt(w2m1) * std::cos(phi)); + break; + case 2: + // l = 2, m = -2 + rn[i] = 0.288675134594813 * (-3. * w * w + 3.) * std::sin(2. * phi); + // l = 2, m = -1 + rn[i + 1] = -(1.73205080756888 * w*std::sqrt(w2m1) * std::sin(phi)); + // l = 2, m = 0 + rn[i + 2] = 1.5 * w * w - 0.5; + // l = 2, m = 1 + rn[i + 3] = -(1.73205080756888 * w*std::sqrt(w2m1) * std::cos(phi)); + // l = 2, m = 2 + rn[i + 4] = 0.288675134594813 * (-3. * w * w + 3.) * std::cos(2. * phi); + break; + case 3: + // l = 3, m = -3 + rn[i] = -(0.790569415042095 * std::pow(w2m1, 1.5) * std::sin(3. * phi)); + // l = 3, m = -2 + rn[i + 1] = 1.93649167310371 * w*(w2m1) * std::sin(2.*phi); + // l = 3, m = -1 + rn[i + 2] = -(0.408248290463863*std::sqrt(w2m1)*((7.5)*w * w - 3./2.) * + std::sin(phi)); + // l = 3, m = 0 + rn[i + 3] = 2.5 * std::pow(w, 3) - 1.5 * w; + // l = 3, m = 1 + rn[i + 4] = -(0.408248290463863*std::sqrt(w2m1)*((7.5)*w * w - 3./2.) * + std::cos(phi)); + // l = 3, m = 2 + rn[i + 5] = 1.93649167310371 * w*(w2m1) * std::cos(2.*phi); + // l = 3, m = 3 + rn[i + 6] = -(0.790569415042095 * std::pow(w2m1, 1.5) * std::cos(3.* phi)); + break; + case 4: + // l = 4, m = -4 + rn[i] = 0.739509972887452 * (w2m1 * w2m1) * std::sin(4.0*phi); + // l = 4, m = -3 + rn[i + 1] = -(2.09165006633519 * w * std::pow(w2m1, 1.5) * std::sin(3.* phi)); + // l = 4, m = -2 + rn[i + 2] = 0.074535599249993 * (w2m1)*(52.5 * w * w - 7.5) * std::sin(2. *phi); + // l = 4, m = -1 + rn[i + 3] = -(0.316227766016838*std::sqrt(w2m1)*(17.5 * std::pow(w, 3) - 7.5 * w) * + std::sin(phi)); + // l = 4, m = 0 + rn[i + 4] = 4.375 * std::pow(w, 4) - 3.75 * w * w + 0.375; + // l = 4, m = 1 + rn[i + 5] = -(0.316227766016838*std::sqrt(w2m1)*(17.5 * std::pow(w, 3) - 7.5*w) * + std::cos(phi)); + // l = 4, m = 2 + rn[i + 6] = 0.074535599249993 * (w2m1)*(52.5*w * w - 7.5) * std::cos(2.*phi); + // l = 4, m = 3 + rn[i + 7] = -(2.09165006633519 * w * std::pow(w2m1, 1.5) * std::cos(3.* phi)); + // l = 4, m = 4 + rn[i + 8] = 0.739509972887452 * w2m1 * w2m1 * std::cos(4.0*phi); + break; + case 5: + // l = 5, m = -5 + rn[i] = -(0.701560760020114 * std::pow(w2m1, 2.5) * std::sin(5.0 * phi)); + // l = 5, m = -4 + rn[i + 1] = 2.21852991866236 * w * w2m1 * w2m1 * std::sin(4.0 * phi); + // l = 5, m = -3 + rn[i + 2] = -(0.00996023841111995 * std::pow(w2m1, 1.5) * + ((945.0 /2.)* w * w - 52.5) * std::sin(3.*phi)); + // l = 5, m = -2 + rn[i + 3] = 0.0487950036474267 * (w2m1) + * ((315.0/2.)* std::pow(w, 3) - 52.5 * w) * std::sin(2.*phi); + // l = 5, m = -1 + rn[i + 4] = -(0.258198889747161*std::sqrt(w2m1) * + (39.375 * std::pow(w, 4) - 105.0/4.0 * w * w + 15.0/8.0) * std::sin(phi)); + // l = 5, m = 0 + rn[i + 5] = 7.875 * std::pow(w, 5) - 8.75 * std::pow(w, 3) + 1.875 * w; + // l = 5, m = 1 + rn[i + 6] = -(0.258198889747161 * std::sqrt(w2m1)* + (39.375 * std::pow(w, 4) - 105.0/4.0 * w * w + 15.0/8.0) * std::cos(phi)); + // l = 5, m = 2 + rn[i + 7] = 0.0487950036474267 * (w2m1) * + ((315.0 / 2.) * std::pow(w, 3) - 52.5*w) * std::cos(2.*phi); + // l = 5, m = 3 + rn[i + 8] = -(0.00996023841111995 * std::pow(w2m1, 1.5) * + ((945.0 / 2.) * w * w - 52.5) * std::cos(3.*phi)); + // l = 5, m = 4 + rn[i + 9] = 2.21852991866236 * w * w2m1 * w2m1 * std::cos(4.0*phi); + // l = 5, m = 5 + rn[i + 10] = -(0.701560760020114 * std::pow(w2m1, 2.5) * std::cos(5.0* phi)); + break; + case 6: + // l = 6, m = -6 + rn[i] = 0.671693289381396 * std::pow(w2m1, 3) * std::sin(6.0*phi); + // l = 6, m = -5 + rn[i + 1] = -(2.32681380862329 * w*std::pow(w2m1, 2.5) * std::sin(5.0*phi)); + // l = 6, m = -4 + rn[i + 2] = 0.00104990131391452 * w2m1 * w2m1 * + ((10395.0/2.) * w * w - 945.0/2.) * std::sin(4.0 * phi); + // l = 6, m = -3 + rn[i + 3] = -(0.00575054632785295 * std::pow(w2m1, 1.5) * + ((3465.0/2.) * std::pow(w, 3) - 945.0/2.*w) * std::sin(3.*phi)); + // l = 6, m = -2 + rn[i + 4] = 0.0345032779671177 * (w2m1) * + ((3465.0/8.0)* std::pow(w, 4) - 945.0/4.0 * w * w + 105.0/8.0) * + std::sin(2. * phi); + // l = 6, m = -1 + rn[i + 5] = -(0.218217890235992*std::sqrt(w2m1) * + ((693.0/8.0)* std::pow(w, 5)- 315.0/4.0 * std::pow(w, 3) + (105.0/8.0)*w) * + std::sin(phi)); + // l = 6, m = 0 + rn[i + 6] = 14.4375 * std::pow(w, 6) - 19.6875 * std::pow(w, 4) + 6.5625 * w * w - + 0.3125; + // l = 6, m = 1 + rn[i + 7] = -(0.218217890235992*std::sqrt(w2m1) * + ((693.0/8.0)* std::pow(w, 5)- 315.0/4.0 * std::pow(w, 3) + (105.0/8.0)*w) * + std::cos(phi)); + // l = 6, m = 2 + rn[i + 8] = 0.0345032779671177 * w2m1 * + ((3465.0/8.0)* std::pow(w, 4) -945.0/4.0 * w * w + 105.0/8.0) * + std::cos(2.*phi); + // l = 6, m = 3 + rn[i + 9] = -(0.00575054632785295 * std::pow(w2m1, 1.5) * + ((3465.0/2.) * std::pow(w, 3) - 945.0/2.*w) * std::cos(3.*phi)); + // l = 6, m = 4 + rn[i + 10] = 0.00104990131391452 * w2m1 * w2m1 * + ((10395.0/2.)*w * w - 945.0/2.) * std::cos(4.0*phi); + // l = 6, m = 5 + rn[i + 11] = -(2.32681380862329 * w * std::pow(w2m1, 2.5) * std::cos(5.0*phi)); + // l = 6, m = 6 + rn[i + 12] = 0.671693289381396 * std::pow(w2m1, 3) * std::cos(6.0*phi); + break; + case 7: + // l = 7, m = -7 + rn[i] = -(0.647259849287749 * std::pow(w2m1, 3.5) * std::sin(7.0*phi)); + // l = 7, m = -6 + rn[i + 1] = 2.42182459624969 * w*std::pow(w2m1, 3) * std::sin(6.0*phi); + // l = 7, m = -5 + rn[i + 2] = -(9.13821798555235e-5*std::pow(w2m1, 2.5) * + ((135135.0/2.)*w * w - 10395.0/2.) * std::sin(5.0*phi)); + // l = 7, m = -4 + rn[i + 3] = 0.000548293079133141 * w2m1 * w2m1 * + ((45045.0/2.)*std::pow(w, 3) - 10395.0/2.*w) * std::sin(4.0*phi); + // l = 7, m = -3 + rn[i + 4] = -(0.00363696483726654 * std::pow(w2m1, 1.5) * + ((45045.0/8.0)* std::pow(w, 4) - 10395.0/4.0 * w * w + 945.0/8.0) * + std::sin(3.*phi)); + // l = 7, m = -2 + rn[i + 5] = 0.025717224993682 * (w2m1) * + ((9009.0/8.0)* std::pow(w, 5) -3465.0/4.0 * std::pow(w, 3) + (945.0/8.0)*w) * + std::sin(2.*phi); + // l = 7, m = -1 + rn[i + 6] = -(0.188982236504614*std::sqrt(w2m1) * + ((3003.0/16.0)* std::pow(w, 6) - 3465.0/16.0 * std::pow(w, 4) + + (945.0/16.0)*w * w - 35.0/16.0) * std::sin(phi)); + // l = 7, m = 0 + rn[i + 7] = 26.8125 * std::pow(w, 7) - 43.3125 * std::pow(w, 5) + 19.6875 * std::pow(w, 3) - + 2.1875 * w; + // l = 7, m = 1 + rn[i + 8] = -(0.188982236504614*std::sqrt(w2m1) * ((3003.0/16.0) * std::pow(w, 6) - + 3465.0/16.0 * std::pow(w, 4) + (945.0/16.0)*w * w - 35.0/16.0) * std::cos(phi)); + // l = 7, m = 2 + rn[i + 9] = 0.025717224993682 * (w2m1) * ((9009.0/8.0)* std::pow(w, 5) - + 3465.0/4.0 * std::pow(w, 3) + (945.0/8.0)*w) * std::cos(2.*phi); + // l = 7, m = 3 + rn[i + 10] = -(0.00363696483726654 * std::pow(w2m1, 1.5) * + ((45045.0/8.0)* std::pow(w, 4) - 10395.0/4.0 * w * w + 945.0/8.0) * + std::cos(3.*phi)); + // l = 7, m = 4 + rn[i + 11] = 0.000548293079133141 * w2m1 * w2m1 * + ((45045.0/2.)*std::pow(w, 3) - 10395.0/2.*w) * std::cos(4.0*phi); + // l = 7, m = 5 + rn[i + 12] = -(9.13821798555235e-5*std::pow(w2m1, 2.5) * + ((135135.0/2.)*w * w - 10395.0/2.) * std::cos(5.0*phi)); + // l = 7, m = 6 + rn[i + 13] = 2.42182459624969 * w*std::pow(w2m1, 3) * std::cos(6.0*phi); + // l = 7, m = 7 + rn[i + 14] = -(0.647259849287749 * std::pow(w2m1, 3.5) * std::cos(7.0*phi)); + break; + case 8: + // l = 8, m = -8 + rn[i] = 0.626706654240044 * std::pow(w2m1, 4) * std::sin(8.0*phi); + // l = 8, m = -7 + rn[i + 1] = -(2.50682661696018 * w*std::pow(w2m1, 3.5) * std::sin(7.0*phi)); + // l = 8, m = -6 + rn[i + 2] = 6.77369783729086e-6*std::pow(w2m1, 3)* + ((2027025.0/2.)*w * w - 135135.0/2.) * std::sin(6.0*phi); + // l = 8, m = -5 + rn[i + 3] = -(4.38985792528482e-5*std::pow(w2m1, 2.5) * + ((675675.0/2.)*std::pow(w, 3) - 135135.0/2.*w) * std::sin(5.0*phi)); + // l = 8, m = -4 + rn[i + 4] = 0.000316557156832328 * w2m1 * w2m1 * + ((675675.0/8.0)* std::pow(w, 4) - 135135.0/4.0 * w * w + 10395.0/8.0) * + std::sin(4.0*phi); + // l = 8, m = -3 + rn[i + 5] = -(0.00245204119306875 * std::pow(w2m1, 1.5) * ((135135.0/8.0) * + std::pow(w, 5) - 45045.0/4.0 * std::pow(w, 3) + (10395.0/8.0)*w) * std::sin(3.*phi)); + // l = 8, m = -2 + rn[i + 6] = 0.0199204768222399 * (w2m1) * + ((45045.0/16.0)* std::pow(w, 6)- 45045.0/16.0 * std::pow(w, 4) + + (10395.0/16.0)*w * w - 315.0/16.0) * std::sin(2.*phi); + // l = 8, m = -1 + rn[i + 7] = -(0.166666666666667*std::sqrt(w2m1) * + ((6435.0/16.0)* std::pow(w, 7) - 9009.0/16.0 * std::pow(w, 5) + + (3465.0/16.0)*std::pow(w, 3) - 315.0/16.0 * w) * std::sin(phi)); + // l = 8, m = 0 + rn[i + 8] = 50.2734375 * std::pow(w, 8) - 93.84375 * std::pow(w, 6) + 54.140625 * + std::pow(w, 4) - 9.84375 * w * w + 0.2734375; + // l = 8, m = 1 + rn[i + 9] = -(0.166666666666667*std::sqrt(w2m1) * + ((6435.0/16.0)* std::pow(w, 7) - 9009.0/16.0 * std::pow(w, 5) + + (3465.0/16.0)*std::pow(w, 3) - 315.0/16.0 * w) * std::cos(phi)); + // l = 8, m = 2 + rn[i + 10] = 0.0199204768222399 * (w2m1)*((45045.0/16.0)* std::pow(w, 6)- + 45045.0/16.0 * std::pow(w, 4) + (10395.0/16.0)*w * w - + 315.0/16.0) * std::cos(2.*phi); + // l = 8, m = 3 + rn[i + 11] = -(0.00245204119306875 * std::pow(w2m1, 1.5)* + ((135135.0/8.0) * std::pow(w, 5) - 45045.0/4.0 * std::pow(w, 3) + + (10395.0/8.0)*w) * std::cos(3.*phi)); + // l = 8, m = 4 + rn[i + 12] = 0.000316557156832328 * w2m1 * w2m1*((675675.0/8.0)* std::pow(w, 4) - + 135135.0/4.0 * w * w + 10395.0/8.0) * std::cos(4.0*phi); + // l = 8, m = 5 + rn[i + 13] = -(4.38985792528482e-5*std::pow(w2m1, 2.5)*((675675.0/2.)*std::pow(w, 3) - + 135135.0/2.*w) * std::cos(5.0*phi)); + // l = 8, m = 6 + rn[i + 14] = 6.77369783729086e-6*std::pow(w2m1, 3)*((2027025.0/2.)*w * w - + 135135.0/2.) * std::cos(6.0*phi); + // l = 8, m = 7 + rn[i + 15] = -(2.50682661696018 * w*std::pow(w2m1, 3.5) * std::cos(7.0*phi)); + // l = 8, m = 8 + rn[i + 16] = 0.626706654240044 * std::pow(w2m1, 4) * std::cos(8.0*phi); + break; + case 9: + // l = 9, m = -9 + rn[i] = -(0.609049392175524 * std::pow(w2m1, 4.5) * std::sin(9.0 * phi)); + // l = 9, m = -8 + rn[i + 1] = 2.58397773170915 * w*std::pow(w2m1, 4) * std::sin(8.0 * phi); + // l = 9, m = -7 + rn[i + 2] = -(4.37240315267812e-7*std::pow(w2m1, 3.5) * + ((34459425.0/2.)*w * w - 2027025.0/2.) * std::sin(7.0 * phi)); + // l = 9, m = -6 + rn[i + 3] = 3.02928976464514e-6*std::pow(w2m1, 3)* + ((11486475.0/2.)*std::pow(w, 3) - 2027025.0/2.*w) * std::sin(6.0 * phi); + // l = 9, m = -5 + rn[i + 4] = -(2.34647776186144e-5*std::pow(w2m1, 2.5) * + ((11486475.0/8.0)* std::pow(w, 4) - 2027025.0 / 4.0 * w * w + + 135135.0/8.0) * std::sin(5.0 * phi)); + // l = 9, m = -4 + rn[i + 5] = 0.000196320414650061 * w2m1 * w2m1*((2297295.0/8.0)* std::pow(w, 5) - + 675675.0/4.0 * std::pow(w, 3) + (135135.0/8.0)*w) * std::sin(4.0*phi); + // l = 9, m = -3 + rn[i + 6] = -(0.00173385495536766 * std::pow(w2m1, 1.5) * + ((765765.0/16.0)* std::pow(w, 6) - 675675.0/16.0 * std::pow(w, 4) + + (135135.0/16.0)*w * w - 3465.0/16.0) * std::sin(3. * phi)); + // l = 9, m = -2 + rn[i + 7] = 0.0158910431540932 * (w2m1)*((109395.0/16.0)* std::pow(w, 7)- + 135135.0/16.0 * std::pow(w, 5) + (45045.0/16.0)*std::pow(w, 3) - + 3465.0/16.0 * w) * std::sin(2. * phi); + // l = 9, m = -1 + rn[i + 8] = -(0.149071198499986*std::sqrt(w2m1)*((109395.0/128.0)* std::pow(w, 8) - + 45045.0/32.0 * std::pow(w, 6) + (45045.0/64.0)* std::pow(w, 4) - + 3465.0/32.0 * w * w + 315.0/128.0) * std::sin(phi)); + // l = 9, m = 0 + rn[i + 9] = 94.9609375 * std::pow(w, 9) - 201.09375 * std::pow(w, 7) + + 140.765625 * std::pow(w, 5)- 36.09375 * std::pow(w, 3) + 2.4609375 * w; + // l = 9, m = 1 + rn[i + 10] = -(0.149071198499986*std::sqrt(w2m1)*((109395.0/128.0)* std::pow(w, 8) - + 45045.0/32.0 * std::pow(w, 6) + (45045.0/64.0)* std::pow(w, 4) - + 3465.0/32.0 * w * w + 315.0/128.0) * std::cos(phi)); + // l = 9, m = 2 + rn[i + 11] = 0.0158910431540932 * (w2m1)*((109395.0/16.0)* std::pow(w, 7) - + 135135.0/16.0 * std::pow(w, 5) + (45045.0/16.0)*std::pow(w, 3) - + 3465.0/ 16.0 * w) * std::cos(2. * phi); + // l = 9, m = 3 + rn[i + 12] = -(0.00173385495536766 * std::pow(w2m1, 1.5)*((765765.0/16.0) * + std::pow(w, 6) - 675675.0/16.0 * std::pow(w, 4) + + (135135.0/16.0)* w * w - 3465.0/16.0)* std::cos(3. * phi)); + // l = 9, m = 4 + rn[i + 13] = 0.000196320414650061 * w2m1 * w2m1*((2297295.0/8.0) * std::pow(w, 5) - + 675675.0/4.0 * std::pow(w, 3) + (135135.0/8.0)*w) * std::cos(4.0 * phi); + // l = 9, m = 5 + rn[i + 14] = -(2.34647776186144e-5*std::pow(w2m1, 2.5)*((11486475.0/8.0) * + std::pow(w, 4) - 2027025.0/4.0 * w * w + 135135.0/8.0) * + std::cos(5.0 * phi)); + // l = 9, m = 6 + rn[i + 15] = 3.02928976464514e-6*std::pow(w2m1, 3)*((11486475.0/2.)*std::pow(w, 3) - + 2027025.0/2. * w) * std::cos(6.0 * phi); + // l = 9, m = 7 + rn[i + 16] = -(4.37240315267812e-7*std::pow(w2m1, 3.5)* + ((34459425.0/2.) * w * w - 2027025.0/2.) * std::cos(7.0 * phi)); + // l = 9, m = 8 + rn[i + 17] = 2.58397773170915 * w*std::pow(w2m1, 4) * std::cos(8.0 * phi); + // l = 9, m = 9 + rn[i + 18] = -(0.609049392175524 * std::pow(w2m1, 4.5) * std::cos(9.0 * phi)); + break; + case 10: + // l = 10, m = -10 + rn[i] = 0.593627917136573 * std::pow(w2m1, 5) * std::sin(10.0 * phi); + // l = 10, m = -9 + rn[i + 1] = -(2.65478475211798 * w * std::pow(w2m1, 4.5) * std::sin(9.0 * phi)); + // l = 10, m = -8 + rn[i + 2] = 2.49953651452314e-8 * std::pow(w2m1, 4) * + ((654729075.0/2.) * w * w - 34459425.0/2.) * std::sin(8.0 * phi); + // l = 10, m = -7 + rn[i + 3] = -(1.83677671621093e-7*std::pow(w2m1, 3.5)* + ((218243025.0/2.)*std::pow(w, 3) - 34459425.0/2.*w) * + std::sin(7.0 * phi)); + // l = 10, m = -6 + rn[i + 4] = 1.51464488232257e-6*std::pow(w2m1, 3)*((218243025.0/8.0)* std::pow(w, 4) - + 34459425.0/4.0 * w * w + 2027025.0/8.0) * std::sin(6.0 * phi); + // l = 10, m = -5 + rn[i + 5] = -(1.35473956745817e-5*std::pow(w2m1, 2.5)* + ((43648605.0/8.0)* std::pow(w, 5) - 11486475.0/4.0 * std::pow(w, 3) + + (2027025.0/8.0)*w) * std::sin(5.0 * phi)); + // l = 10, m = -4 + rn[i + 6] = 0.000128521880085575 * w2m1 * w2m1*((14549535.0/16.0)* std::pow(w, 6) - + 11486475.0/16.0 * std::pow(w, 4) + (2027025.0/16.0)*w * w - + 45045.0/16.0) * std::sin(4.0 * phi); + // l = 10, m = -3 + rn[i + 7] = -(0.00127230170115096 * std::pow(w2m1, 1.5)* + ((2078505.0/16.0)* std::pow(w, 7) - 2297295.0/16.0 * std::pow(w, 5) + + (675675.0/16.0)*std::pow(w, 3) - 45045.0/16.0 * w) * std::sin(3. * phi)); + // l = 10, m = -2 + rn[i + 8] = 0.012974982402692 * (w2m1)*((2078505.0/128.0)* std::pow(w, 8) - + 765765.0/32.0 * std::pow(w, 6) + (675675.0/64.0)* std::pow(w, 4) - + 45045.0/32.0 * w * w + 3465.0/128.0) * std::sin(2. * phi); + // l = 10, m = -1 + rn[i + 9] = -(0.134839972492648*std::sqrt(w2m1)*((230945.0/128.0)* std::pow(w, 9) - + 109395.0/32.0 * std::pow(w, 7) + (135135.0/64.0)* std::pow(w, 5) - + 15015.0/32.0 * std::pow(w, 3) + (3465.0/128.0)*w) * std::sin(phi)); + // l = 10, m = 0 + rn[i + 10] = 180.42578125 * std::pow(w, 10) - 427.32421875 * std::pow(w, 8) +351.9140625 + * std::pow(w, 6) - 117.3046875 * std::pow(w, 4) + 13.53515625 * w * w -0.24609375; + // l = 10, m = 1 + rn[i + 11] = -(0.134839972492648*std::sqrt(w2m1)*((230945.0/128.0)* std::pow(w, 9) - + 109395.0/32.0 * std::pow(w, 7) + (135135.0/64.0)* std::pow(w, 5) -15015.0/ + 32.0 * std::pow(w, 3) + (3465.0/128.0)*w) * std::cos(phi)); + // l = 10, m = 2 + rn[i + 12] = 0.012974982402692 * (w2m1)*((2078505.0/128.0)* std::pow(w, 8) - + 765765.0/32.0 * std::pow(w, 6) + (675675.0/64.0)* std::pow(w, 4) - + 45045.0/32.0 * w * w + 3465.0/128.0) * std::cos(2. * phi); + // l = 10, m = 3 + rn[i + 13] = -(0.00127230170115096 * std::pow(w2m1, 1.5)* + ((2078505.0/16.0)* std::pow(w, 7) - 2297295.0/16.0 * std::pow(w, 5) + + (675675.0/16.0)*std::pow(w, 3) - 45045.0/16.0 * w) * std::cos(3. * phi)); + // l = 10, m = 4 + rn[i + 14] = 0.000128521880085575 * w2m1 * w2m1*((14549535.0/16.0)* std::pow(w, 6) - + 11486475.0/16.0 * std::pow(w, 4) + (2027025.0/16.0) * w * w - + 45045.0/16.0) * std::cos(4.0 * phi); + // l = 10, m = 5 + rn[i + 15] = -(1.35473956745817e-5*std::pow(w2m1, 2.5)* + ((43648605.0/8.0)* std::pow(w, 5) - 11486475.0/4.0 * std::pow(w, 3) + + (2027025.0/8.0)*w) * std::cos(5.0 * phi)); + // l = 10, m = 6 + rn[i + 16] = 1.51464488232257e-6*std::pow(w2m1, 3)*((218243025.0/8.0)* std::pow(w, 4) - + 34459425.0/4.0 * w * w + 2027025.0/8.0) * std::cos(6.0 * phi); + // l = 10, m = 7 + rn[i + 17] = -(1.83677671621093e-7*std::pow(w2m1, 3.5) * + ((218243025.0/2.)*std::pow(w, 3) - 34459425.0/2.*w) * std::cos(7.0 * phi)); + // l = 10, m = 8 + rn[i + 18] = 2.49953651452314e-8*std::pow(w2m1, 4)* + ((654729075.0/2.)*w * w - 34459425.0/2.) * std::cos(8.0 * phi); + // l = 10, m = 9 + rn[i + 19] = -(2.65478475211798 * w*std::pow(w2m1, 4.5) * std::cos(9.0 * phi)); + // l = 10, m = 10 + rn[i + 20] = 0.593627917136573 * std::pow(w2m1, 5) * std::cos(10.0 * phi); + } + } +} + + +void calc_zn_c(int n, double rho, double phi, double zn[]) { + // =========================================================================== + // 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) + + double sin_phi = std::sin(phi); + double cos_phi = std::cos(phi); + + double sin_phi_vec[n + 1]; // Sin[n * phi] + double cos_phi_vec[n + 1]; // Cos[n * 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 (int 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 (int i = 0; i <= n; i++) { + sin_phi_vec[i] *= sin_phi; + } + + // =========================================================================== + // Calculate R_pq(rho) + double zn_mat[n + 1][n + 1]; // Matrix forms of the coefficients which are + // easier to work with + + // Fill the main diagonal first (Eq 3.9 in Chong) + for (int p = 0; p <= n; p++) { + zn_mat[p][p] = std::pow(rho, p); + } + + // Fill the 2nd diagonal (Eq 3.10 in Chong) + for (int 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 (int p = 4; p <= n; p++) { + double k2 = 2 * p * (p - 1) * (p - 2); + for (int q = p - 4; q >= 0; q -= 2) { + double k1 = ((p + q) * (p - q) * (p - 2)) / 2.; + double k3 = -q * q * (p - 1) - p * (p - 1) * (p - 2); + double k4 = (-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)... ] + int i = 0; + for (int p = 0; p <= n; p++) { + for (int q = -p; q <= p; q += 2) { + if (q < 0) { + zn[i] = zn_mat[std::abs(q)][p] * sin_phi_vec[std::abs(q) - 1]; + } else if (q == 0) { + zn[i] = zn_mat[q][p]; + } else { + zn[i] = zn_mat[q][p] * cos_phi_vec[q]; + } + i++; + } + } + +} + + +void rotate_angle_c(double uvw[3], double mu, double* phi) { + // Copy original directional cosines + double u0 = uvw[0]; // original cosine in x direction + double v0 = uvw[1]; // original cosine in y direction + double w0 = uvw[2]; // original cosine in z direction + + // Sample azimuthal angle in [0,2pi) if none provided + double phi_; + if (phi != nullptr) { + phi_ = (*phi); + } else { + phi_ = 2. * PI * prn(); + } + + // Precompute factors to save flops + double sinphi = std::sin(phi_); + double cosphi = std::cos(phi_); + double a = std::sqrt(std::fmax(0., 1. - mu * mu)); + double b = std::sqrt(std::fmax(0., 1. - 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) { + uvw[0] = mu * u0 + a * (u0 * w0 * cosphi - v0 * sinphi) / b; + uvw[1] = mu * v0 + a * (v0 * w0 * cosphi + u0 * sinphi) / b; + uvw[2] = mu * w0 - a * b * cosphi; + } else { + b = std::sqrt(1. - v0 * v0); + uvw[0] = mu * u0 + a * (u0 * v0 * cosphi + w0 * sinphi) / b; + uvw[1] = mu * v0 - a * b * cosphi; + uvw[2] = mu * w0 + a * (v0 * w0 * cosphi - u0 * sinphi) / b; + } +} + + +double maxwell_spectrum_c(double T) { + // Set the random numbers + double r1 = prn(); + double r2 = prn(); + double r3 = prn(); + + // determine cosine of pi/2*r + double c = std::cos(PI / 2. * r3); + + // Determine outgoing energy + double E_out = -T * (std::log(r1) + std::log(r2) * c * c); + + return E_out; +} + + +double watt_spectrum_c(double a, double b) { + double w = maxwell_spectrum_c(a); + double E_out = w + 0.25 * a * a * b + (2. * prn() - 1.) * std::sqrt(a * a * b * w); + + return E_out; +} + + +void broaden_wmp_polynomials_c(double E, double dopp, int n, double factors[]) { + // Factors is already pre-allocated + double sqrtE = std::sqrt(E); + double beta = sqrtE * dopp; + double half_inv_dopp2 = 0.5 / (dopp * dopp); + double quarter_inv_dopp4 = half_inv_dopp2 * half_inv_dopp2; + + double erf_beta; // error function of beta + double exp_m_beta2; // exp(-beta**2) + if (beta > 6.0) { + // Save time, ERF(6) is 1 to machine precision. + // beta/sqrtpi*exp(-beta**2) is also approximately 1 machine epsilon. + erf_beta = 1.; + exp_m_beta2 = 0.; + } else { + erf_beta = std::erf(beta); + exp_m_beta2 = std::exp(-beta * beta); + } + + // Assume that, for sure, we'll use a second order (1/E, 1/V, const) + // fit, and no less. + + factors[0] = erf_beta / E; + factors[1] = 1. / sqrtE; + factors[2] = factors[0] * (half_inv_dopp2 + E) + exp_m_beta2 / + (beta * SQRT_PI); + + // Perform recursive broadening of high order components + for (int i = 0; i < n - 3; i++) { + double ip1_dbl = i + 1; + if (i != 0) { + factors[i + 3] = -factors[i - 1] * (ip1_dbl - 1.) * ip1_dbl * + quarter_inv_dopp4 + factors[i + 1] * + (E + (1. + 2. * ip1_dbl) * half_inv_dopp2); + } else { + // Although it's mathematically identical, factors[0] will contain + // nothing, and we don't want to have to worry about memory. + factors[i + 3] = factors[i + 1] * + (E + (1. + 2. * ip1_dbl) * half_inv_dopp2); + } + } +} + +} // namespace openmc diff --git a/src/math_functions.h b/src/math_functions.h new file mode 100644 index 000000000..68e89251e --- /dev/null +++ b/src/math_functions.h @@ -0,0 +1,163 @@ +//! \file math_functions.h +//! A collection of elementary math functions. + +#ifndef MATH_FUNCTIONS_H +#define MATH_FUNCTIONS_H + +#include +#include + +#include "random_lcg.h" + + +namespace openmc { + +//============================================================================== +// Module constants. +//============================================================================== + +// TODO: cmath::M_PI has 3 more digits precision than the Fortran constant we +// use so for now we will reuse the Fortran constant until we are OK with +// modifying test results +extern "C" constexpr double PI {3.1415926535898}; + +extern "C" const double SQRT_PI {std::sqrt(PI)}; + +//============================================================================== +//! Calculate the percentile of the standard normal distribution with a +//! specified probability level. +//! +//! @param p The probability level +//! @return The requested percentile +//============================================================================== + +extern "C" double normal_percentile_c(double p); + +//============================================================================== +//! Calculate the percentile of the Student's t distribution with a specified +//! probability level and number of degrees of freedom. +//! +//! @param p The probability level +//! @param df The degrees of freedom +//! @return The requested percentile +//============================================================================== + +extern "C" double t_percentile_c(double p, int df); + +//============================================================================== +//! Calculate the n-th order Legendre polynomials at the value of x. +//! +//! @param n The maximum order requested +//! @param x The value to evaluate at; x is expected to be within [-1,1] +//! @param pnx The requested Legendre polynomials of order 0 to n (inclusive) +//! evaluated at x. +//============================================================================== + +extern "C" void calc_pn_c(int n, double x, double pnx[]); + +//============================================================================== +//! Find the value of f(x) given a set of Legendre coefficients and the value +//! of x. +//! +//! @param n The maximum order of the expansion +//! @param data The polynomial expansion coefficient data; without the (2l+1)/2 +//! factor. +//! @param x The value to evaluate at; x is expected to be within [-1,1] +//! @return The requested Legendre polynomials of order 0 to n (inclusive) +//! evaluated at x +//============================================================================== + +extern "C" double evaluate_legendre_c(int n, const double data[], double x); + +//============================================================================== +//! Calculate the n-th order real spherical harmonics for a given angle (in +//! terms of (u,v,w)) for all 0<=n and -m<=n<=n. +//! +//! @param n The maximum order requested +//! @param uvw[3] The direction the harmonics are requested at +//! @param rn The requested harmonics of order 0 to n (inclusive) +//! evaluated at uvw. +//============================================================================== + +extern "C" void calc_rn_c(int n, const double uvw[3], double rn[]); + +//============================================================================== +//! Calculate the n-th order modified Zernike polynomial moment for a given +//! angle (rho, theta) location on the unit disk. +//! +//! 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 The requested moments of order 0 to n (inclusive) +//! evaluated at rho and phi. +//============================================================================== + +extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]); + +//============================================================================== +//! Rotate the direction cosines through a polar angle whose cosine is mu and +//! through an azimuthal angle sampled uniformly. +//! +//! This is done with direct sampling rather than rejection sampling as is done +//! in MCNP and Serpent. +//! +//! @param uvw[3] The initial, and final, direction vector +//! @param mu The cosine of angle in lab or CM +//! @param phi The azimuthal angle; will randomly chosen angle if a nullptr +//! is passed +//============================================================================== + +extern "C" void rotate_angle_c(double uvw[3], double mu, double* phi); + +//============================================================================== +//! Samples an energy from the Maxwell fission distribution based on a direct +//! sampling scheme. +//! +//! The probability distribution function for a Maxwellian is given as +//! p(x) = 2/(T*sqrt(pi))*sqrt(x/T)*exp(-x/T). This PDF can be sampled using +//! rule C64 in the Monte Carlo Sampler LA-9721-MS. +//! +//! @param T The tabulated function of the incoming energy +//! @result The sampled outgoing energy +//============================================================================== + +extern "C" double maxwell_spectrum_c(double T); + +//============================================================================== +//! Samples an energy from a Watt energy-dependent fission distribution. +//! +//! Although fitted parameters exist for many nuclides, generally the +//! continuous tabular distributions (LAW 4) should be used in lieu of the Watt +//! spectrum. This direct sampling scheme is an unpublished scheme based on the +//! original Watt spectrum derivation (See F. Brown's MC lectures). +//! +//! @param a Watt parameter a +//! @param b Watt parameter b +//! @result The sampled outgoing energy +//============================================================================== + +extern "C" double watt_spectrum_c(double a, double b); + +//============================================================================== +//! Doppler broadens the windowed multipole curvefit. +//! +//! The curvefit is a polynomial of the form a/E + b/sqrt(E) + c + d sqrt(E)... +//! +//! @param E The energy to evaluate the broadening at +//! @param dopp sqrt(atomic weight ratio / kT) with kT given in eV +//! @param n The number of components to the polynomial +//! @param factors The output leading coefficient +//============================================================================== + +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/src/tallies/tally.F90 b/src/tallies/tally.F90 index 2e1d6184c..ef8b5915c 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -7,7 +7,7 @@ module tally use dict_header, only: EMPTY use error, only: fatal_error use geometry_header - use math, only: t_percentile, calc_pn, calc_rn + use math, only: t_percentile use mesh_header, only: RegularMesh, meshes use message_passing use mgxs_header diff --git a/src/tallies/tally_filter_legendre.F90 b/src/tallies/tally_filter_legendre.F90 index b4ee7b06b..4663df651 100644 --- a/src/tallies/tally_filter_legendre.F90 +++ b/src/tallies/tally_filter_legendre.F90 @@ -51,13 +51,12 @@ contains type(TallyFilterMatch), intent(inout) :: match integer :: i - real(8) :: wgt + real(C_DOUBLE) :: wgt(this % n_bins) - ! TODO: Use recursive formula to calculate higher orders - do i = 0, this % order - wgt = calc_pn(i, p % mu) - call match % bins % push_back(i + 1) - call match % weights % push_back(wgt) + call calc_pn(this % order, p % mu, wgt) + do i = 1, this % n_bins + call match % bins % push_back(i) + call match % weights % push_back(wgt(i)) end do end subroutine get_all_bins diff --git a/src/tallies/tally_filter_sph_harm.F90 b/src/tallies/tally_filter_sph_harm.F90 index bc95baaa3..4a1f43273 100644 --- a/src/tallies/tally_filter_sph_harm.F90 +++ b/src/tallies/tally_filter_sph_harm.F90 @@ -74,30 +74,32 @@ contains integer :: i, j, n integer :: num_nm - real(8) :: wgt - real(8) :: rn(2*this % order + 1) + real(C_DOUBLE) :: wgt(this % order + 1) + real(C_DOUBLE) :: rn(this % n_bins) + + ! Determine cosine term for scatter expansion if necessary + if (this % cosine == COSINE_SCATTER) then + call calc_pn(this % order, p % mu, wgt) + else + wgt = ONE + end if + + ! Find the Rn,m values + call calc_rn(this % order, p % last_uvw, rn) - ! TODO: Use recursive formula to calculate higher orders j = 0 do n = 0, this % order - ! Determine cosine term for scatter expansion if necessary - if (this % cosine == COSINE_SCATTER) then - wgt = calc_pn(n, p % mu) - else - wgt = ONE - end if - ! Calculate n-th order spherical harmonics for (u,v,w) num_nm = 2*n + 1 - rn(1:num_nm) = calc_rn(n, p % last_uvw) ! Append matching (bin,weight) for each moment do i = 1, num_nm j = j + 1 call match % bins % push_back(j) - call match % weights % push_back(wgt * rn(i)) + call match % weights % push_back(wgt(n + 1) * rn(j)) end do end do + end subroutine get_all_bins subroutine to_statepoint(this, filter_group) diff --git a/src/tallies/tally_filter_sptl_legendre.F90 b/src/tallies/tally_filter_sptl_legendre.F90 index 83db9a864..1bfbd0e3b 100644 --- a/src/tallies/tally_filter_sptl_legendre.F90 +++ b/src/tallies/tally_filter_sptl_legendre.F90 @@ -75,20 +75,19 @@ contains type(TallyFilterMatch), intent(inout) :: match integer :: i - real(8) :: wgt - real(8) :: x ! Position on specified axis - real(8) :: x_norm ! Normalized position + real(C_DOUBLE) :: wgt(this % n_bins) + real(C_DOUBLE) :: x ! Position on specified axis + real(C_DOUBLE) :: x_norm ! Normalized position x = p % coord(1) % xyz(this % axis) if (this % min <= x .and. x <= this % max) then ! Calculate normalized position between min and max x_norm = TWO*(x - this % min)/(this % max - this % min) - ONE - ! TODO: Use recursive formula to calculate higher orders - do i = 0, this % order - wgt = calc_pn(i, x_norm) - call match % bins % push_back(i + 1) - call match % weights % push_back(wgt) + call calc_pn(this % order, x_norm, wgt) + do i = 1, this % n_bins + call match % bins % push_back(i) + call match % weights % push_back(wgt(i)) end do end if end subroutine get_all_bins diff --git a/src/tallies/tally_filter_zernike.F90 b/src/tallies/tally_filter_zernike.F90 index fa205b180..d00541613 100644 --- a/src/tallies/tally_filter_zernike.F90 +++ b/src/tallies/tally_filter_zernike.F90 @@ -61,7 +61,7 @@ contains integer :: i real(8) :: x, y, r, theta - real(8) :: zn(this % n_bins) + real(C_DOUBLE) :: zn(this % n_bins) ! Determine normalized (r,theta) positions x = p % coord(1) % xyz(1) - this % x diff --git a/tests/unit_tests/test_math.py b/tests/unit_tests/test_math.py new file mode 100644 index 000000000..6825f6b93 --- /dev/null +++ b/tests/unit_tests/test_math.py @@ -0,0 +1,212 @@ +import numpy as np +import scipy as sp + +import openmc +import openmc.capi + + +def test_t_percentile(): + # Permutations include 1 DoF, 2 DoF, and > 2 DoF + # We will test 5 p-values at 3-DoF values + test_ps = [0.02, 0.4, 0.5, 0.6, 0.98] + test_dfs = [1, 2, 5] + + # The reference solutions come from Scipy + ref_ts = [[sp.stats.t.ppf(p, df) for p in test_ps] for df in test_dfs] + + test_ts = [[openmc.capi.math.t_percentile(p, df) for p in test_ps] + for df in test_dfs] + + # The 5 DoF approximation in openmc.capi.math.t_percentile is off by up to + # 8e-3 from the scipy solution, so test that one separately with looser + # tolerance + assert np.allclose(ref_ts[:-1], test_ts[:-1]) + assert np.allclose(ref_ts[-1], test_ts[-1], atol=1e-2) + + +def test_calc_pn(): + max_order = 10 + test_xs = np.linspace(-1., 1., num=5, endpoint=True) + + # Reference solutions from scipy + ref_vals = np.array([sp.special.eval_legendre(n, test_xs) + for n in range(0, max_order + 1)]) + + test_vals = [] + for x in test_xs: + test_vals.append(openmc.capi.math.calc_pn(max_order, x).tolist()) + + test_vals = np.swapaxes(np.array(test_vals), 0, 1) + + assert np.allclose(ref_vals, test_vals) + + +def test_evaluate_legendre(): + max_order = 10 + # Coefficients are set to 1, but will incorporate the (2l+1)/2 norm factor + # for the reference solution + test_coeffs = [0.5 * (2. * l + 1.) for l in range(max_order + 1)] + test_xs = np.linspace(-1., 1., num=5, endpoint=True) + + ref_vals = np.polynomial.legendre.legval(test_xs, test_coeffs) + + # Set the coefficients back to 1s for the test values since + # evaluate legendre incorporates the (2l+1)/2 term on its own + test_coeffs = [1. for l in range(max_order + 1)] + + test_vals = np.array([openmc.capi.math.evaluate_legendre(test_coeffs, x) + for x in test_xs]) + + assert np.allclose(ref_vals, test_vals) + + +def test_calc_rn(): + max_order = 10 + test_ns = np.array([i for i in range(0, max_order + 1)]) + azi = 0.1 # Longitude + pol = 0.2 # Latitude + test_uvw = np.array([np.sin(pol) * np.cos(azi), + np.sin(pol) * np.sin(azi), + np.cos(pol)]) + + # Reference solutions from the equations + ref_vals = [] + + def coeff(n, m): + return np.sqrt((2. * n + 1) * sp.special.factorial(n - m) / + (sp.special.factorial(n + m))) + + def pnm_bar(n, m, mu): + val = coeff(n, m) + if m != 0: + val *= np.sqrt(2.) + val *= sp.special.lpmv([m], [n], [mu]) + return val[0] + + ref_vals = [] + for n in test_ns: + for m in range(-n, n + 1): + if m < 0: + ylm = pnm_bar(n, np.abs(m), np.cos(pol)) * \ + np.sin(np.abs(m) * azi) + else: + ylm = pnm_bar(n, m, np.cos(pol)) * np.cos(m * azi) + + # Un-normalize for comparison + ylm /= np.sqrt(2. * n + 1.) + ref_vals.append(ylm) + + test_vals = [] + test_vals = openmc.capi.math.calc_rn(max_order, test_uvw) + + assert np.allclose(ref_vals, test_vals) + + +def test_calc_zn(): + n = 10 + rho = 0.5 + phi = 0.5 + + # Reference solution from running the Fortran implementation + ref_vals = np.array([ + 1.00000000e+00, 2.39712769e-01, 4.38791281e-01, + 2.10367746e-01, -5.00000000e-01, 1.35075576e-01, + 1.24686873e-01, -2.99640962e-01, -5.48489101e-01, + 8.84215021e-03, 5.68310892e-02, -4.20735492e-01, + -1.25000000e-01, -2.70151153e-01, -2.60091773e-02, + 1.87022545e-02, -3.42888902e-01, 1.49820481e-01, + 2.74244551e-01, -2.43159131e-02, -2.50357380e-02, + 2.20500013e-03, -1.98908812e-01, 4.07587508e-01, + 4.37500000e-01, 2.61708929e-01, 9.10321205e-02, + -1.54686328e-02, -2.74049397e-03, -7.94845816e-02, + 4.75368705e-01, 7.11647284e-02, 1.30266162e-01, + 3.37106977e-02, 1.06401886e-01, -7.31606787e-03, + -2.95625975e-03, -1.10250006e-02, 3.55194307e-01, + -1.44627826e-01, -2.89062500e-01, -9.28644588e-02, + -1.62557358e-01, 7.73431638e-02, -2.55329539e-03, + -1.90923851e-03, 1.57578403e-02, 1.72995854e-01, + -3.66267690e-01, -1.81657333e-01, -3.32521518e-01, + -2.59738162e-02, -2.31580576e-01, 4.20673902e-02, + -4.11710546e-04, -9.36449487e-04, 1.92156884e-02, + 2.82515641e-02, -3.90713738e-01, -1.69280296e-01, + -8.98437500e-02, -1.08693628e-01, 1.78813094e-01, + -1.98191857e-01, 1.65964201e-02, 2.77013853e-04]) + + test_vals = openmc.capi.math.calc_zn(n, rho, phi) + + assert np.allclose(ref_vals, test_vals) + + +def test_rotate_angle(): + uvw0 = np.array([1., 0., 0.]) + phi = 0. + mu = 0. + + # reference: mu of 0 pulls the vector the bottom, so: + ref_uvw = np.array([0., 0., -1.]) + + test_uvw = openmc.capi.math.rotate_angle(uvw0, mu, phi) + + assert np.array_equal(ref_uvw, test_uvw) + + # Repeat for mu = 1 (no change) + mu = 1. + ref_uvw = np.array([1., 0., 0.]) + + test_uvw = openmc.capi.math.rotate_angle(uvw0, mu, phi) + + assert np.array_equal(ref_uvw, test_uvw) + + # Now to test phi is None + mu = 0.9 + settings = openmc.capi.settings + settings.seed = 1 + + # When seed = 1, phi will be sampled as 1.9116495709698769 + # The resultant reference is from hand-calculations given the above + ref_uvw = [0.9, 0.410813051297112, 0.1457142302040] + test_uvw = openmc.capi.math.rotate_angle(uvw0, mu) + + assert np.allclose(ref_uvw, test_uvw) + + +def test_maxwell_spectrum(): + settings = openmc.capi.settings + settings.seed = 1 + T = 0.5 + ref_val = 0.6129982175261098 + test_val = openmc.capi.math.maxwell_spectrum(T) + + assert ref_val == test_val + + +def test_watt_spectrum(): + settings = openmc.capi.settings + settings.seed = 1 + a = 0.5 + b = 0.75 + ref_val = 0.6247242713640233 + test_val = openmc.capi.math.watt_spectrum(a, b) + + assert ref_val == test_val + + +def test_broaden_wmp_polynomials(): + # Two branches of the code to worry about, beta > 6 and otherwise + # beta = sqrtE * dopp + # First lets do beta > 6 + test_E = 0.5 + test_dopp = 100. # approximately U235 at room temperature + n = 6 + + ref_val = [2., 1.41421356, 1.0001, 0.70731891, 0.50030001, 0.353907] + test_val = openmc.capi.math.broaden_wmp_polynomials(test_E, test_dopp, n) + + assert np.allclose(ref_val, test_val) + + # now beta < 6 + test_dopp = 5. + ref_val = [1.99999885, 1.41421356, 1.04, 0.79195959, 0.6224, 0.50346003] + test_val = openmc.capi.math.broaden_wmp_polynomials(test_E, test_dopp, n) + + assert np.allclose(ref_val, test_val)