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

Merge pull request #1003 from lindsayad/consistent-norms

Browse source 
Use consistent normalization for Zernikes
This commit is contained in:
Paul Romano 2018-05-09 06:34:05 -05:00 • committed by GitHub
commit 6fa9b02b4d
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
2 changed files with 10 additions and 15 deletions
Download patch file Download diff file Expand all files Collapse all files

13
openmc/filter_expansion.py
View file

@ -325,15 +325,15 @@ class ZernikeFilter(ExpansionFilter):
This filter allows scores to be multiplied by Zernike polynomials of the
particle's position normalized to a given unit circle, up to a
user-specified order. The Zernike polynomials follow the definition by `Noll
<https://doi.org/10.1364/JOSA.66.000207>`_ and are defined as
user-specified order. The standard Zernike polynomials follow the definition by
Born and Wolf, *Principles of Optics* and are defined as
.. math::
Z_n^m(\rho, \theta) = \sqrt{2n + 2} R_n^m(\rho) \cos (m\theta), \quad m > 0
Z_n^m(\rho, \theta) = R_n^m(\rho) \cos (m\theta), \quad m > 0
Z_n^{m}(\rho, \theta) = \sqrt{2n + 2} R_n^{m}(\rho) \sin (m\theta), \quad m < 0
Z_n^{m}(\rho, \theta) = R_n^{m}(\rho) \sin (m\theta), \quad m < 0
Z_n^{m}(\rho, \theta) = \sqrt{n + 1} R_n^{m}(\rho), \quad m = 0
Z_n^{m}(\rho, \theta) = R_n^{m}(\rho), \quad m = 0
where the radial polynomials are
@ -342,7 +342,8 @@ class ZernikeFilter(ExpansionFilter):
\frac{n+m}{2} - k)! (\frac{n-m}{2} - k)!} \rho^{n-2k}.
With this definition, the integral of :math:`(Z_n^m)^2` over the unit disk
is exactly :math:`\pi` for each polynomial.
is :math:`\frac{\epsilon_m\pi}{2n+2}` for each polynomial where :math:`\epsilon_m` is
2 if :math:`m` equals 0 and 1 otherwise.
Specifying a filter with order N tallies moments for all :math:`n` from 0 to
N and each value of :math:`m`. The ordering of the Zernike polynomial

12
src/math.F90
View file

@ -600,12 +600,6 @@ contains
real(8) :: k1, k2, k3, k4 ! Variables for R_m_n calculation
integer :: i,p,q ! Loop counters
real(8), parameter :: SQRT_N_1(0:10) = [&
sqrt(1.0_8), sqrt(2.0_8), sqrt(3.0_8), sqrt(4.0_8), &
sqrt(5.0_8), sqrt(6.0_8), sqrt(7.0_8), sqrt(8.0_8), &
sqrt(9.0_8), sqrt(10.0_8), sqrt(11.0_8)]
real(8), parameter :: SQRT_2N_2(0:10) = SQRT_N_1*sqrt(2.0_8)
! n == radial degree
! m == azimuthal frequency
@ -669,11 +663,11 @@ contains
do p = 0, n
do q = -p, p, 2
if (q < 0) then
zn(i) = zn_mat(p+1, abs(q)+1) * sin_phi_vec(abs(q)) * SQRT_2N_2(p)
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) * SQRT_N_1(p)
zn(i) = zn_mat(p+1, q+1)
else
zn(i) = zn_mat(p+1, q+1) * cos_phi_vec(abs(q)+1) * SQRT_2N_2(p)
zn(i) = zn_mat(p+1, q+1) * cos_phi_vec(abs(q)+1)
end if
i = i + 1
end do