Merge pull request #1054 from liangjg/wmp_v1.0

Accommodate the new WMP Library

.travis.yml

@@ -17,6 +17,7 @@ cache:
   directories:
   - $HOME/nndc_hdf5
   - $HOME/endf-b-vii.1
+  - $HOME/WMP_Library
 env:
   global:
     - FC=gfortran
@@ -25,7 +26,7 @@ env:
     - OMP_NUM_THREADS=2
     - OPENMC_CROSS_SECTIONS=$HOME/nndc_hdf5/cross_sections.xml
     - OPENMC_ENDF_DATA=$HOME/endf-b-vii.1
-    - OPENMC_MULTIPOLE_LIBRARY=$HOME/multipole_lib
+    - OPENMC_MULTIPOLE_LIBRARY=$HOME/WMP_Library
     - PATH=$PATH:$HOME/NJOY2016/build
     - DISPLAY=:99.0
     - COVERALLS_PARALLEL=true

docs/source/io_formats/data_wmp.rst

@@ -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

docs/source/methods/cross_sections.rst

@@ -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

examples/jupyter/nuclear-data.ipynb

@@ -1466,7 +1466,7 @@
    },
    "outputs": [],
    "source": [
-    "url = 'https://anl.box.com/shared/static/ulhcoohm12gduwdalknmf8dpnepzkxj0.h5'\n",
+    "url = 'https://github.com/mit-crpg/WMP_Library/releases/download/v1.0/092238.h5'\n",
     "filename, headers = urllib.request.urlretrieve(url, '092238.h5')"
    ]
   },
@@ -1485,7 +1485,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The `WindowedMultipole` object can be called with energy and temperature values.  Calling the object gives a tuple of 3 cross sections: total, radiative capture, and fission."
+    "The `WindowedMultipole` object can be called with energy and temperature values.  Calling the object gives a tuple of 3 cross sections: elastic scattering, radiative capture, and fission."
    ]
   },
   {
@@ -1498,9 +1498,7 @@
     {
      "data": {
       "text/plain": [
-       "(array(9.638243132516015),\n",
-       " array(0.5053244245010787),\n",
-       " array(2.931753364280356e-06))"
+       "(array(9.13284265), array(0.50530278), array(2.9316765e-06))"
       ]
      },
      "execution_count": 43,

openmc/data/__init__.py

@@ -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 *

openmc/data/multipole.py

@@ -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,20 @@ 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,12 +221,6 @@ 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:
@@ -300,19 +235,19 @@ class WindowedMultipole(EqualityMixin):
             cv.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:
+            cv.check_type('E_min', E_min, Real)
+            cv.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:
+            cv.check_type('E_max', E_max, Real)
+            cv.check_greater_than('E_max', E_max, 0.0, equality=False)
+        self._E_max = E_max
 
     @data.setter
     def data(self, data):
@@ -320,71 +255,25 @@ class WindowedMultipole(EqualityMixin):
             cv.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 not np.issubdtype(data.dtype, complex):
+            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, np.complexfloating):
                 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:
+            cv.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, np.integer):
+                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):
@@ -392,7 +281,7 @@ class WindowedMultipole(EqualityMixin):
             cv.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):
+            if not np.issubdtype(broaden_poly.dtype, np.bool_):
                 raise TypeError('Multipole broaden_poly arrays must be boolean'
                                 ' dtype')
         self._broaden_poly = broaden_poly
@@ -403,10 +292,10 @@ class WindowedMultipole(EqualityMixin):
             cv.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):
+            if not np.issubdtype(curvefit.dtype, np.floating):
                 raise TypeError('Multipole curvefit arrays must be float dtype')
         self._curvefit = curvefit
 
@@ -443,19 +332,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 +347,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 +365,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 +382,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 +397,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 +415,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 +425,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 +439,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 +450,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 +500,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)

scripts/openmc-get-multipole-data

