mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
update for new wmp library v1.0 - python api and docs
This commit is contained in:
parent
fa9a6610af
commit
3b8d1ffb4b
4 changed files with 111 additions and 320 deletions
|
|
@ -5,7 +5,7 @@ Windowed Multipole Library Format
|
|||
=================================
|
||||
|
||||
**/version** (*char[]*)
|
||||
The format version of the file. The current version is "v0.2"
|
||||
The format version of the file. The current version is "v1.0"
|
||||
|
||||
**/nuclide/**
|
||||
- **broaden_poly** (*int[]*)
|
||||
|
|
@ -23,55 +23,25 @@ Windowed Multipole Library Format
|
|||
\text{data}[:,i] = [\text{pole},~\text{residue}_1,~\text{residue}_2,
|
||||
~\ldots]
|
||||
|
||||
The residues are in the order: total, competitive if present,
|
||||
absorption, fission. Complex numbers are stored by forming a type with
|
||||
":math:`r`" and ":math:`i`" identifiers, similar to how `h5py`_ does it.
|
||||
- **end_E** (*double*)
|
||||
The residues are in the order: scattering, absorption, fission. Complex
|
||||
numbers are stored by forming a type with ":math:`r`" and ":math:`i`"
|
||||
identifiers, similar to how `h5py`_ does it.
|
||||
- **E_max** (*double*)
|
||||
Highest energy the windowed multipole part of the library is valid for.
|
||||
- **formalism** (*int*)
|
||||
The formalism of the underlying data. Uses the `ENDF-6`_ format
|
||||
formalism numbers.
|
||||
|
||||
.. table:: Table of supported formalisms.
|
||||
|
||||
+-------------+------------------+
|
||||
| Formalism | Formalism number |
|
||||
+=============+==================+
|
||||
| MLBW | 2 |
|
||||
+-------------+------------------+
|
||||
| Reich-Moore | 3 |
|
||||
+-------------+------------------+
|
||||
|
||||
- **l_value** (*int[]*)
|
||||
The index for a corresponding pole. Equivalent to the :math:`l` quantum
|
||||
number of the resonance the pole comes from :math:`+1`.
|
||||
- **pseudo_K0RS** (*double[]*)
|
||||
:math:`l` dependent value of
|
||||
|
||||
.. math::
|
||||
\sqrt{\frac{2 m_n}{\hbar}}\frac{AWR}{AWR + 1} r_{s,l}
|
||||
|
||||
Where :math:`m_n` is mass of neutron, :math:`AWR` is the atomic weight
|
||||
ratio of the target to the neutron, and :math:`r_{s,l}` is the
|
||||
scattering radius for a given :math:`l`.
|
||||
- **E_min** (*double*)
|
||||
Lowest energy the windowed multipole part of the library is valid for.
|
||||
- **spacing** (*double*)
|
||||
.. math::
|
||||
\frac{\sqrt{E_{max}}- \sqrt{E_{min}}}{n_w}
|
||||
\frac{\sqrt{E_{max}} - \sqrt{E_{min}}}{n_w}
|
||||
|
||||
Where :math:`E_{max}` is the maximum energy the windows go up to. This
|
||||
is not equivalent to the maximum energy for which the windowed multipole
|
||||
data is valid for. It is slightly higher to ensure an integer number of
|
||||
windows. :math:`E_{min}` is the minimum energy and equivalent to
|
||||
``start_E``, and :math:`n_w` is the number of windows, given by
|
||||
``windows``.
|
||||
Where :math:`E_{max}` is the maximum energy the windows go up to.
|
||||
:math:`E_{min}` is the minimum energy, and :math:`n_w` is the number of
|
||||
windows, given by ``windows``.
|
||||
- **sqrtAWR** (*double*)
|
||||
Square root of the atomic weight ratio.
|
||||
- **start_E** (*double*)
|
||||
Lowest energy the windowed multipole part of the library is valid for.
|
||||
- **w_start** (*int[]*)
|
||||
The pole to start from for each window.
|
||||
- **w_end** (*int[]*)
|
||||
The pole to end at for each window.
|
||||
- **windows** (*int[][]*)
|
||||
The poles to start from and end at for each window. windows[i, 0] and
|
||||
windows[i, 1] are, respectively, the indexes (1-based) of the first and
|
||||
last pole in window i.
|
||||
|
||||
.. _h5py: http://docs.h5py.org/en/latest/
|
||||
.. _ENDF-6: https://www.oecd-nea.org/dbdata/data/manual-endf/endf102.pdf
|
||||
|
|
|
|||
|
|
@ -90,7 +90,7 @@ Assuming free-gas thermal motion, cross sections in the multipole form can be
|
|||
analytically Doppler broadened to give the form:
|
||||
|
||||
.. math::
|
||||
\sigma(E, T) = \frac{1}{2 E \sqrt{\xi}} \sum_j \text{Re} \left[i r_j
|
||||
\sigma(E, T) = \frac{1}{2 E \sqrt{\xi}} \sum_j \text{Re} \left[r_j
|
||||
\sqrt{\pi} W_i(z) - \frac{r_j}{\sqrt{\pi}} C \left(\frac{p_j}{\sqrt{\xi}},
|
||||
\frac{u}{2 \sqrt{\xi}}\right)\right]
|
||||
.. math::
|
||||
|
|
@ -141,7 +141,7 @@ scattering does not occur in the resolved resonance region. This is usually,
|
|||
but not always the case. Future library versions may eliminate this issue.
|
||||
|
||||
The data format used by OpenMC to represent windowed multipole data is specified
|
||||
in :ref:`io_data_wmp`.
|
||||
in :ref:`io_data_wmp` with a publicly available `WMP library`_.
|
||||
|
||||
.. _temperature_treatment:
|
||||
|
||||
|
|
@ -270,6 +270,7 @@ or even isotropic scattering.
|
|||
https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-14-24530.pdf
|
||||
.. _Hwang: http://www.ans.org/pubs/journals/nse/a_16381
|
||||
.. _Josey: http://dx.doi.org/10.1016/j.jcp.2015.08.013
|
||||
.. _WMP Library: https://github.com/mit-crpg/WMP_Library
|
||||
.. _MCNP: http://mcnp.lanl.gov
|
||||
.. _Serpent: http://montecarlo.vtt.fi
|
||||
.. _NJOY: http://t2.lanl.gov/codes.shtml
|
||||
|
|
|
|||
|
|
@ -4,7 +4,7 @@ HDF5_VERSION_MINOR = 0
|
|||
HDF5_VERSION = (HDF5_VERSION_MAJOR, HDF5_VERSION_MINOR)
|
||||
|
||||
# Version of WMP nuclear data format
|
||||
WMP_VERSION = 'v0.2'
|
||||
WMP_VERSION = 'v1.0'
|
||||
|
||||
|
||||
from .data import *
|
||||
|
|
|
|||
|
|
@ -10,26 +10,16 @@ import openmc.checkvalue as cv
|
|||
from openmc.mixin import EqualityMixin
|
||||
|
||||
|
||||
# Formalisms
|
||||
_FORM_MLBW = 2
|
||||
_FORM_RM = 3
|
||||
|
||||
# Constants that determine which value to access
|
||||
_MP_EA = 0 # Pole
|
||||
|
||||
# Reich-Moore indices
|
||||
_RM_RT = 1 # Residue total
|
||||
_RM_RA = 2 # Residue absorption
|
||||
_RM_RF = 3 # Residue fission
|
||||
|
||||
# Multi-level Breit Wigner indices
|
||||
_MLBW_RT = 1 # Residue total
|
||||
_MLBW_RX = 2 # Residue competitive
|
||||
_MLBW_RA = 3 # Residue absorption
|
||||
_MLBW_RF = 4 # Residue fission
|
||||
# Residue indices
|
||||
_MP_RS = 1 # Residue scattering
|
||||
_MP_RA = 2 # Residue absorption
|
||||
_MP_RF = 3 # Residue fission
|
||||
|
||||
# Polynomial fit indices
|
||||
_FIT_T = 0 # Total
|
||||
_FIT_S = 0 # Scattering
|
||||
_FIT_A = 1 # Absorption
|
||||
_FIT_F = 2 # Fission
|
||||
|
||||
|
|
@ -143,98 +133,61 @@ class WindowedMultipole(EqualityMixin):
|
|||
|
||||
Parameters
|
||||
----------
|
||||
formalism : {'MLBW', 'RM'}
|
||||
The R-matrix formalism used to reconstruct resonances. Either 'MLBW'
|
||||
for multi-level Breit Wigner or 'RM' for Reich-Moore.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
num_l : Integral
|
||||
Number of possible l quantum states for this nuclide.
|
||||
fit_order : Integral
|
||||
Order of the windowed curvefit.
|
||||
fissionable : bool
|
||||
Whether or not the target nuclide has fission data.
|
||||
formalism : {'MLBW', 'RM'}
|
||||
The R-matrix formalism used to reconstruct resonances. Either 'MLBW'
|
||||
for multi-level Breit Wigner or 'RM' for Reich-Moore.
|
||||
spacing : Real
|
||||
The width of each window in sqrt(E)-space. For example, the frst window
|
||||
will end at (sqrt(start_E) + spacing)**2 and the second window at
|
||||
(sqrt(start_E) + 2*spacing)**2.
|
||||
will end at (sqrt(E_min) + spacing)**2 and the second window at
|
||||
(sqrt(E_min) + 2*spacing)**2.
|
||||
sqrtAWR : Real
|
||||
Square root of the atomic weight ratio of the target nuclide.
|
||||
start_E : Real
|
||||
E_min : Real
|
||||
Lowest energy in eV the library is valid for.
|
||||
end_E : Real
|
||||
E_max : Real
|
||||
Highest energy in eV the library is valid for.
|
||||
data : np.ndarray
|
||||
A 2D array of complex poles and residues. data[i, 0] gives the energy
|
||||
at which pole i is located. data[i, 1:] gives the residues associated
|
||||
with the i-th pole. There are 3 residues for Reich-Moore data, one each
|
||||
for the total, absorption, and fission channels. Multi-level
|
||||
Breit Wigner data has an additional residue for the competitive channel.
|
||||
pseudo_k0RS : np.ndarray
|
||||
A 1D array of Real values. There is one value for each valid l
|
||||
quantum number. The values are equal to
|
||||
sqrt(2 m / hbar) * AWR / (AWR + 1) * r
|
||||
where m is the neutron mass, AWR is the atomic weight ratio, and r
|
||||
is the l-dependent scattering radius.
|
||||
l_value : np.ndarray
|
||||
A 1D array of Integral values equal to the l quantum number for each
|
||||
pole + 1.
|
||||
w_start : np.ndarray
|
||||
A 1D array of Integral values. w_start[i] - 1 is the index of the first
|
||||
pole in window i.
|
||||
w_end : np.ndarray
|
||||
A 1D array of Integral values. w_end[i] - 1 is the index of the last
|
||||
pole in window i.
|
||||
with the i-th pole. There are 3 residues, one each for the scattering,
|
||||
absorption, and fission channels.
|
||||
windows : np.ndarray
|
||||
A 2D array of Integral values. windows[i, 0] - 1 is the index of the
|
||||
first pole in window i. windows[i, 1] - 1 is the index of the last pole
|
||||
in window i.
|
||||
broaden_poly : np.ndarray
|
||||
A 1D array of boolean values indicating whether or not the polynomial
|
||||
curvefit in that window should be Doppler broadened.
|
||||
curvefit : np.ndarray
|
||||
A 3D array of Real curvefit polynomial coefficients. curvefit[i, 0, :]
|
||||
gives coefficients for the total cross section in window i.
|
||||
gives coefficients for the scattering cross section in window i.
|
||||
curvefit[i, 1, :] gives absorption coefficients and curvefit[i, 2, :]
|
||||
gives fission coefficients. The polynomial terms are increasing powers
|
||||
of sqrt(E) starting with 1/E e.g:
|
||||
a/E + b/sqrt(E) + c + d sqrt(E) + ...
|
||||
|
||||
"""
|
||||
def __init__(self, formalism):
|
||||
self._num_l = None
|
||||
self.formalism = formalism
|
||||
def __init__(self):
|
||||
self.spacing = None
|
||||
self.sqrtAWR = None
|
||||
self.start_E = None
|
||||
self.end_E = None
|
||||
self.E_min = None
|
||||
self.E_max = None
|
||||
self.data = None
|
||||
self.pseudo_k0RS = None
|
||||
self.l_value = None
|
||||
self.w_start = None
|
||||
self.w_end = None
|
||||
self.windows = None
|
||||
self.broaden_poly = None
|
||||
self.curvefit = None
|
||||
|
||||
@property
|
||||
def num_l(self):
|
||||
return self._num_l
|
||||
|
||||
@property
|
||||
def fit_order(self):
|
||||
return self.curvefit.shape[1] - 1
|
||||
|
||||
@property
|
||||
def fissionable(self):
|
||||
if self.formalism == 'RM':
|
||||
return self.data.shape[1] == 4
|
||||
else:
|
||||
# Assume self.formalism == 'MLBW'
|
||||
return self.data.shape[1] == 5
|
||||
|
||||
@property
|
||||
def formalism(self):
|
||||
return self._formalism
|
||||
return self.data.shape[1] == 4
|
||||
|
||||
@property
|
||||
def spacing(self):
|
||||
|
|
@ -245,32 +198,24 @@ class WindowedMultipole(EqualityMixin):
|
|||
return self._sqrtAWR
|
||||
|
||||
@property
|
||||
def start_E(self):
|
||||
return self._start_E
|
||||
def E_min(self):
|
||||
return self._E_min
|
||||
|
||||
@property
|
||||
def end_E(self):
|
||||
return self._end_E
|
||||
def E_max(self):
|
||||
return self._E_max
|
||||
|
||||
@property
|
||||
def data(self):
|
||||
return self._data
|
||||
|
||||
@property
|
||||
def pseudo_k0RS(self):
|
||||
return self._pseudo_k0RS
|
||||
|
||||
@property
|
||||
def l_value(self):
|
||||
return self._l_value
|
||||
|
||||
@property
|
||||
def w_start(self):
|
||||
return self._w_start
|
||||
|
||||
@property
|
||||
def w_end(self):
|
||||
return self._w_end
|
||||
def windows(self):
|
||||
return self._windows
|
||||
|
||||
@property
|
||||
def broaden_poly(self):
|
||||
|
|
@ -280,116 +225,64 @@ class WindowedMultipole(EqualityMixin):
|
|||
def curvefit(self):
|
||||
return self._curvefit
|
||||
|
||||
@formalism.setter
|
||||
def formalism(self, formalism):
|
||||
cv.check_type('formalism', formalism, str)
|
||||
cv.check_value('formalism', formalism, ('MLBW', 'RM'))
|
||||
self._formalism = formalism
|
||||
|
||||
@spacing.setter
|
||||
def spacing(self, spacing):
|
||||
if spacing is not None:
|
||||
cv.check_type('spacing', spacing, Real)
|
||||
cv.check_greater_than('spacing', spacing, 0.0, equality=False)
|
||||
check_type('spacing', spacing, Real)
|
||||
check_greater_than('spacing', spacing, 0.0, equality=False)
|
||||
self._spacing = spacing
|
||||
|
||||
@sqrtAWR.setter
|
||||
def sqrtAWR(self, sqrtAWR):
|
||||
if sqrtAWR is not None:
|
||||
cv.check_type('sqrtAWR', sqrtAWR, Real)
|
||||
cv.check_greater_than('sqrtAWR', sqrtAWR, 0.0, equality=False)
|
||||
check_type('sqrtAWR', sqrtAWR, Real)
|
||||
check_greater_than('sqrtAWR', sqrtAWR, 0.0, equality=False)
|
||||
self._sqrtAWR = sqrtAWR
|
||||
|
||||
@start_E.setter
|
||||
def start_E(self, start_E):
|
||||
if start_E is not None:
|
||||
cv.check_type('start_E', start_E, Real)
|
||||
cv.check_greater_than('start_E', start_E, 0.0, equality=True)
|
||||
self._start_E = start_E
|
||||
@E_min.setter
|
||||
def E_min(self, E_min):
|
||||
if E_min is not None:
|
||||
check_type('E_min', E_min, Real)
|
||||
check_greater_than('E_min', E_min, 0.0, equality=True)
|
||||
self._E_min = E_min
|
||||
|
||||
@end_E.setter
|
||||
def end_E(self, end_E):
|
||||
if end_E is not None:
|
||||
cv.check_type('end_E', end_E, Real)
|
||||
cv.check_greater_than('end_E', end_E, 0.0, equality=False)
|
||||
self._end_E = end_E
|
||||
@E_max.setter
|
||||
def E_max(self, E_max):
|
||||
if E_max is not None:
|
||||
check_type('E_max', E_max, Real)
|
||||
check_greater_than('E_max', E_max, 0.0, equality=False)
|
||||
self._E_max = E_max
|
||||
|
||||
@data.setter
|
||||
def data(self, data):
|
||||
if data is not None:
|
||||
cv.check_type('data', data, np.ndarray)
|
||||
check_type('data', data, np.ndarray)
|
||||
if len(data.shape) != 2:
|
||||
raise ValueError('Multipole data arrays must be 2D')
|
||||
if self.formalism == 'RM':
|
||||
if data.shape[1] not in (3, 4):
|
||||
raise ValueError('For the Reich-Moore formalism, '
|
||||
'data.shape[1] must be 3 or 4. One value for the pole.'
|
||||
' One each for the total and absorption residues. '
|
||||
'Possibly one more for a fission residue.')
|
||||
else:
|
||||
# Assume self.formalism == 'MLBW'
|
||||
if data.shape[1] not in (4, 5):
|
||||
raise ValueError('For the Multi-level Breit-Wigner '
|
||||
'formalism, data.shape[1] must be 4 or 5. One value '
|
||||
'for the pole. One each for the total, competitive, '
|
||||
'and absorption residues. Possibly one more for a '
|
||||
'fission residue.')
|
||||
if data.shape[1] not in (3, 4):
|
||||
raise ValueError(
|
||||
'data.shape[1] must be 3 or 4. One value for the pole.'
|
||||
' One each for the scattering and absorption residues. '
|
||||
'Possibly one more for a fission residue.')
|
||||
if not np.issubdtype(data.dtype, complex):
|
||||
raise TypeError('Multipole data arrays must be complex dtype')
|
||||
self._data = data
|
||||
|
||||
@pseudo_k0RS.setter
|
||||
def pseudo_k0RS(self, pseudo_k0RS):
|
||||
if pseudo_k0RS is not None:
|
||||
cv.check_type('pseudo_k0RS', pseudo_k0RS, np.ndarray)
|
||||
if len(pseudo_k0RS.shape) != 1:
|
||||
raise ValueError('Multipole pseudo_k0RS arrays must be 1D')
|
||||
if not np.issubdtype(pseudo_k0RS.dtype, float):
|
||||
raise TypeError('Multipole data arrays must be float dtype')
|
||||
self._pseudo_k0RS = pseudo_k0RS
|
||||
|
||||
@l_value.setter
|
||||
def l_value(self, l_value):
|
||||
if l_value is not None:
|
||||
cv.check_type('l_value', l_value, np.ndarray)
|
||||
if len(l_value.shape) != 1:
|
||||
raise ValueError('Multipole l_value arrays must be 1D')
|
||||
if not np.issubdtype(l_value.dtype, int):
|
||||
raise TypeError('Multipole l_value arrays must be integer'
|
||||
@windows.setter
|
||||
def windows(self, windows):
|
||||
if windows is not None:
|
||||
check_type('windows', windows, np.ndarray)
|
||||
if len(windows.shape) != 2:
|
||||
raise ValueError('Multipole windows arrays must be 2D')
|
||||
if not np.issubdtype(windows.dtype, int):
|
||||
raise TypeError('Multipole windows arrays must be integer'
|
||||
' dtype')
|
||||
|
||||
self._num_l = len(np.unique(l_value))
|
||||
|
||||
else:
|
||||
self._num_l = None
|
||||
|
||||
self._l_value = l_value
|
||||
|
||||
@w_start.setter
|
||||
def w_start(self, w_start):
|
||||
if w_start is not None:
|
||||
cv.check_type('w_start', w_start, np.ndarray)
|
||||
if len(w_start.shape) != 1:
|
||||
raise ValueError('Multipole w_start arrays must be 1D')
|
||||
if not np.issubdtype(w_start.dtype, int):
|
||||
raise TypeError('Multipole w_start arrays must be integer'
|
||||
' dtype')
|
||||
self._w_start = w_start
|
||||
|
||||
@w_end.setter
|
||||
def w_end(self, w_end):
|
||||
if w_end is not None:
|
||||
cv.check_type('w_end', w_end, np.ndarray)
|
||||
if len(w_end.shape) != 1:
|
||||
raise ValueError('Multipole w_end arrays must be 1D')
|
||||
if not np.issubdtype(w_end.dtype, int):
|
||||
raise TypeError('Multipole w_end arrays must be integer dtype')
|
||||
self._w_end = w_end
|
||||
self._windows = windows
|
||||
|
||||
@broaden_poly.setter
|
||||
def broaden_poly(self, broaden_poly):
|
||||
if broaden_poly is not None:
|
||||
cv.check_type('broaden_poly', broaden_poly, np.ndarray)
|
||||
check_type('broaden_poly', broaden_poly, np.ndarray)
|
||||
if len(broaden_poly.shape) != 1:
|
||||
raise ValueError('Multipole broaden_poly arrays must be 1D')
|
||||
if not np.issubdtype(broaden_poly.dtype, bool):
|
||||
|
|
@ -400,10 +293,10 @@ class WindowedMultipole(EqualityMixin):
|
|||
@curvefit.setter
|
||||
def curvefit(self, curvefit):
|
||||
if curvefit is not None:
|
||||
cv.check_type('curvefit', curvefit, np.ndarray)
|
||||
check_type('curvefit', curvefit, np.ndarray)
|
||||
if len(curvefit.shape) != 3:
|
||||
raise ValueError('Multipole curvefit arrays must be 3D')
|
||||
if curvefit.shape[2] not in (2, 3): # sig_t, sig_a (maybe sig_f)
|
||||
if curvefit.shape[2] not in (2, 3): # sig_s, sig_a (maybe sig_f)
|
||||
raise ValueError('The third dimension of multipole curvefit'
|
||||
' arrays must have a length of 2 or 3')
|
||||
if not np.issubdtype(curvefit.dtype, float):
|
||||
|
|
@ -443,19 +336,14 @@ class WindowedMultipole(EqualityMixin):
|
|||
'Python API expects version ' + WMP_VERSION)
|
||||
group = h5file['nuclide']
|
||||
|
||||
# Read scalars.
|
||||
out = cls()
|
||||
|
||||
if group['formalism'].value == _FORM_MLBW:
|
||||
out = cls('MLBW')
|
||||
elif group['formalism'].value == _FORM_RM:
|
||||
out = cls('RM')
|
||||
else:
|
||||
raise ValueError('Unrecognized/Unsupported R-matrix formalism')
|
||||
# Read scalars.
|
||||
|
||||
out.spacing = group['spacing'].value
|
||||
out.sqrtAWR = group['sqrtAWR'].value
|
||||
out.start_E = group['start_E'].value
|
||||
out.end_E = group['end_E'].value
|
||||
out.E_min = group['E_min'].value
|
||||
out.E_max = group['E_max'].value
|
||||
|
||||
# Read arrays.
|
||||
|
||||
|
|
@ -463,27 +351,15 @@ class WindowedMultipole(EqualityMixin):
|
|||
|
||||
out.data = group['data'].value
|
||||
|
||||
out.l_value = group['l_value'].value
|
||||
if out.l_value.shape[0] != out.data.shape[0]:
|
||||
raise ValueError(err.format('l_value', 'data'))
|
||||
|
||||
out.pseudo_k0RS = group['pseudo_K0RS'].value
|
||||
if out.pseudo_k0RS.shape[0] != out.num_l:
|
||||
raise ValueError(err.format('pseudo_k0RS', 'l_value'))
|
||||
|
||||
out.w_start = group['w_start'].value
|
||||
|
||||
out.w_end = group['w_end'].value
|
||||
if out.w_end.shape[0] != out.w_start.shape[0]:
|
||||
raise ValueError(err.format('w_end', 'w_start'))
|
||||
out.windows = group['windows'].value
|
||||
|
||||
out.broaden_poly = group['broaden_poly'].value.astype(np.bool)
|
||||
if out.broaden_poly.shape[0] != out.w_start.shape[0]:
|
||||
raise ValueError(err.format('broaden_poly', 'w_start'))
|
||||
if out.broaden_poly.shape[0] != out.windows.shape[0]:
|
||||
raise ValueError(err.format('broaden_poly', 'windows'))
|
||||
|
||||
out.curvefit = group['curvefit'].value
|
||||
if out.curvefit.shape[0] != out.w_start.shape[0]:
|
||||
raise ValueError(err.format('curvefit', 'w_start'))
|
||||
if out.curvefit.shape[0] != out.windows.shape[0]:
|
||||
raise ValueError(err.format('curvefit', 'windows'))
|
||||
|
||||
# _broaden_wmp_polynomials assumes the curve fit has at least 3 terms.
|
||||
if out.fit_order < 2:
|
||||
|
|
@ -493,7 +369,7 @@ class WindowedMultipole(EqualityMixin):
|
|||
return out
|
||||
|
||||
def _evaluate(self, E, T):
|
||||
"""Compute total, absorption, and fission cross sections.
|
||||
"""Compute scattering, absorption, and fission cross sections.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
|
|
@ -510,8 +386,8 @@ class WindowedMultipole(EqualityMixin):
|
|||
|
||||
"""
|
||||
|
||||
if E < self.start_E: return (0, 0, 0)
|
||||
if E > self.end_E: return (0, 0, 0)
|
||||
if E < self.E_min: return (0, 0, 0)
|
||||
if E > self.E_max: return (0, 0, 0)
|
||||
|
||||
# ======================================================================
|
||||
# Bookkeeping
|
||||
|
|
@ -525,34 +401,12 @@ class WindowedMultipole(EqualityMixin):
|
|||
# the 1-based vs. 0-based indexing. Similarly startw needs to be
|
||||
# decreased by 1. endw does not need to be decreased because
|
||||
# range(startw, endw) does not include endw.
|
||||
i_window = int(np.floor((sqrtE - sqrt(self.start_E)) / self.spacing))
|
||||
startw = self.w_start[i_window] - 1
|
||||
endw = self.w_end[i_window]
|
||||
|
||||
# Fill in factors. Because of the unique interference dips in scatering
|
||||
# resonances, the total cross section has a special "factor" that does
|
||||
# not appear in the absorption and fission equations.
|
||||
if startw <= endw:
|
||||
twophi = np.zeros(self.num_l, dtype=np.float)
|
||||
sig_t_factor = np.zeros(self.num_l, dtype=np.cfloat)
|
||||
|
||||
for iL in range(self.num_l):
|
||||
twophi[iL] = self.pseudo_k0RS[iL] * sqrtE
|
||||
if iL == 1:
|
||||
twophi[iL] = twophi[iL] - np.arctan(twophi[iL])
|
||||
elif iL == 2:
|
||||
arg = 3.0 * twophi[iL] / (3.0 - twophi[iL]**2)
|
||||
twophi[iL] = twophi[iL] - np.arctan(arg)
|
||||
elif iL == 3:
|
||||
arg = (twophi[iL] * (15.0 - twophi[iL]**2)
|
||||
/ (15.0 - 6.0 * twophi[iL]**2))
|
||||
twophi[iL] = twophi[iL] - np.arctan(arg)
|
||||
|
||||
twophi = 2.0 * twophi
|
||||
sig_t_factor = np.cos(twophi) - 1j*np.sin(twophi)
|
||||
i_window = int(np.floor((sqrtE - sqrt(self.E_min)) / self.spacing))
|
||||
startw = self.windows[i_window, 0] - 1
|
||||
endw = self.windows[i_window, 1]
|
||||
|
||||
# Initialize the ouptut cross sections.
|
||||
sig_t = 0.0
|
||||
sig_s = 0.0
|
||||
sig_a = 0.0
|
||||
sig_f = 0.0
|
||||
|
||||
|
|
@ -565,7 +419,7 @@ class WindowedMultipole(EqualityMixin):
|
|||
broadened_polynomials = _broaden_wmp_polynomials(E, dopp,
|
||||
self.fit_order + 1)
|
||||
for i_poly in range(self.fit_order+1):
|
||||
sig_t += (self.curvefit[i_window, i_poly, _FIT_T]
|
||||
sig_s += (self.curvefit[i_window, i_poly, _FIT_S]
|
||||
* broadened_polynomials[i_poly])
|
||||
sig_a += (self.curvefit[i_window, i_poly, _FIT_A]
|
||||
* broadened_polynomials[i_poly])
|
||||
|
|
@ -575,7 +429,7 @@ class WindowedMultipole(EqualityMixin):
|
|||
else:
|
||||
temp = invE
|
||||
for i_poly in range(self.fit_order+1):
|
||||
sig_t += self.curvefit[i_window, i_poly, _FIT_T] * temp
|
||||
sig_s += self.curvefit[i_window, i_poly, _FIT_S] * temp
|
||||
sig_a += self.curvefit[i_window, i_poly, _FIT_A] * temp
|
||||
if self.fissionable:
|
||||
sig_f += self.curvefit[i_window, i_poly, _FIT_F] * temp
|
||||
|
|
@ -589,22 +443,10 @@ class WindowedMultipole(EqualityMixin):
|
|||
for i_pole in range(startw, endw):
|
||||
psi_chi = -1j / (self.data[i_pole, _MP_EA] - sqrtE)
|
||||
c_temp = psi_chi / E
|
||||
if self.formalism == 'MLBW':
|
||||
sig_t += ((self.data[i_pole, _MLBW_RT] * c_temp *
|
||||
sig_t_factor[self.l_value[i_pole]-1]).real
|
||||
+ (self.data[i_pole, _MLBW_RX] * c_temp).real)
|
||||
sig_a += (self.data[i_pole, _MLBW_RA] * c_temp).real
|
||||
if self.fissionable:
|
||||
sig_f += (self.data[i_pole, _MLBW_RF] * c_temp).real
|
||||
elif self.formalism == 'RM':
|
||||
sig_t += (self.data[i_pole, _RM_RT] * c_temp *
|
||||
sig_t_factor[self.l_value[i_pole]-1]).real
|
||||
sig_a += (self.data[i_pole, _RM_RA] * c_temp).real
|
||||
if self.fissionable:
|
||||
sig_f += (self.data[i_pole, _RM_RF] * c_temp).real
|
||||
else:
|
||||
raise ValueError('Unrecognized/Unsupported R-matrix'
|
||||
' formalism')
|
||||
sig_s += (self.data[i_pole, _MP_RS] * c_temp).real
|
||||
sig_a += (self.data[i_pole, _MP_RA] * c_temp).real
|
||||
if self.fissionable:
|
||||
sig_f += (self.data[i_pole, _MP_RF] * c_temp).real
|
||||
|
||||
else:
|
||||
# At temperature, use Faddeeva function-based form.
|
||||
|
|
@ -612,27 +454,15 @@ class WindowedMultipole(EqualityMixin):
|
|||
for i_pole in range(startw, endw):
|
||||
Z = (sqrtE - self.data[i_pole, _MP_EA]) * dopp
|
||||
w_val = _faddeeva(Z) * dopp * invE * sqrt(pi)
|
||||
if self.formalism == 'MLBW':
|
||||
sig_t += ((self.data[i_pole, _MLBW_RT] *
|
||||
sig_t_factor[self.l_value[i_pole]-1] +
|
||||
self.data[i_pole, _MLBW_RX]) * w_val).real
|
||||
sig_a += (self.data[i_pole, _MLBW_RA] * w_val).real
|
||||
if self.fissionable:
|
||||
sig_f += (self.data[i_pole, _MLBW_RF] * w_val).real
|
||||
elif self.formalism == 'RM':
|
||||
sig_t += (self.data[i_pole, _RM_RT] * w_val *
|
||||
sig_t_factor[self.l_value[i_pole]-1]).real
|
||||
sig_a += (self.data[i_pole, _RM_RA] * w_val).real
|
||||
if self.fissionable:
|
||||
sig_f += (self.data[i_pole, _RM_RF] * w_val).real
|
||||
else:
|
||||
raise ValueError('Unrecognized/Unsupported R-matrix'
|
||||
' formalism')
|
||||
sig_s += (self.data[i_pole, _MP_RS] * w_val).real
|
||||
sig_a += (self.data[i_pole, _MP_RA] * w_val).real
|
||||
if self.fissionable:
|
||||
sig_f += (self.data[i_pole, _MP_RF] * w_val).real
|
||||
|
||||
return sig_t, sig_a, sig_f
|
||||
return sig_s, sig_a, sig_f
|
||||
|
||||
def __call__(self, E, T):
|
||||
"""Compute total, absorption, and fission cross sections.
|
||||
"""Compute scattering, absorption, and fission cross sections.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
|
|
@ -674,24 +504,14 @@ class WindowedMultipole(EqualityMixin):
|
|||
g = f.create_group('nuclide')
|
||||
|
||||
# Write scalars.
|
||||
if self.formalism == 'MLBW':
|
||||
g.create_dataset('formalism',
|
||||
data=np.array(_FORM_MLBW, dtype=np.int32))
|
||||
else:
|
||||
# Assume RM.
|
||||
g.create_dataset('formalism',
|
||||
data=np.array(_FORM_RM, dtype=np.int32))
|
||||
g.create_dataset('spacing', data=np.array(self.spacing))
|
||||
g.create_dataset('sqrtAWR', data=np.array(self.sqrtAWR))
|
||||
g.create_dataset('start_E', data=np.array(self.start_E))
|
||||
g.create_dataset('end_E', data=np.array(self.end_E))
|
||||
g.create_dataset('E_min', data=np.array(self.E_min))
|
||||
g.create_dataset('E_max', data=np.array(self.E_max))
|
||||
|
||||
# Write arrays.
|
||||
g.create_dataset('data', data=self.data)
|
||||
g.create_dataset('l_value', data=self.l_value)
|
||||
g.create_dataset('pseudo_K0RS', data=self.pseudo_k0RS)
|
||||
g.create_dataset('w_start', data=self.w_start)
|
||||
g.create_dataset('w_end', data=self.w_end)
|
||||
g.create_dataset('windows', data=self.windows)
|
||||
g.create_dataset('broaden_poly',
|
||||
data=self.broaden_poly.astype(np.int8))
|
||||
g.create_dataset('curvefit', data=self.curvefit)
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue