mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-26 05:05:30 -04:00
Merge pull request #1004 from nelsonag/cpp_math
Convert math.f90 to C++
This commit is contained in:
Paul Romano
GitHub
commit
04441ca619
15 changed files with 1463 additions and 871 deletions
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -47,3 +47,4 @@ from .mesh import *
|
|||
from .filter import *
|
||||
from .tally import *
|
||||
from .settings import settings
|
||||
from .math import *
|
||||
|
|
|
|||
250
openmc/capi/math.py
Normal file
250
openmc/capi/math.py
Normal file
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
947
src/math.F90
947
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
|
||||
|
|
|
|||
693
src/math_functions.cpp
Normal file
693
src/math_functions.cpp
Normal file
|
|
@ -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
|
||||
163
src/math_functions.h
Normal file
163
src/math_functions.h
Normal file
|
|
@ -0,0 +1,163 @@
|
|||
//! \file math_functions.h
|
||||
//! A collection of elementary math functions.
|
||||
|
||||
#ifndef MATH_FUNCTIONS_H
|
||||
#define MATH_FUNCTIONS_H
|
||||
|
||||
#include <cmath>
|
||||
#include <cstdlib>
|
||||
|
||||
#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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
212
tests/unit_tests/test_math.py
Normal file
212
tests/unit_tests/test_math.py
Normal file
|
|
@ -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)
|
||||
Loading…
Add table
Add a link
Reference in a new issue