@@ -29,9 +29,9 @@ parser.add_argument('-b', '--batch', action='store_true',
 args = parser.parse_args()
 
 
-baseUrl = 'https://github.com/smharper/windowed_multipole_library/blob/master/'
-files = ['multipole_lib.tar.gz?raw=true']
-checksums = ['3985aea96f7162a9419c7ed8352e6abb']
+baseUrl = 'https://github.com/mit-crpg/WMP_Library/releases/download/v1.0/'
+files = ['WMP_Library_v1.0.tar.gz']
+checksums = ['22cb675734cfccb278dffd40dcfbf26a']
 block_size = 16384
 
 # ==============================================================================
@@ -101,12 +101,7 @@ for f in files:
     # Extract files
     with tarfile.open(fname, 'r') as tgz:
         print('Extracting {0}...'.format(fname))
-        tgz.extractall(path='wmp/')
-
-# Move data files down one level
-for filename in glob.glob('wmp/multipole_lib/*'):
-    shutil.move(filename, 'wmp/')
-os.rmdir('wmp/multipole_lib')
+        tgz.extractall(path='')
 
 # ==============================================================================
 # PROMPT USER TO DELETE .TAR.GZ FILES

src/constants.F90

@@ -25,7 +25,7 @@ module constants
   integer, parameter :: VERSION_VOLUME(2)           = [1, 0]
   integer, parameter :: VERSION_VOXEL(2)            = [1, 0]
   integer, parameter :: VERSION_MGXS_LIBRARY(2)     = [1, 0]
-  character(10), parameter :: VERSION_MULTIPOLE     = "v0.2"
+  character(10), parameter :: VERSION_MULTIPOLE     = "v1.0"
 
   ! ============================================================================
   ! ADJUSTABLE PARAMETERS

src/multipole_header.F90

@@ -10,27 +10,16 @@ module multipole_header
   !========================================================================
   ! Multipole related constants
 
-  ! Formalisms
-  integer, parameter :: FORM_MLBW = 2, &
-                        FORM_RM   = 3, &
-                        FORM_RML  = 7
-
   ! Constants that determine which value to access
   integer, parameter :: MP_EA = 1       ! Pole
 
-  ! Reich-Moore indices
-  integer, parameter :: RM_RT = 2, &    ! Residue total
-                        RM_RA = 3, &    ! Residue absorption
-                        RM_RF = 4       ! Residue fission
-
-  ! Multi-level Breit Wigner indices
-  integer, parameter :: MLBW_RT = 2, &  ! Residue total
-                        MLBW_RX = 3, &  ! Residue compettitive
-                        MLBW_RA = 4, &  ! Residue absorption
-                        MLBW_RF = 5     ! Residue fission
+  ! Residue indices
+  integer, parameter :: MP_RS = 2, &    ! Residue scattering
+                        MP_RA = 3, &    ! Residue absorption
+                        MP_RF = 4       ! Residue fission
 
   ! Polynomial fit indices
-  integer, parameter :: FIT_T = 1, &    ! Total
+  integer, parameter :: FIT_S = 1, &    ! Scattering
                         FIT_A = 2, &    ! Absorption
                         FIT_F = 3       ! Fission
 
@@ -46,33 +35,22 @@ module multipole_header
     ! Isotope Properties
 
     logical                 :: fissionable     ! Is this isotope fissionable?
-    integer, allocatable    :: l_value(:)      ! The l index of the pole
-    integer                 :: num_l           ! Number of unique l values
-    real(8), allocatable    :: pseudo_k0RS(:)  ! The value (sqrt(2*mass neutron
-                                               !   /reduced planck constant)
-                                               !   * AWR/(AWR + 1)
-                                               !   * scattering radius for
-                                               !    each l
     complex(8), allocatable :: data(:,:)       ! Poles and residues
     real(8)                 :: sqrtAWR         ! Square root of the atomic
                                                ! weight ratio
-    integer                 :: formalism       ! R-matrix formalism
 
     !=========================================================================
     ! Windows
 
     integer :: fit_order                    ! Order of the fit. 1 linear,
                                             !   2 quadratic, etc.
-    real(8) :: start_E                      ! Start energy for the windows
-    real(8) :: end_E                        ! End energy for the windows
+    real(8) :: E_min                        ! Start energy for the windows
+    real(8) :: E_max                        ! End energy for the windows
     real(8) :: spacing                      ! The actual spacing in sqrt(E)
                                             !   space.
-    ! spacing = sqrt(multipole_w % endE - multipole_w % startE)
-    !           / multipole_w % windows
-    integer, allocatable :: w_start(:)      ! Contains the index of the pole at
-                                            !   the start of the window
-    integer, allocatable :: w_end(:)        ! Contains the index of the pole at
-                                            !   the end of the window
+    integer, allocatable :: windows(:, :)   ! Contains the indexes of the poles
+                                            !   at the start and end of the
+                                            !   window
     real(8), allocatable :: curvefit(:,:,:) ! Contains the fitting function.
                                             !   (reaction type, coeff index,
                                             !   window index)
@@ -96,7 +74,7 @@ contains
     character(len=*),      intent(in)    :: filename
 
     character(len=10) :: version
-    integer :: i, n_poles, n_residue_types, n_windows
+    integer :: i, n_poles, n_residues, n_windows
     integer(HSIZE_T) :: dims_1d(1), dims_2d(2), dims_3d(3)
     integer(HID_T) :: file_id
     integer(HID_T) :: group_id
@@ -114,65 +92,49 @@ contains
          // trim(filename) // " uses version " // trim(version) // ".")
 
     ! Read scalar values.
-    call read_dataset(this % formalism, group_id, "formalism")
     call read_dataset(this % spacing, group_id, "spacing")
     call read_dataset(this % sqrtAWR, group_id, "sqrtAWR")
-    call read_dataset(this % start_E, group_id, "start_E")
-    call read_dataset(this % end_E, group_id, "end_E")
+    call read_dataset(this % E_min, group_id, "E_min")
+    call read_dataset(this % E_max, group_id, "E_max")
 
     ! Read the "data" array.  Use its shape to figure out the number of poles
     ! and residue types in this data.
     dset = open_dataset(group_id, "data")
     call get_shape(dset, dims_2d)
-    n_residue_types = int(dims_2d(1), 4) - 1
+    n_residues = int(dims_2d(1), 4) - 1
     n_poles = int(dims_2d(2), 4)
-    allocate(this % data(n_residue_types+1, n_poles))
-    call read_dataset(this % data, dset)
+    allocate(this % data(n_residues+1, n_poles))
+    if (n_poles > 0) call read_dataset(this % data, dset)
     call close_dataset(dset)
 
     ! Check to see if this data includes fission residues.
-    if (this % formalism == FORM_RM) then
-      this % fissionable = (n_residue_types == 3)
-    else
-      ! Assume FORM_MLBW.
-      this % fissionable = (n_residue_types == 4)
-    end if
+    this % fissionable = (n_residues == 3)
 
-    ! Read the "l_value" array.
-    allocate(this % l_value(n_poles))
-    call read_dataset(this % l_value, group_id, "l_value")
-
-    ! Figure out the number of unique l values in the l_value array.
-    do i = 1, n_poles
-      if (.not. l_val_dict % has(this % l_value(i))) then
-        call l_val_dict % set(this % l_value(i), 0)
-      end if
-    end do
-    this % num_l = l_val_dict % size()
-    call l_val_dict % clear()
-
-    ! Read the "pseudo_K0RS" array.
-    allocate(this % pseudo_k0RS(this % num_l))
-    call read_dataset(this % pseudo_k0RS, group_id, "pseudo_K0RS")
-
-    ! Read the "w_start" array and use its shape to figure out the number of
+    ! Read the "windows" array and use its shape to figure out the number of
     ! windows.
-    dset = open_dataset(group_id, "w_start")
-    call get_shape(dset, dims_1d)
-    n_windows = int(dims_1d(1), 4)
-    allocate(this % w_start(n_windows))
-    call read_dataset(this % w_start, dset)
+    dset = open_dataset(group_id, "windows")
+    call get_shape(dset, dims_2d)
+    n_windows = int(dims_2d(2), 4)
+    allocate(this % windows(dims_2d(1), n_windows))
+    call read_dataset(this % windows, dset)
     call close_dataset(dset)
 
-    ! Read the "w_end" and "broaden_poly" arrays.
-    allocate(this % w_end(n_windows))
-    call read_dataset(this % w_end, group_id, "w_end")
+    ! Read the "broaden_poly" arrays.
+    dset = open_dataset(group_id, "broaden_poly")
+    call get_shape(dset, dims_1d)
+    if (dims_1d(1) /= n_windows) call fatal_error("broaden_poly array shape is&
+         &not consistent with the windows array shape in multipole library"&
+         // trim(filename) // ".")
     allocate(this % broaden_poly(n_windows))
-    call read_dataset(this % broaden_poly, group_id, "broaden_poly")
+    call read_dataset(this % broaden_poly, dset)
+    call close_dataset(dset)
 
     ! Read the "curvefit" array.
     dset = open_dataset(group_id, "curvefit")
     call get_shape(dset, dims_3d)
+    if (dims_3d(3) /= n_windows) call fatal_error("curvefit array shape is not&
+         &consistent with the windows array shape in multipole library"&
+         // trim(filename) // ".")
     allocate(this % curvefit(dims_3d(1), dims_3d(2), dims_3d(3)))
     call read_dataset(this % curvefit, dset)
     call close_dataset(dset)

src/nuclide_header.F90

@@ -12,11 +12,9 @@ module nuclide_header
   use hdf5_interface
   use math,                   only: faddeeva, w_derivative, &
                                     broaden_wmp_polynomials
-  use multipole_header,       only: FORM_RM, FORM_MLBW, MP_EA, RM_RT, RM_RA, &
-                                    RM_RF, MLBW_RT, MLBW_RX, MLBW_RA, MLBW_RF, &
-                                    FIT_T, FIT_A, FIT_F, MultipoleArray
+  use multipole_header,       only: MP_EA, MP_RS, MP_RA, MP_RF, &
+                                    FIT_S, FIT_A, FIT_F, MultipoleArray
   use message_passing
-  use multipole_header,       only: MultipoleArray
   use random_lcg,             only: prn, future_prn, prn_set_stream
   use reaction_header,        only: Reaction
   use sab_header,             only: SAlphaBeta, sab_tables
@@ -849,7 +847,7 @@ contains
     integer :: threshold ! threshold energy index
     real(8) :: f      ! interp factor on nuclide energy grid
     real(8) :: kT     ! temperature in eV
-    real(8) :: sig_t, sig_a, sig_f ! Intermediate multipole variables
+    real(8) :: sig_s, sig_a, sig_f ! Intermediate multipole variables
 
     ! Initialize cached cross sections to zero
     micro_xs % elastic         = CACHE_INVALID
@@ -859,8 +857,8 @@ contains
     ! Check to see if there is multipole data present at this energy
     use_mp = .false.
     if (this % mp_present) then
-      if (E >= this % multipole % start_E .and. &
-           E <= this % multipole % end_E) then
+      if (E >= this % multipole % E_min .and. &
+           E <= this % multipole % E_max) then
         use_mp = .true.
       end if
     end if
@@ -868,9 +866,10 @@ contains
     ! Evaluate multipole or interpolate
     if (use_mp) then
       ! Call multipole kernel
-      call multipole_eval(this % multipole, E, sqrtkT, sig_t, sig_a, sig_f)
+      call multipole_eval(this % multipole, E, sqrtkT, sig_s, sig_a, sig_f)
 
-      micro_xs % total = sig_t
+      micro_xs % total = sig_s + sig_a
+      micro_xs % elastic = sig_s
       micro_xs % absorption = sig_a
       micro_xs % fission = sig_f
 
@@ -1076,9 +1075,6 @@ contains
         micro_xs % elastic = (ONE - f) * rx % xs(i_temp, i_grid) + &
              f * rx % xs(i_temp, i_grid + 1)
       end associate
-    else
-      ! For multipole, elastic is total - absorption
-      micro_xs % elastic = micro_xs % total - micro_xs % absorption
     end if
   end subroutine nuclide_calculate_elastic_xs
 
@@ -1130,7 +1126,7 @@ contains
 ! sections in the resolved resonance regions
 !===============================================================================
 
-  subroutine multipole_eval(multipole, E, sqrtkT, sig_t, sig_a, sig_f)
+  subroutine multipole_eval(multipole, E, sqrtkT, sig_s, sig_a, sig_f)
     type(MultipoleArray), intent(in) :: multipole ! The windowed multipole
                                                   !  object to process.
     real(8), intent(in)              :: E         ! The energy at which to
@@ -1138,7 +1134,7 @@ contains
     real(8), intent(in)              :: sqrtkT    ! The temperature in the form
                                                   !  sqrt(kT), at which
                                                   !  to evaluate the XS.
-    real(8), intent(out)             :: sig_t     ! Total cross section
+    real(8), intent(out)             :: sig_s     ! Scattering cross section
     real(8), intent(out)             :: sig_a     ! Absorption cross section
     real(8), intent(out)             :: sig_f     ! Fission cross section
     complex(8) :: psi_chi  ! The value of the psi-chi function for the
@@ -1146,7 +1142,6 @@ contains
     complex(8) :: c_temp   ! complex temporary variable
     complex(8) :: w_val    ! The faddeeva function evaluated at Z
     complex(8) :: Z        ! sqrt(atomic weight ratio / kT) * (sqrt(E) - pole)
-    complex(8) :: sig_t_factor(multipole % num_l)
     real(8) :: broadened_polynomials(multipole % fit_order + 1)
     real(8) :: sqrtE       ! sqrt(E), eV
     real(8) :: invE        ! 1/E, eV
@@ -1166,18 +1161,13 @@ contains
     invE = ONE / E
 
     ! Locate us.
-    i_window = floor((sqrtE - sqrt(multipole % start_E)) / multipole % spacing &
+    i_window = floor((sqrtE - sqrt(multipole % E_min)) / multipole % spacing &
          + ONE)
-    startw = multipole % w_start(i_window)
-    endw = multipole % w_end(i_window)
-
-    ! Fill in factors.
-    if (startw <= endw) then
-      call compute_sig_t_factor(multipole, sqrtE, sig_t_factor)
-    end if
+    startw = multipole % windows(1, i_window)
+    endw = multipole % windows(2, i_window)
 
     ! Initialize the ouptut cross sections.
-    sig_t = ZERO
+    sig_s = ZERO
     sig_a = ZERO
     sig_f = ZERO
 
@@ -1190,7 +1180,7 @@ contains
       call broaden_wmp_polynomials(E, dopp, multipole % fit_order + 1, &
            broadened_polynomials)
       do i_poly = 1, multipole % fit_order+1
-        sig_t = sig_t + multipole % curvefit(FIT_T, i_poly, i_window) &
+        sig_s = sig_s + multipole % curvefit(FIT_S, i_poly, i_window) &
              * broadened_polynomials(i_poly)
         sig_a = sig_a + multipole % curvefit(FIT_A, i_poly, i_window) &
              * broadened_polynomials(i_poly)
@@ -1202,7 +1192,7 @@ contains
     else ! Evaluate as if it were a polynomial
       temp = invE
       do i_poly = 1, multipole % fit_order+1
-        sig_t = sig_t + multipole % curvefit(FIT_T, i_poly, i_window) * temp
+        sig_s = sig_s + multipole % curvefit(FIT_S, i_poly, i_window) * temp
         sig_a = sig_a + multipole % curvefit(FIT_A, i_poly, i_window) * temp
         if (multipole % fissionable) then
           sig_f = sig_f + multipole % curvefit(FIT_F, i_poly, i_window) * temp
@@ -1219,21 +1209,10 @@ contains
       do i_pole = startw, endw
         psi_chi = -ONEI / (multipole % data(MP_EA, i_pole) - sqrtE)
         c_temp = psi_chi / E
-        if (multipole % formalism == FORM_MLBW) then
-          sig_t = sig_t + real(multipole % data(MLBW_RT, i_pole) * c_temp * &
-                               sig_t_factor(multipole % l_value(i_pole))) &
-                        + real(multipole % data(MLBW_RX, i_pole) * c_temp)
-          sig_a = sig_a + real(multipole % data(MLBW_RA, i_pole) * c_temp)
-          if (multipole % fissionable) then
-            sig_f = sig_f + real(multipole % data(MLBW_RF, i_pole) * c_temp)
-          end if
-        else if (multipole % formalism == FORM_RM) then
-          sig_t = sig_t + real(multipole % data(RM_RT, i_pole) * c_temp * &
-                               sig_t_factor(multipole % l_value(i_pole)))
-          sig_a = sig_a + real(multipole % data(RM_RA, i_pole) * c_temp)
-          if (multipole % fissionable) then
-            sig_f = sig_f + real(multipole % data(RM_RF, i_pole) * c_temp)
-          end if
+        sig_s = sig_s + real(multipole % data(MP_RS, i_pole) * c_temp)
+        sig_a = sig_a + real(multipole % data(MP_RA, i_pole) * c_temp)
+        if (multipole % fissionable) then
+          sig_f = sig_f + real(multipole % data(MP_RF, i_pole) * c_temp)
         end if
       end do
     else
@@ -1243,21 +1222,10 @@ contains
         do i_pole = startw, endw
           Z = (sqrtE - multipole % data(MP_EA, i_pole)) * dopp
           w_val = faddeeva(Z) * dopp * invE * SQRT_PI
-          if (multipole % formalism == FORM_MLBW) then
-            sig_t = sig_t + real((multipole % data(MLBW_RT, i_pole) * &
-                            sig_t_factor(multipole % l_value(i_pole)) + &
-                            multipole % data(MLBW_RX, i_pole)) * w_val)
-            sig_a = sig_a + real(multipole % data(MLBW_RA, i_pole) * w_val)
-            if (multipole % fissionable) then
-              sig_f = sig_f + real(multipole % data(MLBW_RF, i_pole) * w_val)
-            end if
-          else if (multipole % formalism == FORM_RM) then
-            sig_t = sig_t + real(multipole % data(RM_RT, i_pole) * w_val * &
-                                 sig_t_factor(multipole % l_value(i_pole)))
-            sig_a = sig_a + real(multipole % data(RM_RA, i_pole) * w_val)
-            if (multipole % fissionable) then
-              sig_f = sig_f + real(multipole % data(RM_RF, i_pole) * w_val)
-            end if
+          sig_s = sig_s + real(multipole % data(MP_RS, i_pole) * w_val)
+          sig_a = sig_a + real(multipole % data(MP_RA, i_pole) * w_val)
+          if (multipole % fissionable) then
+            sig_f = sig_f + real(multipole % data(MP_RF, i_pole) * w_val)
           end if
         end do
       end if
@@ -1270,7 +1238,7 @@ contains
 ! temperature.
 !===============================================================================
 
-  subroutine multipole_deriv_eval(multipole, E, sqrtkT, sig_t, sig_a, sig_f)
+  subroutine multipole_deriv_eval(multipole, E, sqrtkT, sig_s, sig_a, sig_f)
     type(MultipoleArray), intent(in) :: multipole ! The windowed multipole
                                                   !  object to process.
     real(8), intent(in)              :: E         ! The energy at which to
@@ -1278,12 +1246,11 @@ contains
     real(8), intent(in)              :: sqrtkT    ! The temperature in the form
                                                   !  sqrt(kT), at which to
                                                   !  evaluate the XS.
-    real(8), intent(out)             :: sig_t     ! Total cross section
+    real(8), intent(out)             :: sig_s     ! Scattering cross section
     real(8), intent(out)             :: sig_a     ! Absorption cross section
     real(8), intent(out)             :: sig_f     ! Fission cross section
     complex(8) :: w_val    ! The faddeeva function evaluated at Z
     complex(8) :: Z        ! sqrt(atomic weight ratio / kT) * (sqrt(E) - pole)
-    complex(8) :: sig_t_factor(multipole % num_l)
     real(8) :: sqrtE       ! sqrt(E), eV
     real(8) :: invE        ! 1/E, eV
     real(8) :: dopp        ! sqrt(atomic weight ratio / kT)
@@ -1305,18 +1272,13 @@ contains
          &derivatives are not implemented for 0 Kelvin cross sections.")
 
     ! Locate us
-    i_window = floor((sqrtE - sqrt(multipole % start_E)) / multipole % spacing &
+    i_window = floor((sqrtE - sqrt(multipole % E_min)) / multipole % spacing &
          + ONE)
-    startw = multipole % w_start(i_window)
-    endw = multipole % w_end(i_window)
-
-    ! Fill in factors.
-    if (startw <= endw) then
-      call compute_sig_t_factor(multipole, sqrtE, sig_t_factor)
-    end if
+    startw = multipole % windows(1, i_window)
+    endw = multipole % windows(2, i_window)
 
     ! Initialize the ouptut cross sections.
-    sig_t = ZERO
+    sig_s = ZERO
     sig_a = ZERO
     sig_f = ZERO
 
@@ -1333,61 +1295,18 @@ contains
       do i_pole = startw, endw
         Z = (sqrtE - multipole % data(MP_EA, i_pole)) * dopp
         w_val = -invE * SQRT_PI * HALF * w_derivative(Z, 2)
-        if (multipole % formalism == FORM_MLBW) then
-          sig_t = sig_t + real((multipole % data(MLBW_RT, i_pole) * &
-                          sig_t_factor(multipole%l_value(i_pole)) + &
-                          multipole % data(MLBW_RX, i_pole)) * w_val)
-          sig_a = sig_a + real(multipole % data(MLBW_RA, i_pole) * w_val)
-          if (multipole % fissionable) then
-            sig_f = sig_f + real(multipole % data(MLBW_RF, i_pole) * w_val)
-          end if
-        else if (multipole % formalism == FORM_RM) then
-          sig_t = sig_t + real(multipole % data(RM_RT, i_pole) * w_val * &
-                               sig_t_factor(multipole % l_value(i_pole)))
-          sig_a = sig_a + real(multipole % data(RM_RA, i_pole) * w_val)
-          if (multipole % fissionable) then
-            sig_f = sig_f + real(multipole % data(RM_RF, i_pole) * w_val)
-          end if
+        sig_s = sig_s + real(multipole % data(MP_RS, i_pole) * w_val)
+        sig_a = sig_a + real(multipole % data(MP_RA, i_pole) * w_val)
+        if (multipole % fissionable) then
+          sig_f = sig_f + real(multipole % data(MP_RF, i_pole) * w_val)
         end if
       end do
-      sig_t = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sig_t
+      sig_s = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sig_s
       sig_a = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sig_a
       sig_f = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sig_f
     end if
   end subroutine multipole_deriv_eval
 
-!===============================================================================
-! COMPUTE_SIG_T_FACTOR calculates the sig_t_factor, a factor inside of the sig_t
-! equation not present in the sig_a and sig_f equations.
-!===============================================================================
-
-  subroutine compute_sig_t_factor(multipole, sqrtE, sig_t_factor)
-    type(MultipoleArray), intent(in)  :: multipole
-    real(8),              intent(in)  :: sqrtE
-    complex(8),           intent(out) :: sig_t_factor(multipole % num_l)
-
-    integer :: iL
-    real(8) :: twophi(multipole % num_l)
-    real(8) :: arg
-
-    do iL = 1, multipole % num_l
-      twophi(iL) = multipole % pseudo_k0RS(iL) * sqrtE
-      if (iL == 2) then
-        twophi(iL) = twophi(iL) - atan(twophi(iL))
-      else if (iL == 3) then
-        arg = 3.0_8 * twophi(iL) / (3.0_8 - twophi(iL)**2)
-        twophi(iL) = twophi(iL) - atan(arg)
-      else if (iL == 4) then
-        arg = twophi(iL) * (15.0_8 - twophi(iL)**2) &
-             / (15.0_8 - 6.0_8 * twophi(iL)**2)
-        twophi(iL) = twophi(iL) - atan(arg)
-      end if
-    end do
-
-    twophi = 2.0_8 * twophi
-    sig_t_factor = cmplx(cos(twophi), -sin(twophi), KIND=8)
-  end subroutine compute_sig_t_factor
-
 !===============================================================================
 ! 0K_ELASTIC_XS determines the microscopic 0K elastic cross section
 ! for a given nuclide at the trial relative energy used in resonance scattering

src/physics.F90

@@ -526,8 +526,8 @@ contains
     ! Check to see if we are in a windowed multipole range.  WMP only supports
     ! the first fission reaction.
     if (nuc % mp_present) then
-      if (E >= nuc % multipole % start_E .and. &
-           E <= nuc % multipole % end_E) then
+      if (E >= nuc % multipole % E_min .and. &
+           E <= nuc % multipole % E_max) then
         i_reaction = nuc % index_fission(1)
         return
       end if

src/tallies/tally.F90

@@ -3009,7 +3009,7 @@ contains
     integer :: l
     logical :: scoring_diff_nuclide
     real(8) :: flux_deriv
-    real(8) :: dsig_t, dsig_a, dsig_f, cum_dsig
+    real(8) :: dsig_s, dsig_a, dsig_f, cum_dsig
 
     if (score == ZERO) return
 
@@ -3250,17 +3250,19 @@ contains
                   if (mat % nuclide(l) == p % event_nuclide) exit
                 end do
 
-                dsig_t = ZERO
+                dsig_s = ZERO
+                dsig_a = ZERO
                 associate (nuc => nuclides(p % event_nuclide))
                   if (nuc % mp_present .and. &
-                       p % last_E >= nuc % multipole % start_E .and. &
-                       p % last_E <= nuc % multipole % end_E) then
+                       p % last_E >= nuc % multipole % E_min .and. &
+                       p % last_E <= nuc % multipole % E_max) then
                     call multipole_deriv_eval(nuc % multipole, p % last_E, &
-                         p % sqrtkT, dsig_t, dsig_a, dsig_f)
+                         p % sqrtkT, dsig_s, dsig_a, dsig_f)
                   end if
                 end associate
                 score = score * (flux_deriv &
-                     + dsig_t * mat % atom_density(l) / material_xs % total)
+                     + (dsig_s + dsig_a) * mat % atom_density(l) &
+                     / material_xs % total)
               end associate
             else
               score = score * flux_deriv
@@ -3276,18 +3278,16 @@ contains
                   if (mat % nuclide(l) == p % event_nuclide) exit
                 end do
 
-                dsig_t = ZERO
-                dsig_a = ZERO
+                dsig_s = ZERO
                 associate (nuc => nuclides(p % event_nuclide))
                   if (nuc % mp_present .and. &
-                       p % last_E >= nuc % multipole % start_E .and. &
-                       p % last_E <= nuc % multipole % end_E) then
+                       p % last_E >= nuc % multipole % E_min .and. &
+                       p % last_E <= nuc % multipole % E_max) then
                     call multipole_deriv_eval(nuc % multipole, p % last_E, &
-                         p % sqrtkT, dsig_t, dsig_a, dsig_f)
+                         p % sqrtkT, dsig_s, dsig_a, dsig_f)
                   end if
                 end associate
-                score = score * (flux_deriv + (dsig_t - dsig_a) &
-                     * mat % atom_density(l) / &
+                score = score * (flux_deriv + dsig_s * mat % atom_density(l) / &
                      (material_xs % total - material_xs % absorption))
               end associate
             else
@@ -3306,10 +3306,10 @@ contains
                 dsig_a = ZERO
                 associate (nuc => nuclides(p % event_nuclide))
                   if (nuc % mp_present .and. &
-                       p % last_E >= nuc % multipole % start_E .and. &
-                       p % last_E <= nuc % multipole % end_E) then
+                       p % last_E >= nuc % multipole % E_min .and. &
+                       p % last_E <= nuc % multipole % E_max) then
                     call multipole_deriv_eval(nuc % multipole, p % last_E, &
-                         p % sqrtkT, dsig_t, dsig_a, dsig_f)
+                         p % sqrtkT, dsig_s, dsig_a, dsig_f)
                   end if
                 end associate
                 score = score * (flux_deriv + dsig_a * mat % atom_density(l) &
@@ -3331,10 +3331,10 @@ contains
                 dsig_f = ZERO
                 associate (nuc => nuclides(p % event_nuclide))
                   if (nuc % mp_present .and. &
-                       p % last_E >= nuc % multipole % start_E .and. &
-                       p % last_E <= nuc % multipole % end_E) then
+                       p % last_E >= nuc % multipole % E_min .and. &
+                       p % last_E <= nuc % multipole % E_max) then
                     call multipole_deriv_eval(nuc % multipole, p % last_E, &
-                         p % sqrtkT, dsig_t, dsig_a, dsig_f)
+                         p % sqrtkT, dsig_s, dsig_a, dsig_f)
                   end if
                 end associate
                 score = score * (flux_deriv &
@@ -3356,10 +3356,10 @@ contains
                 dsig_f = ZERO
                 associate (nuc => nuclides(p % event_nuclide))
                   if (nuc % mp_present .and. &
-                       p % last_E >= nuc % multipole % start_E .and. &
-                       p % last_E <= nuc % multipole % end_E) then
+                       p % last_E >= nuc % multipole % E_min .and. &
+                       p % last_E <= nuc % multipole % E_max) then
                     call multipole_deriv_eval(nuc % multipole, p % last_E, &
-                         p % sqrtkT, dsig_t, dsig_a, dsig_f)
+                         p % sqrtkT, dsig_s, dsig_a, dsig_f)
                   end if
                 end associate
                 score = score * (flux_deriv &
@@ -3392,12 +3392,13 @@ contains
                 do l = 1, mat % n_nuclides
                   associate (nuc => nuclides(mat % nuclide(l)))
                     if (nuc % mp_present .and. &
-                         p % last_E >= nuc % multipole % start_E .and. &
-                         p % last_E <= nuc % multipole % end_E .and. &
+                         p % last_E >= nuc % multipole % E_min .and. &
+                         p % last_E <= nuc % multipole % E_max .and. &
                          micro_xs(mat % nuclide(l)) % total > ZERO) then
                       call multipole_deriv_eval(nuc % multipole, p % last_E, &
-                           p % sqrtkT, dsig_t, dsig_a, dsig_f)
-                      cum_dsig = cum_dsig + dsig_t * mat % atom_density(l)
+                           p % sqrtkT, dsig_s, dsig_a, dsig_f)
+                      cum_dsig = cum_dsig + (dsig_s + dsig_a) &
+                           * mat % atom_density(l)
                     end if
                   end associate
                 end do
