crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 21:55:41 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

CPPized the Pn and Rn functions

Browse source
This commit is contained in:
Adam G Nelson 2018-04-30 16:13:11 -04:00
parent 0e3194315a
commit 15b6b1d4ee
5 changed files with 84 additions and 495 deletions
Download patch file Download diff file Expand all files Collapse all files

48
openmc/capi/math.py
View file

@ -98,6 +98,30 @@ def calc_pn(n, x):
return _dll.calc_pn(c_int(n), c_double(x))
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_int(len(data)),
data_arr.ctypes.data_as(POINTER(c_double)),
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 (where -n <= m <= n).
@ -151,30 +175,6 @@ def calc_zn(n, rho, phi):
return zn
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_int(len(data)),
data_arr.ctypes.data_as(POINTER(c_double)),
c_double(x))
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.

450
src/math.F90
View file

@ -148,7 +148,7 @@ contains
! the return value will be 1.0.
!===============================================================================
pure function calc_pn(n,x) result(pnx) bind(C)
pure function calc_pn(n, x) result(pnx) bind(C)
integer(C_INT), intent(in) :: n ! Legendre order requested
real(C_DOUBLE), intent(in) :: x ! Independent variable the Legendre is to
@ -157,39 +157,25 @@ contains
real(C_DOUBLE) :: 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
pnx = calc_pn_cc(n, x)
end function calc_pn
!===============================================================================
! EVALUATE_LEGENDRE Find the value of f(x) given a set of Legendre coefficients
! and the value of x
!===============================================================================
pure function evaluate_legendre(n, data, x) result(val) bind(C)
integer(C_INT), intent(in) :: n
real(C_DOUBLE), intent(in) :: data(n)
real(C_DOUBLE), intent(in) :: x
real(C_DOUBLE) :: val
val = evaluate_legendre_cc(size(data), data, x)
end function evaluate_legendre
!===============================================================================
! 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)
@ -202,387 +188,7 @@ contains
! assumed to be on unit sphere
real(C_DOUBLE) :: rn(2*n + 1) ! The resultant R_n(uvw)
real(C_DOUBLE) :: phi, w ! Azimuthal and Cosine of Polar angles (from uvw)
real(C_DOUBLE) :: 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
call calc_rn_cc(n, uvw, rn)
end subroutine calc_rn
@ -693,26 +299,6 @@ contains
end do
end subroutine calc_zn
!===============================================================================
! EVALUATE_LEGENDRE Find the value of f(x) given a set of Legendre coefficients
! and the value of x
!===============================================================================
pure function evaluate_legendre(n, data, x) result(val) bind(C)
integer(C_INT), intent(in) :: n
real(C_DOUBLE), intent(in) :: data(n)
real(C_DOUBLE), intent(in) :: x
real(C_DOUBLE) :: val
integer(C_INT) :: l
val = HALF * data(1)
do l = 1, n - 1
val = val + (real(l, 8) + HALF) * data(l + 1) * calc_pn(l,x)
end do
end function evaluate_legendre
!===============================================================================
! ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is
! mu and through an azimuthal angle sampled uniformly. Note that this is done

31
src/math_functions.cpp
View file

@ -161,6 +161,22 @@ double __attribute__ ((const)) calc_pn_c(int n, double x) {
return pnx;
}
//==============================================================================
// EVALUATE_LEGENDRE Find the value of f(x) given a set of Legendre coefficients
// and the value of x
//==============================================================================
double __attribute__ ((const)) evaluate_legendre_c(int n, double data[],
double x) {
double val;
val = 0.5 * data[0];
for (int l = 1; l < n; l++) {
val += (static_cast<double>(l) + 0.5) * data[l] * calc_pn_c(l, x);
}
return val;
}
//==============================================================================
// CALC_RN calculates the n-th order spherical harmonics for a given angle
// (in terms of (u,v,w)). All Rn,m values are provided (where -n<=m<=n)
@ -561,21 +577,6 @@ void calc_rn_c(int n, double uvw[3], double rn[]){
}
}
//==============================================================================
// EVALUATE_LEGENDRE Find the value of f(x) given a set of Legendre coefficients
// and the value of x
//==============================================================================
double __attribute__ ((const)) evaluate_legendre_c(int n, double data[],
double x) {
double val;
val = 0.5 * data[0];
for (int l = 1; l < n; l++) {
val += (static_cast<double>(l) + 0.5) * data[l] * calc_pn_c(l, x);
}
}
//==============================================================================
// ROTATE_ANGLE rotates direction std::cosines through a polar angle whose

14
src/math_functions.h
View file

@ -32,13 +32,6 @@ extern "C" double t_percentile_c(double p, int df) __attribute__ ((const));
extern "C" double calc_pn_c(int n, double x) __attribute__ ((const));
//==============================================================================
// CALC_RN calculates the n-th order spherical harmonics for a given angle
// (in terms of (u,v,w)). All Rn,m values are provided (where -n<=m<=n)
//==============================================================================
extern "C" void calc_rn_c(int n, double uvw[3], double rn[]);
//==============================================================================
// EVALUATE_LEGENDRE Find the value of f(x) given a set of Legendre coefficients
// and the value of x
@ -47,6 +40,13 @@ extern "C" void calc_rn_c(int n, double uvw[3], double rn[]);
extern "C" double evaluate_legendre_c(int n, double data[], double x)
__attribute__ ((const));
//==============================================================================
// CALC_RN calculates the n-th order spherical harmonics for a given angle
// (in terms of (u,v,w)). All Rn,m values are provided (where -n<=m<=n)
//==============================================================================
extern "C" void calc_rn_c(int n, double uvw[3], double rn[]);
//==============================================================================
// ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is
// mu and through an azimuthal angle sampled uniformly. Note that this is done

36
tests/unit_tests/test_math.py
View file

@ -51,6 +51,25 @@ def test_calc_pn():
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 includes the (2l+1)/2 term
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)])
@ -99,23 +118,6 @@ def test_calc_zn():
pass
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
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_rotate_angle():
uvw0 = np.array([1., 0., 0.])
phi = 0.