@@ -3406,17 +3407,18 @@ contains
                    + cum_dsig / material_xs % total)
             else if (materials(p % material) % id() == deriv % diff_material &
                  .and. material_xs % total > ZERO) then
-              dsig_t = ZERO
+              dsig_s = ZERO
+              dsig_a = ZERO
               associate (nuc => nuclides(i_nuclide))
                 if (nuc % mp_present .and. &
-                     p % last_E >= nuc % multipole % start_E .and. &
-                     p % last_E <= nuc % multipole % end_E) then
+                     p % last_E >= nuc % multipole % E_min .and. &
+                     p % last_E <= nuc % multipole % E_max) then
                   call multipole_deriv_eval(nuc % multipole, p % last_E, &
-                       p % sqrtkT, dsig_t, dsig_a, dsig_f)
+                       p % sqrtkT, dsig_s, dsig_a, dsig_f)
                 end if
               end associate
               score = score * (flux_deriv &
-                   + dsig_t / micro_xs(i_nuclide) % total)
+                   + (dsig_s + dsig_a) / micro_xs(i_nuclide) % total)
             else
               score = score * flux_deriv
             end if
@@ -3430,14 +3432,13 @@ contains
                 do l = 1, mat % n_nuclides
                   associate (nuc => nuclides(mat % nuclide(l)))
                     if (nuc % mp_present .and. &
-                         p % last_E >= nuc % multipole % start_E .and. &
-                         p % last_E <= nuc % multipole % end_E .and. &
+                         p % last_E >= nuc % multipole % E_min .and. &
+                         p % last_E <= nuc % multipole % E_max .and. &
                          (micro_xs(mat % nuclide(l)) % total &
                          - micro_xs(mat % nuclide(l)) % absorption) > ZERO) then
                       call multipole_deriv_eval(nuc % multipole, p % last_E, &
-                           p % sqrtkT, dsig_t, dsig_a, dsig_f)
-                      cum_dsig = cum_dsig &
-                           + (dsig_t - dsig_a) * mat % atom_density(l)
+                           p % sqrtkT, dsig_s, dsig_a, dsig_f)
+                      cum_dsig = cum_dsig + dsig_s * mat % atom_density(l)
                     end if
                   end associate
                 end do
@@ -3447,17 +3448,16 @@ contains
             else if ( materials(p % material) % id() == deriv % diff_material &
                  .and. (material_xs % total - material_xs % absorption) > ZERO)&
                  then
-              dsig_t = ZERO
-              dsig_a = ZERO
+              dsig_s = ZERO
               associate (nuc => nuclides(i_nuclide))
                 if (nuc % mp_present .and. &
-                     p % last_E >= nuc % multipole % start_E .and. &
-                     p % last_E <= nuc % multipole % end_E) then
+                     p % last_E >= nuc % multipole % E_min .and. &
+                     p % last_E <= nuc % multipole % E_max) then
                   call multipole_deriv_eval(nuc % multipole, p % last_E, &
-                       p % sqrtkT, dsig_t, dsig_a, dsig_f)
+                       p % sqrtkT, dsig_s, dsig_a, dsig_f)
                 end if
               end associate
-              score = score * (flux_deriv + (dsig_t - dsig_a) &
+              score = score * (flux_deriv + dsig_s &
                    / (micro_xs(i_nuclide) % total &
                    - micro_xs(i_nuclide) % absorption))
             else
@@ -3473,11 +3473,11 @@ contains
                 do l = 1, mat % n_nuclides
                   associate (nuc => nuclides(mat % nuclide(l)))
                     if (nuc % mp_present .and. &
-                         p % last_E >= nuc % multipole % start_E .and. &
-                         p % last_E <= nuc % multipole % end_E .and. &
+                         p % last_E >= nuc % multipole % E_min .and. &
+                         p % last_E <= nuc % multipole % E_max .and. &
                          micro_xs(mat % nuclide(l)) % absorption > ZERO) then
                       call multipole_deriv_eval(nuc % multipole, p % last_E, &
-                           p % sqrtkT, dsig_t, dsig_a, dsig_f)
+                           p % sqrtkT, dsig_s, dsig_a, dsig_f)
                       cum_dsig = cum_dsig + dsig_a * mat % atom_density(l)
                     end if
                   end associate
@@ -3490,10 +3490,10 @@ contains
               dsig_a = ZERO
               associate (nuc => nuclides(i_nuclide))
                 if (nuc % mp_present .and. &
-                     p % last_E >= nuc % multipole % start_E .and. &
-                     p % last_E <= nuc % multipole % end_E) then
+                     p % last_E >= nuc % multipole % E_min .and. &
+                     p % last_E <= nuc % multipole % E_max) then
                   call multipole_deriv_eval(nuc % multipole, p % last_E, &
-                       p % sqrtkT, dsig_t, dsig_a, dsig_f)
+                       p % sqrtkT, dsig_s, dsig_a, dsig_f)
                 end if
               end associate
               score = score * (flux_deriv &
@@ -3511,11 +3511,11 @@ contains
                 do l = 1, mat % n_nuclides
                   associate (nuc => nuclides(mat % nuclide(l)))
                     if (nuc % mp_present .and. &
-                         p % last_E >= nuc % multipole % start_E .and. &
-                         p % last_E <= nuc % multipole % end_E .and. &
+                         p % last_E >= nuc % multipole % E_min .and. &
+                         p % last_E <= nuc % multipole % E_max .and. &
                          micro_xs(mat % nuclide(l)) % fission > ZERO) then
                       call multipole_deriv_eval(nuc % multipole, p % last_E, &
-                           p % sqrtkT, dsig_t, dsig_a, dsig_f)
+                           p % sqrtkT, dsig_s, dsig_a, dsig_f)
                       cum_dsig = cum_dsig + dsig_f * mat % atom_density(l)
                     end if
                   end associate
@@ -3528,10 +3528,10 @@ contains
               dsig_f = ZERO
               associate (nuc => nuclides(i_nuclide))
                 if (nuc % mp_present .and. &
-                     p % last_E >= nuc % multipole % start_E .and. &
-                     p % last_E <= nuc % multipole % end_E) then
+                     p % last_E >= nuc % multipole % E_min .and. &
+                     p % last_E <= nuc % multipole % E_max) then
                   call multipole_deriv_eval(nuc % multipole, p % last_E, &
-                       p % sqrtkT, dsig_t, dsig_a, dsig_f)
+                       p % sqrtkT, dsig_s, dsig_a, dsig_f)
                 end if
               end associate
               score = score * (flux_deriv &
@@ -3549,11 +3549,11 @@ contains
                 do l = 1, mat % n_nuclides
                   associate (nuc => nuclides(mat % nuclide(l)))
                     if (nuc % mp_present .and. &
-                         p % last_E >= nuc % multipole % start_E .and. &
-                         p % last_E <= nuc % multipole % end_E .and. &
+                         p % last_E >= nuc % multipole % E_min .and. &
+                         p % last_E <= nuc % multipole % E_max .and. &
                          micro_xs(mat % nuclide(l)) % nu_fission > ZERO) then
                       call multipole_deriv_eval(nuc % multipole, p % last_E, &
-                           p % sqrtkT, dsig_t, dsig_a, dsig_f)
+                           p % sqrtkT, dsig_s, dsig_a, dsig_f)
                       cum_dsig = cum_dsig + dsig_f * mat % atom_density(l) &
                            * micro_xs(mat % nuclide(l)) % nu_fission &
                            / micro_xs(mat % nuclide(l)) % fission
@@ -3568,10 +3568,10 @@ contains
               dsig_f = ZERO
               associate (nuc => nuclides(i_nuclide))
                 if (nuc % mp_present .and. &
-                     p % last_E >= nuc % multipole % start_E .and. &
-                     p % last_E <= nuc % multipole % end_E) then
+                     p % last_E >= nuc % multipole % E_min .and. &
+                     p % last_E <= nuc % multipole % E_max) then
                   call multipole_deriv_eval(nuc % multipole, p % last_E, &
-                       p % sqrtkT, dsig_t, dsig_a, dsig_f)
+                       p % sqrtkT, dsig_s, dsig_a, dsig_f)
                 end if
               end associate
               score = score * (flux_deriv &
@@ -3603,7 +3603,7 @@ contains
     real(8),        intent(in) :: distance ! Neutron flight distance
 
     integer :: i, l
-    real(8) :: dsig_t, dsig_a, dsig_f
+    real(8) :: dsig_s, dsig_a, dsig_f
 
     ! A void material cannot be perturbed so it will not affect flux derivatives
     if (p % material == MATERIAL_VOID) return
@@ -3640,15 +3640,15 @@ contains
               do l=1, mat % n_nuclides
                 associate (nuc => nuclides(mat % nuclide(l)))
                   if (nuc % mp_present .and. &
-                       p % E >= nuc % multipole % start_E .and. &
-                       p % E <= nuc % multipole % end_E) then
+                       p % E >= nuc % multipole % E_min .and. &
+                       p % E <= nuc % multipole % E_max) then
                     ! phi is proportional to e^(-Sigma_tot * dist)
                     ! (1 / phi) * (d_phi / d_T) = - (d_Sigma_tot / d_T) * dist
                     ! (1 / phi) * (d_phi / d_T) = - N (d_sigma_tot / d_T) * dist
                     call multipole_deriv_eval(nuc % multipole, p % E, &
-                         p % sqrtkT, dsig_t, dsig_a, dsig_f)
+                         p % sqrtkT, dsig_s, dsig_a, dsig_f)
                     deriv % flux_deriv = deriv % flux_deriv &
-                         - distance * dsig_t * mat % atom_density(l)
+                         - distance * (dsig_s + dsig_a) * mat % atom_density(l)
                   end if
                 end associate
               end do
@@ -3679,7 +3679,7 @@ contains
     type(Particle), intent(in) :: p
 
     integer :: i, j, l
-    real(8) :: dsig_t, dsig_a, dsig_f
+    real(8) :: dsig_s, dsig_a, dsig_f
 
     ! A void material cannot be perturbed so it will not affect flux derivatives
     if (p % material == MATERIAL_VOID) return
@@ -3727,14 +3727,14 @@ contains
                 associate (nuc => nuclides(mat % nuclide(l)))
                   if (mat % nuclide(l) == p % event_nuclide .and. &
                        nuc % mp_present .and. &
-                       p % last_E >= nuc % multipole % start_E .and. &
-                       p % last_E <= nuc % multipole % end_E) then
+                       p % last_E >= nuc % multipole % E_min .and. &
+                       p % last_E <= nuc % multipole % E_max) then
                     ! phi is proportional to Sigma_s
                     ! (1 / phi) * (d_phi / d_T) = (d_Sigma_s / d_T) / Sigma_s
                     ! (1 / phi) * (d_phi / d_T) = (d_sigma_s / d_T) / sigma_s
                     call multipole_deriv_eval(nuc % multipole, p % last_E, &
-                         p % sqrtkT, dsig_t, dsig_a, dsig_f)
-                    deriv % flux_deriv = deriv % flux_deriv + (dsig_t - dsig_a)&
+                         p % sqrtkT, dsig_s, dsig_a, dsig_f)
+                    deriv % flux_deriv = deriv % flux_deriv + dsig_s&
                          / (micro_xs(mat % nuclide(l)) % total &
                          - micro_xs(mat % nuclide(l)) % absorption)
                     ! Note that this is an approximation!  The real scattering

tests/regression_tests/diff_tally/results_true.dat

@@ -1,27 +1,27 @@
 d_material,d_nuclide,d_variable,score,mean,std. dev.
-3,,density,flux,-4.5951022e+00,4.0201889e-01
-3,,density,flux,-9.9630192e+00,1.8897688e+00
-1,,density,flux,-4.2074596e-01,6.8569557e-02
-1,,density,flux,-2.8861314e-01,1.0032157e-01
-1,O16,nuclide_density,flux,-1.7146771e+01,1.8740496e+01
-1,O16,nuclide_density,flux,-1.4411733e+01,1.5986415e+01
-1,U235,nuclide_density,flux,-3.1082590e+03,3.5930199e+02
-1,U235,nuclide_density,flux,-2.5816702e+03,2.4275334e+02
-1,,temperature,flux,-2.2847162e-06,3.5489476e-04
-1,,temperature,flux,-8.5566667e-05,4.3023777e-04
-3,,density,total,-1.5582276e+00,2.8496759e-01
-3,,density,absorption,1.5711173e-01,1.7341244e-01
-3,,density,scatter,-1.7153393e+00,1.9431809e-01
-3,,density,fission,2.0152323e-01,7.4786095e-02
-3,,density,nu-fission,4.8687131e-01,1.8265878e-01
-3,,density,total,2.2943608e-01,8.6539576e-02
-3,,density,absorption,2.4646330e-01,8.6186628e-02
-3,,density,scatter,-1.7027217e-02,3.1732345e-03
-3,,density,fission,2.1668870e-01,7.5183313e-02
-3,,density,nu-fission,5.2774399e-01,1.8320780e-01
-3,,density,total,1.3300225e+01,4.2626034e+00
-3,,density,absorption,2.3903751e-01,8.9355264e-02
-3,,density,scatter,1.3061188e+01,4.1751261e+00
+3,,density,flux,-4.7291290e+00,8.8503902e-01
+3,,density,flux,-1.0533184e+01,3.0256001e+00
+1,,density,flux,-4.9634223e-01,1.3190338e-01
+1,,density,flux,-4.7458622e-01,5.6426916e-02
+1,O16,nuclide_density,flux,-1.4897399e+01,1.3583122e+01
+1,O16,nuclide_density,flux,-2.2389753e+01,1.5574833e+01
+1,U235,nuclide_density,flux,-2.8858701e+03,4.7033287e+02
+1,U235,nuclide_density,flux,-2.4203468e+03,3.2197255e+02
+1,,temperature,flux,-1.1195282e-04,3.5426553e-04
+1,,temperature,flux,-4.4294257e-04,5.4687257e-04
+3,,density,total,-1.5555916e+00,5.3353204e-01
+3,,density,absorption,1.4925823e-01,2.3501173e-01
+3,,density,scatter,-1.7048498e+00,3.3363896e-01
+3,,density,fission,1.3210039e-01,1.6377610e-01
+3,,density,nu-fission,3.1732288e-01,3.9893191e-01
+3,,density,total,1.4744448e-01,2.0435075e-01
+3,,density,absorption,1.6665433e-01,1.9696442e-01
+3,,density,scatter,-1.9209851e-02,7.8974770e-03
+3,,density,fission,1.4878594e-01,1.6616312e-01
+3,,density,nu-fission,3.6225815e-01,4.0486140e-01
+3,,density,total,1.2662462e+01,6.2596594e+00
+3,,density,absorption,2.2864288e-01,1.3031728e-01
+3,,density,scatter,1.2433820e+01,6.1294912e+00
 3,,density,fission,0.0000000e+00,0.0000000e+00
 3,,density,nu-fission,0.0000000e+00,0.0000000e+00
 3,,density,total,0.0000000e+00,0.0000000e+00
@@ -29,19 +29,19 @@ d_material,d_nuclide,d_variable,score,mean,std. dev.
 3,,density,scatter,0.0000000e+00,0.0000000e+00
 3,,density,fission,0.0000000e+00,0.0000000e+00
 3,,density,nu-fission,0.0000000e+00,0.0000000e+00
-1,,density,total,3.3373495e-01,4.6827421e-02
-1,,density,absorption,6.7487092e-04,6.5232977e-03
-1,,density,scatter,3.3306008e-01,4.0307583e-02
-1,,density,fission,-1.8402784e-03,4.2206212e-03
-1,,density,nu-fission,-3.8312038e-03,1.0337874e-02
-1,,density,total,-2.4831544e-04,5.5020059e-03
-1,,density,absorption,-4.0311629e-03,5.0266899e-03
-1,,density,scatter,3.7828475e-03,4.8043964e-04
-1,,density,fission,-3.6379968e-03,4.1821819e-03
-1,,density,nu-fission,-8.8266852e-03,1.0193373e-02
-1,,density,total,-3.6848185e-01,1.9314902e-01
-1,,density,absorption,-7.3640407e-03,4.0184508e-03
-1,,density,scatter,-3.6111780e-01,1.8936612e-01
+1,,density,total,2.9404312e-01,7.0116815e-02
+1,,density,absorption,-1.1851195e-03,1.6943078e-03
+1,,density,scatter,2.9522824e-01,6.8519873e-02
+1,,density,fission,-4.4865588e-03,2.3660449e-03
+1,,density,nu-fission,-1.0266898e-02,5.7956130e-03
+1,,density,total,-3.6990351e-03,3.5138725e-03
+1,,density,absorption,-7.0064347e-03,2.7205092e-03
+1,,density,scatter,3.3073996e-03,7.9802045e-04
+1,,density,fission,-6.3630708e-03,2.3832611e-03
+1,,density,nu-fission,-1.5466339e-02,5.8086647e-03
+1,,density,total,-5.5743331e-01,6.1625011e-02
+1,,density,absorption,-9.1430498e-03,4.0831822e-03
+1,,density,scatter,-5.4829026e-01,5.7756112e-02
 1,,density,fission,0.0000000e+00,0.0000000e+00
 1,,density,nu-fission,0.0000000e+00,0.0000000e+00
 1,,density,total,0.0000000e+00,0.0000000e+00
@@ -49,19 +49,19 @@ d_material,d_nuclide,d_variable,score,mean,std. dev.
 1,,density,scatter,0.0000000e+00,0.0000000e+00
 1,,density,fission,0.0000000e+00,0.0000000e+00
 1,,density,nu-fission,0.0000000e+00,0.0000000e+00
-1,O16,nuclide_density,total,3.9050701e+01,8.5131265e+00
-1,O16,nuclide_density,absorption,-7.7573286e-01,7.1182570e-01
-1,O16,nuclide_density,scatter,3.9826434e+01,7.8536529e+00
-1,O16,nuclide_density,fission,-4.2193390e-01,6.0091792e-01
-1,O16,nuclide_density,nu-fission,-1.0411912e+00,1.4662352e+00
-1,O16,nuclide_density,total,-5.8976625e-01,7.9559349e-01
-1,O16,nuclide_density,absorption,-5.0047331e-01,7.1666177e-01
-1,O16,nuclide_density,scatter,-8.9292946e-02,8.9937898e-02
-1,O16,nuclide_density,fission,-3.6634959e-01,6.1270658e-01
-1,O16,nuclide_density,nu-fission,-8.9338374e-01,1.4932014e+00
-1,O16,nuclide_density,total,-3.3511352e+00,2.2365725e+01
-1,O16,nuclide_density,absorption,1.5473875e-01,4.3692662e-01
-1,O16,nuclide_density,scatter,-3.5058740e+00,2.1928948e+01
+1,O16,nuclide_density,total,3.9360816e+01,6.1390337e+00
+1,O16,nuclide_density,absorption,-6.8169613e-01,4.5302455e-01
+1,O16,nuclide_density,scatter,4.0042512e+01,6.4069240e+00
+1,O16,nuclide_density,fission,-6.1275244e-01,2.5764187e-01
+1,O16,nuclide_density,nu-fission,-1.5090099e+00,6.3048791e-01
+1,O16,nuclide_density,total,-7.5872288e-01,3.4843447e-01
+1,O16,nuclide_density,absorption,-6.7382098e-01,3.1002267e-01
+1,O16,nuclide_density,scatter,-8.4901894e-02,6.5286818e-02
+1,O16,nuclide_density,fission,-5.5597011e-01,2.7085418e-01
+1,O16,nuclide_density,nu-fission,-1.3555055e+00,6.6014209e-01
+1,O16,nuclide_density,total,-1.5539605e+01,2.3345398e+01
+1,O16,nuclide_density,absorption,-8.5880168e-02,5.0677339e-01
+1,O16,nuclide_density,scatter,-1.5453725e+01,2.2838892e+01
 1,O16,nuclide_density,fission,0.0000000e+00,0.0000000e+00
 1,O16,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00
 1,O16,nuclide_density,total,0.0000000e+00,0.0000000e+00
@@ -69,19 +69,19 @@ d_material,d_nuclide,d_variable,score,mean,std. dev.
 1,O16,nuclide_density,scatter,0.0000000e+00,0.0000000e+00
 1,O16,nuclide_density,fission,0.0000000e+00,0.0000000e+00
 1,O16,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00
-1,U235,nuclide_density,total,-9.1643805e+02,1.8814046e+02
-1,U235,nuclide_density,absorption,1.8114000e+02,4.6345847e+01
-1,U235,nuclide_density,scatter,-1.0975781e+03,1.4310981e+02
-1,U235,nuclide_density,fission,2.8765342e+02,2.2586319e+01
-1,U235,nuclide_density,nu-fission,7.0175585e+02,5.5024716e+01
-1,U235,nuclide_density,total,4.6786282e+02,2.8271731e+01
-1,U235,nuclide_density,absorption,3.6134508e+02,2.8453765e+01
-1,U235,nuclide_density,scatter,1.0651774e+02,1.3120589e+00
-1,U235,nuclide_density,fission,2.8859249e+02,2.2522591e+01
-1,U235,nuclide_density,nu-fission,7.0430623e+02,5.4858487e+01
-1,U235,nuclide_density,total,-5.0586963e+03,6.7151365e+02
-1,U235,nuclide_density,absorption,-1.2315731e+02,1.8905184e+01
-1,U235,nuclide_density,scatter,-4.9355390e+03,6.5260939e+02
+1,U235,nuclide_density,total,-7.9766848e+02,2.3476486e+02
+1,U235,nuclide_density,absorption,2.2136536e+02,5.7727802e+01
+1,U235,nuclide_density,scatter,-1.0190338e+03,1.7827380e+02
+1,U235,nuclide_density,fission,3.0782358e+02,2.4835876e+01
+1,U235,nuclide_density,nu-fission,7.5106928e+02,6.0668638e+01
+1,U235,nuclide_density,total,4.9643664e+02,3.6425359e+01
+1,U235,nuclide_density,absorption,3.8860811e+02,3.5290173e+01
+1,U235,nuclide_density,scatter,1.0782853e+02,1.9206965e+00
+1,U235,nuclide_density,fission,3.0827647e+02,2.4266512e+01
+1,U235,nuclide_density,nu-fission,7.5228268e+02,5.9109738e+01
+1,U235,nuclide_density,total,-4.7231732e+03,7.5353357e+02
+1,U235,nuclide_density,absorption,-1.1580370e+02,1.8531790e+01
+1,U235,nuclide_density,scatter,-4.6073695e+03,7.3502263e+02
 1,U235,nuclide_density,fission,0.0000000e+00,0.0000000e+00
 1,U235,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00
 1,U235,nuclide_density,total,0.0000000e+00,0.0000000e+00
@@ -89,19 +89,19 @@ d_material,d_nuclide,d_variable,score,mean,std. dev.
 1,U235,nuclide_density,scatter,0.0000000e+00,0.0000000e+00
 1,U235,nuclide_density,fission,0.0000000e+00,0.0000000e+00
 1,U235,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00
-1,,temperature,total,8.6368265e-05,1.7363899e-04
-1,,temperature,absorption,1.9220262e-05,3.4988977e-05
-1,,temperature,scatter,6.7148003e-05,1.4635162e-04
-1,,temperature,fission,-5.2927287e-06,1.2158338e-05
-1,,temperature,nu-fission,-1.2897086e-05,2.9622912e-05
-1,,temperature,total,-4.2009323e-06,1.8263855e-05
-1,,temperature,absorption,-4.3037914e-06,1.6219907e-05
-1,,temperature,scatter,1.0285902e-07,2.0676492e-06
-1,,temperature,fission,-5.2953783e-06,1.2158786e-05
-1,,temperature,nu-fission,-1.2903531e-05,2.9624000e-05
-1,,temperature,total,-1.1160261e-04,6.1705782e-04
-1,,temperature,absorption,-2.0394397e-06,8.4055062e-06
-1,,temperature,scatter,-1.0956317e-04,6.0869070e-04
+1,,temperature,total,5.7228698e-05,1.7465295e-04
+1,,temperature,absorption,2.9495471e-05,3.1637786e-05
+1,,temperature,scatter,2.7733227e-05,1.4323068e-04
+1,,temperature,fission,-5.6710689e-06,1.2800905e-05
+1,,temperature,nu-fission,-1.3819116e-05,3.1189136e-05
+1,,temperature,total,-5.8315815e-06,1.8705738e-05
+1,,temperature,absorption,-5.3204426e-06,1.6688104e-05
+1,,temperature,scatter,-5.1113883e-07,2.0213875e-06
+1,,temperature,fission,-5.6723577e-06,1.2800214e-05
+1,,temperature,nu-fission,-1.3822264e-05,3.1187458e-05
+1,,temperature,total,-5.5283703e-04,7.6408422e-04
+1,,temperature,absorption,-6.2179157e-06,1.0365194e-05
+1,,temperature,scatter,-5.4661911e-04,7.5373133e-04
 1,,temperature,fission,0.0000000e+00,0.0000000e+00
 1,,temperature,nu-fission,0.0000000e+00,0.0000000e+00
 1,,temperature,total,0.0000000e+00,0.0000000e+00
@@ -109,68 +109,68 @@ d_material,d_nuclide,d_variable,score,mean,std. dev.
 1,,temperature,scatter,0.0000000e+00,0.0000000e+00
 1,,temperature,fission,0.0000000e+00,0.0000000e+00
 1,,temperature,nu-fission,0.0000000e+00,0.0000000e+00
-3,,density,absorption,1.3594219e-01,1.0122046e-01
-3,,density,absorption,2.8656499e-01,3.7923985e-02
-1,,density,absorption,-2.2106085e-03,9.8387075e-03
-1,,density,absorption,-3.3239716e-03,8.7768141e-03
-1,O16,nuclide_density,absorption,-1.2640173e+00,9.1596515e-01
-1,O16,nuclide_density,absorption,1.3146714e-01,9.9380272e-01
-1,U235,nuclide_density,absorption,1.5677452e+02,8.1085465e+01
-1,U235,nuclide_density,absorption,-1.4334804e+02,3.2100956e+01
-1,,temperature,absorption,-8.9633351e-06,3.5750907e-05
-1,,temperature,absorption,7.0140498e-07,1.1694324e-05
+3,,density,absorption,3.1512672e-02,2.4627803e-01
+3,,density,absorption,1.1452915e-01,1.5733776e-01
+1,,density,absorption,-1.4010827e-03,5.7414887e-03
+1,,density,absorption,-8.9116087e-03,7.6054344e-03
+1,O16,nuclide_density,absorption,-1.0238428e+00,3.5596763e-01
+1,O16,nuclide_density,absorption,-5.4386638e-01,6.3275232e-01
+1,U235,nuclide_density,absorption,1.9740455e+02,1.1770740e+02
+1,U235,nuclide_density,absorption,-1.3873202e+02,3.0030411e+01
+1,,temperature,absorption,5.0340072e-06,2.9911538e-05
+1,,temperature,absorption,-3.7625680e-05,2.4883845e-05
 3,,density,nu-fission,0.0000000e+00,0.0000000e+00
-3,,density,scatter,4.7213241e-01,1.6556808e-01
+3,,density,scatter,2.5895299e-01,3.3093831e-01
 3,,density,nu-fission,0.0000000e+00,0.0000000e+00
-3,,density,scatter,3.4067110e-02,2.1267271e-02
-3,,density,nu-fission,4.1265284e-01,1.3891654e-01
-3,,density,scatter,-2.1663022e+00,3.2444877e-01
-3,,density,nu-fission,4.4524898e-01,1.4098968e-01
-3,,density,scatter,-2.1357153e-02,1.7423446e-02
+3,,density,scatter,2.4380993e-02,2.4380993e-02
+3,,density,nu-fission,1.3982705e-01,4.0427820e-01
+3,,density,scatter,-1.8460573e+00,1.6126524e-01
+3,,density,nu-fission,1.7650022e-01,4.0560218e-01
+3,,density,scatter,3.9711501e-02,4.4405692e-02
 3,,density,nu-fission,0.0000000e+00,0.0000000e+00
-3,,density,scatter,8.6343785e+00,2.9292930e+00
+3,,density,scatter,7.8960553e+00,4.4763896e+00
 3,,density,nu-fission,0.0000000e+00,0.0000000e+00
 3,,density,scatter,0.0000000e+00,0.0000000e+00
 3,,density,nu-fission,0.0000000e+00,0.0000000e+00
-3,,density,scatter,4.3792818e+00,1.9915652e+00
+3,,density,scatter,4.6518780e+00,1.6862698e+00
 3,,density,nu-fission,0.0000000e+00,0.0000000e+00
 3,,density,scatter,0.0000000e+00,0.0000000e+00
 1,,density,nu-fission,0.0000000e+00,0.0000000e+00
-1,,density,scatter,-5.0479210e-03,9.5688761e-03
+1,,density,scatter,-1.4155379e-02,7.8989014e-03
 1,,density,nu-fission,0.0000000e+00,0.0000000e+00
-1,,density,scatter,4.4418327e-04,7.8121620e-04
-1,,density,nu-fission,-1.9240114e-02,1.4935352e-02
-1,,density,scatter,3.4099348e-01,2.8182893e-02
-1,,density,nu-fission,-2.4226858e-02,1.4481188e-02
-1,,density,scatter,1.1073847e-03,2.7478584e-04
+1,,density,scatter,1.0449664e-03,5.7043392e-04
+1,,density,nu-fission,-1.3723982e-02,1.5246898e-02
+1,,density,scatter,3.0959958e-01,5.7651589e-02
+1,,density,nu-fission,-1.9529210e-02,1.4804175e-02
+1,,density,scatter,5.1063052e-04,1.7748731e-03
 1,,density,nu-fission,0.0000000e+00,0.0000000e+00
-1,,density,scatter,-2.7695506e-01,1.3529978e-01
+1,,density,scatter,-3.4116700e-01,1.4433252e-01
 1,,density,nu-fission,0.0000000e+00,0.0000000e+00
 1,,density,scatter,0.0000000e+00,0.0000000e+00
 1,,density,nu-fission,0.0000000e+00,0.0000000e+00
-1,,density,scatter,-8.8202815e-02,7.8863942e-02
+1,,density,scatter,-2.0735470e-01,8.6865144e-02
 1,,density,nu-fission,0.0000000e+00,0.0000000e+00
 1,,density,scatter,0.0000000e+00,0.0000000e+00
 1,O16,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00
-1,O16,nuclide_density,scatter,1.1279170e-01,6.9836765e-02
-1,O16,nuclide_density,nu-fission,-1.6695105e+00,1.6777445e+00
-1,O16,nuclide_density,scatter,-1.0671543e-01,1.9127943e-01
+1,O16,nuclide_density,scatter,1.4017611e-01,9.5507870e-02
+1,O16,nuclide_density,nu-fission,-8.8199870e-01,1.6570982e+00
+1,O16,nuclide_density,scatter,-1.5442745e-01,1.3136597e-01
 1,O16,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00
 1,O16,nuclide_density,scatter,0.0000000e+00,0.0000000e+00
 1,O16,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00
 1,O16,nuclide_density,scatter,0.0000000e+00,0.0000000e+00
 1,U235,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00
-1,U235,nuclide_density,scatter,6.3235790e+00,3.9733167e+00
-1,U235,nuclide_density,nu-fission,6.2368225e+02,1.0172850e+02
-1,U235,nuclide_density,scatter,3.2601903e+01,6.2672072e+00
+1,U235,nuclide_density,scatter,8.1567193e+00,5.7077716e+00
+1,U235,nuclide_density,nu-fission,7.1482113e+02,1.4806072e+02
+1,U235,nuclide_density,scatter,5.2848809e+01,1.3486619e+01
 1,U235,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00
 1,U235,nuclide_density,scatter,0.0000000e+00,0.0000000e+00
 1,U235,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00
 1,U235,nuclide_density,scatter,0.0000000e+00,0.0000000e+00
 1,,temperature,nu-fission,0.0000000e+00,0.0000000e+00
-1,,temperature,scatter,-2.0884127e-06,1.6684645e-06
-1,,temperature,nu-fission,-1.5670247e-05,4.5939944e-05
-1,,temperature,scatter,3.9749436e-06,3.9749436e-06
+1,,temperature,scatter,-2.9956385e-07,2.8883252e-07
+1,,temperature,nu-fission,-2.1023382e-05,4.9883940e-05
+1,,temperature,scatter,5.9452821e-09,8.7655524e-09
 1,,temperature,nu-fission,0.0000000e+00,0.0000000e+00
 1,,temperature,scatter,0.0000000e+00,0.0000000e+00
 1,,temperature,nu-fission,0.0000000e+00,0.0000000e+00

tests/regression_tests/multipole/results_true.dat

@@ -1,36 +1,36 @@
 k-combined:
-1.363786E+00 1.103929E-02
+1.342579E+00 1.221176E-02
 tally 1:
-3.960375E+00
-3.138356E+00
-2.875799E+00
-1.655329E+00
-5.521904E-01
-6.105083E-02
-4.617974E-01
-4.274130E-02
+3.839794E+00
+2.951721E+00
+2.785273E+00
+1.554128E+00
+5.349716E-01
+5.732174E-02
+4.499834E-01
+4.055011E-02
 0.000000E+00
 0.000000E+00
-2.254165E+01
-1.016444E+02
+2.258507E+01
+1.020294E+02
 0.000000E+00
 0.000000E+00
-6.839351E-04
-9.359599E-08
-2.251829E+01
-1.014341E+02
-4.903575E-05
-1.199780E-09
-3.595002E+02
-2.586079E+04
-2.875799E+00
-1.655329E+00
-2.174747E+00
-9.465589E-01
-3.543564E+02
-2.512619E+04
-4.903575E-05
-1.199780E-09
+6.862242E-04
+9.421377E-08
+2.256396E+01
+1.018388E+02
+1.146136E-04
+3.296980E-09
+3.624627E+02
+2.628739E+04
+2.785273E+00
+1.554128E+00
+2.176478E+00
+9.480405E-01
+3.574110E+02
+2.555993E+04
+1.146136E-04
+3.296980E-09
 Cell
 	ID             =	11
 	Name           =	

tests/unit_tests/test_data_multipole.py

@@ -3,8 +3,6 @@ import os
 import numpy as np
 import pytest
 import openmc.data
-
-
 pytestmark = pytest.mark.skipif(
     'OPENMC_MULTIPOLE_LIBRARY' not in os.environ,
     reason='OPENMC_MULTIPOLE_LIBRARY environment variable must be set')
@@ -18,44 +16,32 @@ def u235():
 
 
 @pytest.fixture(scope='module')
-def u234():
+def b10():
     directory = os.environ['OPENMC_MULTIPOLE_LIBRARY']
-    filename = os.path.join(directory, '092234.h5')
+    filename = os.path.join(directory, '005010.h5')
     return openmc.data.WindowedMultipole.from_hdf5(filename)
 
 
-@pytest.fixture(scope='module')
-def fe56():
-    directory = os.environ['OPENMC_MULTIPOLE_LIBRARY']
-    filename = os.path.join(directory, '026056.h5')
-    return openmc.data.WindowedMultipole.from_hdf5(filename)
-
-
-def test_evaluate_rm(u235):
-    """Make sure a Reich-Moore multipole object can be called."""
+def test_evaluate(u235):
+    """Test the cross section evaluation of a library."""
     energies = [1e-3, 1.0, 10.0, 50.]
-    total, absorption, fission = u235(energies, 0.0)
-    assert total[1] == pytest.approx(90.64895383)
-    total, absorption, fission = u235(energies, 300.0)
-    assert total[1] == pytest.approx(91.12534964)
+    scattering, absorption, fission = u235(energies, 0.0)
+    assert (scattering[1], absorption[1], fission[1]) == \
+           pytest.approx((13.09, 77.56, 67.36), rel=1e-3)
+    scattering, absorption, fission = u235(energies, 300.0)
+    assert (scattering[2], absorption[2], fission[2]) == \
+           pytest.approx((11.24, 21.26, 15.50), rel=1e-3)
 
 
-def test_evaluate_mlbw(u234):
-    """Make sure a Multi-Level Breit-Wigner multipole object can be called."""
-    energies = [1e-3, 1.0, 10.0, 50.]
-    total, absorption, fission = u234(energies, 0.0)
-    assert total[3] == pytest.approx(15.02827953)
-    total, absorption, fission = u234(energies, 300.0)
-    assert total[3] == pytest.approx(15.08269143)
-
-
-def test_high_l(fe56):
-    """Test a nuclide (Fe56) with a high l-value (4)."""
+def test_evaluate_none_poles(b10):
+    """Test a library with no poles, i.e., purely polynomials."""
     energies = [1e-3, 1.0, 10.0, 1e3, 1e5]
-    total, absorption, fission = fe56(energies, 0.0)
-    assert total[0] == pytest.approx(25.072619556789267)
-    total, absorption, fission = fe56(energies, 300.0)
-    assert total[0] == pytest.approx(27.85535792368082)
+    scattering, absorption, fission = b10(energies, 0.0)
+    assert (scattering[0], absorption[0], fission[0]) == \
+           pytest.approx((2.201, 19330., 0.), rel=1e-3)
+    scattering, absorption, fission = b10(energies, 300.0)
+    assert (scattering[-1], absorption[-1], fission[-1]) == \
+           pytest.approx((2.878, 1.982, 0.), rel=1e-3)
 
 
 def test_export_to_hdf5(tmpdir, u235):

tools/ci/travis-before-script.sh

@@ -17,5 +17,7 @@ if [[ ! -d $ENDF/neutrons || ! -d $ENDF/photoat || ! -d $ENDF/atomic_relax ]]; t
 fi
 
 # Download multipole library
-git clone --branch=master git://github.com/smharper/windowed_multipole_library.git wmp_lib
-tar -C $HOME -xzvf wmp_lib/multipole_lib.tar.gz
+if [[ ! -e $HOME/WMP_Library/092235.h5 ]]; then
+    wget https://github.com/mit-crpg/WMP_Library/releases/download/v1.0/WMP_Library_v1.0.tar.gz
+    tar -C $HOME -xzvf WMP_Library_v1.0.tar.gz
+fi