From 8e10d9e00fefb104692a9fc21978d8ae467b5cb2 Mon Sep 17 00:00:00 2001
From: Adam Nelson <nelsonag@umich.edu>
Date: Mon, 17 Oct 2016 19:15:11 -0400
Subject: [PATCH 01/12] seem to have the input side for the MuFilter working,
 now on to the library generation side

---
 openmc/filter.py       | 313 +++++++++++++++++++++++++++++++----------
 openmc/mgxs/library.py |  82 +++++++++--
 openmc/mgxs/mgxs.py    | 210 ++++++++++++++++++---------
 openmc/mgxs_library.py |  39 +++--
 openmc/tallies.py      |   3 +-
 5 files changed, 474 insertions(+), 173 deletions(-)

diff --git a/openmc/filter.py b/openmc/filter.py
index d11af96f5..91ceb2931 100644
--- a/openmc/filter.py
+++ b/openmc/filter.py
@@ -105,7 +105,8 @@ class Filter(with_metaclass(FilterMeta, object)):
         """Return all subclasses and their subclasses, etc."""
         subs = cls.__subclasses__()
         subsubs = [grand for s in subs for grand in s.__subclasses__()]
-        return subs + subsubs
+        subsubsubs = [grand for s in subsubs for grand in s.__subclasses__()]
+        return subs + subsubs + subsubsubs
 
     @classmethod
     def from_hdf5(cls, group, **kwargs):
@@ -777,7 +778,116 @@ class MeshFilter(Filter):
         return df
 
 
-class EnergyFilter(Filter):
+class RealFilter(Filter):
+    """Tally modifier that describes phase-space and other characteristics
+
+    Parameters
+    ----------
+    bins : Iterable of Real
+        A grid of bin values.
+
+    Attributes
+    ----------
+    bins : Iterable of Real
+        A grid of bin values.
+    num_bins : Integral
+        The number of filter bins
+    stride : Integral
+        The number of filter, nuclide and score bins within each of this
+        filter's bins.
+
+    """
+
+    def __gt__(self, other):
+        if type(self) is type(other):
+            # Compare largest/smallest bin edges in filters
+            # This logic is used when merging tallies with real filters
+            return self.bins[0] >= other.bins[-1]
+        else:
+            return super(RealFilter, self).__gt__(other)
+
+    @property
+    def num_bins(self):
+        return len(self.bins) - 1
+
+    @num_bins.setter
+    def num_bins(self, num_bins):
+        cv.check_type('filter num_bins', num_bins, Integral)
+        cv.check_greater_than('filter num_bins', num_bins, 0, equality=True)
+        self._num_bins = num_bins
+
+    def can_merge(self, other):
+        if type(self) is not type(other):
+            return False
+
+        if self.bins[0] == other.bins[-1]:
+            # This low energy edge coincides with other's high edge
+            return True
+        elif self.bins[-1] == other.bins[0]:
+            # This high energy edge coincides with other's low edge
+            return True
+        else:
+            return False
+
+    def merge(self, other):
+        if not self.can_merge(other):
+            msg = 'Unable to merge "{0}" with "{1}" ' \
+                  'filters'.format(self.type, other.type)
+            raise ValueError(msg)
+
+        # Merge unique filter bins
+        merged_bins = np.concatenate((self.bins, other.bins))
+        merged_bins = np.unique(merged_bins)
+
+        # Create a new filter with these bins
+        return type(self)(sorted(merged_bins))
+
+    def is_subset(self, other):
+        """Determine if another filter is a subset of this filter.
+
+        If all of the bins in the other filter are included as bins in this
+        filter, then it is a subset of this filter.
+
+        Parameters
+        ----------
+        other : openmc.Filter
+            The filter to query as a subset of this filter
+
+        Returns
+        -------
+        bool
+            Whether or not the other filter is a subset of this filter
+
+        """
+
+        if type(self) is not type(other):
+            return False
+        elif len(self.bins) != len(other.bins):
+            return False
+        else:
+            return np.allclose(self.bins, other.bins)
+
+    def get_bin_index(self, filter_bin):
+        # Use lower energy bound to find index for RealFilters
+        deltas = np.abs(self.bins - filter_bin[1]) / filter_bin[1]
+        min_delta = np.min(deltas)
+        if min_delta < 1E-3:
+            return deltas.argmin() - 1
+        else:
+            msg = 'Unable to get the bin index for Filter since "{0}" ' \
+                  'is not one of the bins'.format(filter_bin)
+            raise ValueError(msg)
+
+    def get_bin(self, bin_index):
+        cv.check_type('bin_index', bin_index, Integral)
+        cv.check_greater_than('bin_index', bin_index, 0, equality=True)
+        cv.check_less_than('bin_index', bin_index, self.num_bins)
+
+        # Construct 2-tuple of lower, upper bins for real-valued filters
+        return (self.bins[bin_index], self.bins[bin_index + 1])
+
+
+class EnergyFilter(RealFilter):
     """Bins tally events based on incident particle energy.
 
     Parameters
@@ -797,24 +907,6 @@ class EnergyFilter(Filter):
 
     """
 
-    def __gt__(self, other):
-        if type(self) is type(other):
-            # Compare largest/smallest energy bin edges in energy filters
-            # This logic is used when merging tallies with energy filters
-            return self.bins[0] >= other.bins[-1]
-        else:
-            return super(EnergyFilter, self).__gt__(other)
-
-    @property
-    def num_bins(self):
-        return len(self.bins) - 1
-
-    @num_bins.setter
-    def num_bins(self, num_bins):
-        cv.check_type('filter num_bins', num_bins, Integral)
-        cv.check_greater_than('filter num_bins', num_bins, 0, equality=True)
-        self._num_bins = num_bins
-
     def check_bins(self, bins):
         for edge in bins:
             if not isinstance(edge, Real):
@@ -835,59 +927,6 @@ class EnergyFilter(Filter):
                       'increasing'.format(bins, self.type)
                 raise ValueError(msg)
 
-    def can_merge(self, other):
-        if type(self) is not type(other):
-            return False
-
-        if self.bins[0] == other.bins[-1]:
-            # This low energy edge coincides with other's high energy edge
-            return True
-        elif self.bins[-1] == other.bins[0]:
-            # This high energy edge coincides with other's low energy edge
-            return True
-        else:
-            return False
-
-    def merge(self, other):
-        if not self.can_merge(other):
-            msg = 'Unable to merge "{0}" with "{1}" ' \
-                  'filters'.format(self.type, other.type)
-            raise ValueError(msg)
-
-        # Merge unique filter bins
-        merged_bins = np.concatenate((self.bins, other.bins))
-        merged_bins = np.unique(merged_bins)
-
-        # Create a new filter with these bins
-        return type(self)(sorted(merged_bins))
-
-    def is_subset(self, other):
-        if type(self) is not type(other):
-            return False
-        elif len(self.bins) != len(other.bins):
-            return False
-        else:
-            return np.allclose(self.bins, other.bins)
-
-    def get_bin_index(self, filter_bin):
-        # Use lower energy bound to find index for energy Filters
-        deltas = np.abs(self.bins - filter_bin[1]) / filter_bin[1]
-        min_delta = np.min(deltas)
-        if min_delta < 1E-3:
-            return deltas.argmin() - 1
-        else:
-            msg = 'Unable to get the bin index for Filter since "{0}" ' \
-                  'is not one of the bins'.format(filter_bin)
-            raise ValueError(msg)
-
-    def get_bin(self, bin_index):
-        cv.check_type('bin_index', bin_index, Integral)
-        cv.check_greater_than('bin_index', bin_index, 0, equality=True)
-        cv.check_less_than('bin_index', bin_index, self.num_bins)
-
-        # Construct 2-tuple of lower, upper energies for energy(out) filters
-        return (self.bins[bin_index], self.bins[bin_index+1])
-
     def get_pandas_dataframe(self, data_size, **kwargs):
         """Builds a Pandas DataFrame for the Filter's bins.
 
@@ -1231,7 +1270,7 @@ class DistribcellFilter(Filter):
         return df
 
 
-class MuFilter(Filter):
+class MuFilter(RealFilter):
     """Bins tally events based on particle scattering angle.
 
     Parameters
@@ -1259,8 +1298,82 @@ class MuFilter(Filter):
 
     """
 
+    def check_bins(self, bins):
+        for edge in bins:
+            if not isinstance(edge, Real):
+                msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
+                      'since it is a non-integer or floating point ' \
+                      'value'.format(edge, type(self))
+                raise ValueError(msg)
+            elif edge < -1.:
+                msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
+                      'since it is less than -1'.format(edge, type(self))
+                raise ValueError(msg)
+            elif edge > 1.:
+                msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
+                      'since it is greater than 1'.format(edge, type(self))
+                raise ValueError(msg)
 
-class PolarFilter(Filter):
+        # Check that bin edges are monotonically increasing
+        for index in range(1, len(bins)):
+            if bins[index] < bins[index-1]:
+                msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \
+                      'since they are not monotonically ' \
+                      'increasing'.format(bins, self.type)
+                raise ValueError(msg)
+
+    def get_pandas_dataframe(self, data_size, **kwargs):
+        """Builds a Pandas DataFrame for the Filter's bins.
+
+        This method constructs a Pandas DataFrame object for the filter with
+        columns annotated by filter bin information. This is a helper method
+        for :meth:`Tally.get_pandas_dataframe`.
+
+        Parameters
+        ----------
+        data_size : Integral
+            The total number of bins in the tally corresponding to this filter
+
+        Returns
+        -------
+        pandas.DataFrame
+            A Pandas DataFrame with one column of the lower energy bound and one
+            column of upper energy bound for each filter bin.  The number of
+            rows in the DataFrame is the same as the total number of bins in the
+            corresponding tally, with the filter bin appropriately tiled to map
+            to the corresponding tally bins.
+
+        Raises
+        ------
+        ImportError
+            When Pandas is not installed
+
+        See also
+        --------
+        Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe()
+
+        """
+
+        # Initialize Pandas DataFrame
+        import pandas as pd
+        df = pd.DataFrame()
+
+        # Extract the lower and upper energy bounds, then repeat and tile
+        # them as necessary to account for other filters.
+        lo_bins = np.repeat(self.bins[:-1], self.stride)
+        hi_bins = np.repeat(self.bins[1:], self.stride)
+        tile_factor = data_size / len(lo_bins)
+        lo_bins = np.tile(lo_bins, tile_factor)
+        hi_bins = np.tile(hi_bins, tile_factor)
+
+        # Add the new energy columns to the DataFrame.
+        df.loc[:, self.short_name.lower() + ' low'] = lo_bins
+        df.loc[:, self.short_name.lower() + ' high'] = hi_bins
+
+        return df
+
+
+class PolarFilter(RealFilter):
     """Bins tally events based on the incident particle's direction.
 
     Parameters
@@ -1288,6 +1401,30 @@ class PolarFilter(Filter):
 
     """
 
+    def check_bins(self, bins):
+        for edge in bins:
+            if not isinstance(edge, Real):
+                msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
+                      'since it is a non-integer or floating point ' \
+                      'value'.format(edge, type(self))
+                raise ValueError(msg)
+            elif edge < 0.:
+                msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
+                      'since it is less than 0'.format(edge, type(self))
+                raise ValueError(msg)
+            elif edge > np.pi:
+                msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
+                      'since it is greater than pi'.format(edge, type(self))
+                raise ValueError(msg)
+
+        # Check that bin edges are monotonically increasing
+        for index in range(1, len(bins)):
+            if bins[index] < bins[index-1]:
+                msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \
+                      'since they are not monotonically ' \
+                      'increasing'.format(bins, self.type)
+                raise ValueError(msg)
+
     def get_pandas_dataframe(self, data_size, **kwargs):
         """Builds a Pandas DataFrame for the Filter's bins.
 
@@ -1338,7 +1475,7 @@ class PolarFilter(Filter):
         return df
 
 
-class AzimuthalFilter(Filter):
+class AzimuthalFilter(RealFilter):
     """Bins tally events based on the incident particle's direction.
 
     Parameters
@@ -1366,6 +1503,30 @@ class AzimuthalFilter(Filter):
 
     """
 
+    def check_bins(self, bins):
+        for edge in bins:
+            if not isinstance(edge, Real):
+                msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
+                      'since it is a non-integer or floating point ' \
+                      'value'.format(edge, type(self))
+                raise ValueError(msg)
+            elif edge < -np.pi:
+                msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
+                      'since it is less than -pi'.format(edge, type(self))
+                raise ValueError(msg)
+            elif edge > np.pi:
+                msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
+                      'since it is greater than pi'.format(edge, type(self))
+                raise ValueError(msg)
+
+        # Check that bin edges are monotonically increasing
+        for index in range(1, len(bins)):
+            if bins[index] < bins[index-1]:
+                msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \
+                      'since they are not monotonically ' \
+                      'increasing'.format(bins, self.type)
+                raise ValueError(msg)
+
     def get_pandas_dataframe(self, data_size, distribcell_paths=True):
         """Builds a Pandas DataFrame for the Filter's bins.
 
diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py
index c03aa7877..96bfa6fab 100644
--- a/openmc/mgxs/library.py
+++ b/openmc/mgxs/library.py
@@ -62,15 +62,23 @@ class Library(object):
         The spatial domain(s) for which MGXS in the Library are computed
     correction : {'P0', None}
         Apply the P0 correction to scattering matrices if set to 'P0'
+    scatter_format : {'legendre', or 'histogram'}
+        Representation of the angular scattering distribution (default is
+        'legendre')
     legendre_order : int
-        The highest legendre moment in the scattering matrices (default is 0)
+        The highest Legendre moment in the scattering matrix; this is used if
+        :attr:`ScatterMatrixXS.scatter_format` is 'legendre'. (default is 0)
+    histogram_bins : int
+        The number of equally-spaced bins for the histogram representation of
+        the angular scattering distribution; this is used if
+        :attr:`ScatterMatrixXS.scatter_format` is 'histogram'. (default is 16)
     energy_groups : openmc.mgxs.EnergyGroups
         Energy group structure for energy condensation
     delayed_groups : list of int
         Delayed groups to filter out the xs
     estimator : str or None
-        The tally estimator used to compute multi-group cross sections. If None,
-        the default for each MGXS type is used.
+        The tally estimator used to compute multi-group cross sections.
+        If None, the default for each MGXS type is used.
     tally_trigger : openmc.Trigger
         An (optional) tally precision trigger given to each tally used to
         compute the cross section
@@ -104,7 +112,9 @@ class Library(object):
         self._energy_groups = None
         self._delayed_groups = None
         self._correction = 'P0'
+        self._scatter_format = 'legendre'
         self._legendre_order = 0
+        self._histogram_bins = 16
         self._tally_trigger = None
         self._all_mgxs = OrderedDict()
         self._sp_filename = None
@@ -133,7 +143,9 @@ class Library(object):
             clone._domain_type = self.domain_type
             clone._domains = copy.deepcopy(self.domains)
             clone._correction = self.correction
+            clone._scatter_format = self.scatter_format
             clone._legendre_order = self.legendre_order
+            clone._histogram_bins = self.histogram_bins
             clone._energy_groups = copy.deepcopy(self.energy_groups, memo)
             clone._delayed_groups = copy.deepcopy(self.delayed_groups, memo)
             clone._tally_trigger = copy.deepcopy(self.tally_trigger, memo)
@@ -212,10 +224,18 @@ class Library(object):
     def correction(self):
         return self._correction
 
+    @property
+    def scatter_format(self):
+        return self._scatter_format
+
     @property
     def legendre_order(self):
         return self._legendre_order
 
+    @property
+    def histogram_bins(self):
+        return self._histogram_bins
+
     @property
     def tally_trigger(self):
         return self._tally_trigger
@@ -355,26 +375,57 @@ class Library(object):
     def correction(self, correction):
         cv.check_value('correction', correction, ('P0', None))
 
-        if correction == 'P0' and self.legendre_order > 0:
-            warn('The P0 correction will be ignored since the scattering '
-                 'order "{}" is greater than zero'.format(self.legendre_order))
+        if self.scatter_format == 'legendre':
+            if correction == 'P0' and self.legendre_order > 0:
+                msg = 'The P0 correction will be ignored since the ' \
+                      'scattering order {} is greater than '\
+                      'zero'.format(self.legendre_order)
+                warn(msg)
+        elif self.scatter_format == 'histogram':
+            msg = 'The P0 correction will be ignored since the ' \
+                  'scatter format is set to histogram'
+            warn(msg)
 
         self._correction = correction
 
+    @scatter_format.setter
+    def scatter_format(self, scatter_format):
+        cv.check_value('scatter_format', scatter_format, openmc.mgxs.MU_TREATMENTS)
+        self._scatter_format = scatter_format
+
     @legendre_order.setter
     def legendre_order(self, legendre_order):
         cv.check_type('legendre_order', legendre_order, Integral)
-        cv.check_greater_than('legendre_order', legendre_order, 0, equality=True)
+        cv.check_greater_than('legendre_order', legendre_order, 0,
+                              equality=True)
         cv.check_less_than('legendre_order', legendre_order, 10, equality=True)
 
-        if self.correction == 'P0' and legendre_order > 0:
-            msg = 'The P0 correction will be ignored since the scattering ' \
-                  'order {} is greater than zero'.format(self.legendre_order)
-            warn(msg, RuntimeWarning)
-            self.correction = None
+        if self.scatter_format == 'legendre':
+            if self.correction == 'P0' and legendre_order > 0:
+                msg = 'The P0 correction will be ignored since the ' \
+                      'scattering order {} is greater than '\
+                      'zero'.format(self.legendre_order)
+                warn(msg, RuntimeWarning)
+                self.correction = None
+        elif self.scatter_format == 'histogram':
+            msg = 'The legendre order will be ignored since the ' \
+                  'scatter format is set to histogram'
+            warn(msg)
 
         self._legendre_order = legendre_order
 
+    @histogram_bins.setter
+    def histogram_bins(self, histogram_bins):
+        cv.check_type('histogram_bins', histogram_bins, Integral)
+        cv.check_greater_than('histogram_bins', histogram_bins, 0)
+
+        if self.scatter_format == 'legendre':
+            msg = 'The histogram bins will be ignored since the ' \
+                  'scatter format is set to legendre'
+            warn(msg)
+
+        self._histogram_bins = histogram_bins
+
     @tally_trigger.setter
     def tally_trigger(self, tally_trigger):
         cv.check_type('tally trigger', tally_trigger, openmc.Trigger)
@@ -443,7 +494,9 @@ class Library(object):
                 # Specify whether to use a transport ('P0') correction
                 if isinstance(mgxs, openmc.mgxs.ScatterMatrixXS):
                     mgxs.correction = self.correction
+                    mgxs.scatter_format = self.scatter_format
                     mgxs.legendre_order = self.legendre_order
+                    mgxs.histogram_bins = self.histogram_bins
 
                 self.all_mgxs[domain.id][mgxs_type] = mgxs
 
@@ -901,15 +954,14 @@ class Library(object):
         xsdata = openmc.XSdata(name, self.energy_groups)
 
         if order is None:
-            # Set the order to the Library's order (the defualt behavior)
+            # Set the order to the Library's order (the default behavior)
             xsdata.order = self.legendre_order
         else:
             # Set the order of the xsdata object to the minimum of
             # the provided order or the Library's order.
             xsdata.order = min(order, self.legendre_order)
 
-        # Right now only 'legendre' data and isotropic weighting is supported
-        self.scatter_format = 'legendre'
+        # Right now only isotropic weighting is supported
         self.representation = 'isotropic'
 
         if nuclide != 'total':
diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py
index 65802f2ef..cafda1ecf 100644
--- a/openmc/mgxs/mgxs.py
+++ b/openmc/mgxs/mgxs.py
@@ -60,6 +60,9 @@ _DOMAINS = (openmc.Cell,
             openmc.Material,
             openmc.Mesh)
 
+# Supported ScatterMatrixXS and NuScatterMatrixXS angular distribution types
+MU_TREATMENTS = ('legendre', 'histogram')
+
 
 class MGXS(object):
     """An abstract multi-group cross section for some energy group structure
@@ -3212,8 +3215,8 @@ class NuScatterXS(MGXS):
 
 
 class ScatterMatrixXS(MatrixMGXS):
-    r"""A scattering matrix multi-group cross section for one or more Legendre
-    moments.
+    r"""A scattering matrix multi-group cross section with the cosine of the
+    change-in-angle represented as one or more Legendre moments or a histogram.
 
     This class can be used for both OpenMC input generation and tally data
     post-processing to compute spatially-homogenized and energy-integrated
@@ -3231,7 +3234,7 @@ class ScatterMatrixXS(MatrixMGXS):
 
     For a spatial domain :math:`V`, incoming energy group
     :math:`[E_{g'},E_{g'-1}]`, and outgoing energy group :math:`[E_g,E_{g-1}]`,
-    the scattering moments are calculated as:
+    the Legendre scattering moments are calculated as:
 
     .. math::
 
@@ -3271,9 +3274,18 @@ class ScatterMatrixXS(MatrixMGXS):
     Attributes
     ----------
     correction : 'P0' or None
-        Apply the P0 correction to scattering matrices if set to 'P0'
+        Apply the P0 correction to scattering matrices if set to 'P0'; this is
+        used only if :attr:`ScatterMatrixXS.scatter_format` is 'legendre'
+    scatter_format : {'legendre', or 'histogram'}
+        Representation of the angular scattering distribution (default is
+        'legendre')
     legendre_order : int
-        The highest Legendre moment in the scattering matrix (default is 0)
+        The highest Legendre moment in the scattering matrix; this is used if
+        :attr:`ScatterMatrixXS.scatter_format` is 'legendre'. (default is 0)
+    histogram_bins : int
+        The number of equally-spaced bins for the histogram representation of
+        the angular scattering distribution; this is used if
+        :attr:`ScatterMatrixXS.scatter_format` is 'histogram'. (default is 16)
     name : str, optional
         Name of the multi-group cross section
     rxn_type : str
@@ -3340,7 +3352,9 @@ class ScatterMatrixXS(MatrixMGXS):
                                               groups, by_nuclide, name)
         self._rxn_type = 'scatter'
         self._correction = 'P0'
+        self._scatter_format = 'legendre'
         self._legendre_order = 0
+        self._histogram_bins = 16
         self._hdf5_key = 'scatter matrix'
         self._estimator = 'analog'
         self._valid_estimators = ['analog']
@@ -3348,26 +3362,39 @@ class ScatterMatrixXS(MatrixMGXS):
     def __deepcopy__(self, memo):
         clone = super(ScatterMatrixXS, self).__deepcopy__(memo)
         clone._correction = self.correction
+        clone._scatter_format = self.scatter_format
         clone._legendre_order = self.legendre_order
+        clone._histogram_bins = self.histogram_bins
         return clone
 
     @property
     def correction(self):
         return self._correction
 
+    @property
+    def scatter_format(self):
+        return self._scatter_format
+
     @property
     def legendre_order(self):
         return self._legendre_order
 
+    @property
+    def histogram_bins(self):
+        return self._histogram_bins
+
     @property
     def scores(self):
         scores = ['flux']
 
-        if self.correction == 'P0' and self.legendre_order == 0:
-            scores += ['{}-0'.format(self.rxn_type),
-                       '{}-1'.format(self.rxn_type)]
-        else:
-            scores += ['{}-P{}'.format(self.rxn_type, self.legendre_order)]
+        if self.scatter_format == 'legendre':
+            if self.correction == 'P0' and self.legendre_order == 0:
+                scores += ['{}-0'.format(self.rxn_type),
+                           '{}-1'.format(self.rxn_type)]
+            else:
+                scores += ['{}-P{}'.format(self.rxn_type, self.legendre_order)]
+        elif self.scatter_format == 'histogram':
+            scores += [self.rxn_type]
 
         return scores
 
@@ -3377,10 +3404,14 @@ class ScatterMatrixXS(MatrixMGXS):
         energy = openmc.EnergyFilter(group_edges)
         energyout = openmc.EnergyoutFilter(group_edges)
 
-        if self.correction == 'P0' and self.legendre_order == 0:
-            filters = [[energy], [energy, energyout], [energyout]]
-        else:
-            filters = [[energy], [energy, energyout]]
+        if self.scatter_format == 'legendre':
+            if self.correction == 'P0' and self.legendre_order == 0:
+                filters = [[energy], [energy, energyout], [energyout]]
+            else:
+                filters = [[energy], [energy, energyout]]
+        elif self.scatter_format == 'histogram':
+            bins = np.linspace(-1., 1., num=self.histogram_bins, endpoint=True)
+            filters = [[energy], [energy, energyout, openmc.MuFilter(bins)]]
 
         return filters
 
@@ -3388,20 +3419,24 @@ class ScatterMatrixXS(MatrixMGXS):
     def rxn_rate_tally(self):
 
         if self._rxn_rate_tally is None:
+            if self.scatter_format == 'legendre':
+                # If using P0 correction subtract scatter-1 from the diagonal
+                if self.correction == 'P0' and self.legendre_order == 0:
+                    scatter_p0 = self.tallies['{}-0'.format(self.rxn_type)]
+                    scatter_p1 = self.tallies['{}-1'.format(self.rxn_type)]
+                    energy_filter = scatter_p0.find_filter(openmc.EnergyFilter)
+                    energy_filter = copy.deepcopy(energy_filter)
+                    scatter_p1 = scatter_p1.diagonalize_filter(energy_filter)
+                    self._rxn_rate_tally = scatter_p0 - scatter_p1
 
-            # If using P0 correction subtract scatter-1 from the diagonal
-            if self.correction == 'P0' and self.legendre_order == 0:
-                scatter_p0 = self.tallies['{}-0'.format(self.rxn_type)]
-                scatter_p1 = self.tallies['{}-1'.format(self.rxn_type)]
-                energy_filter = scatter_p0.find_filter(openmc.EnergyFilter)
-                energy_filter = copy.deepcopy(energy_filter)
-                scatter_p1 = scatter_p1.diagonalize_filter(energy_filter)
-                self._rxn_rate_tally = scatter_p0 - scatter_p1
-
-            # Extract scattering moment reaction rate Tally
-            else:
-                tally_key = '{}-P{}'.format(self.rxn_type, self.legendre_order)
-                self._rxn_rate_tally = self.tallies[tally_key]
+                # Extract scattering moment reaction rate Tally
+                else:
+                    tally_key = '{}-P{}'.format(self.rxn_type,
+                                                self.legendre_order)
+                    self._rxn_rate_tally = self.tallies[tally_key]
+            elif self.scatter_format == 'histogram':
+                # Extract scattering rate distribution tally
+                self._rxn_rate_tally = self.tallies[self.rxn_type]
 
             self._rxn_rate_tally.sparse = self.sparse
 
@@ -3411,27 +3446,52 @@ class ScatterMatrixXS(MatrixMGXS):
     def correction(self, correction):
         cv.check_value('correction', correction, ('P0', None))
 
-        if correction == 'P0' and self.legendre_order > 0:
-            msg = 'The P0 correction will be ignored since the scattering ' \
-                  'order {} is greater than zero'.format(self.legendre_order)
+        if self.scatter_format == 'legendre':
+            if correction == 'P0' and self.legendre_order > 0:
+                msg = 'The P0 correction will be ignored since the ' \
+                      'scattering order {} is greater than '\
+                      'zero'.format(self.legendre_order)
+                warnings.warn(msg)
+        elif self.scatter_format == 'histogram':
+            msg = 'The P0 correction will be ignored since the ' \
+                  'scatter format is set to histogram'
             warnings.warn(msg)
 
         self._correction = correction
 
+    @scatter_format.setter
+    def scatter_format(self, scatter_format):
+        cv.check_value('scatter_format', scatter_format, MU_TREATMENTS)
+        self._scatter_format = scatter_format
+
     @legendre_order.setter
     def legendre_order(self, legendre_order):
         cv.check_type('legendre_order', legendre_order, Integral)
-        cv.check_greater_than('legendre_order', legendre_order, 0, equality=True)
+        cv.check_greater_than('legendre_order', legendre_order, 0,
+                              equality=True)
         cv.check_less_than('legendre_order', legendre_order, 10, equality=True)
 
-        if self.correction == 'P0' and legendre_order > 0:
-            msg = 'The P0 correction will be ignored since the scattering ' \
-                  'order {} is greater than zero'.format(self.legendre_order)
-            warnings.warn(msg, RuntimeWarning)
-            self.correction = None
+        if self.scatter_format == 'legendre':
+            if self.correction == 'P0' and legendre_order > 0:
+                msg = 'The P0 correction will be ignored since the ' \
+                      'scattering order {} is greater than '\
+                      'zero'.format(self.legendre_order)
+                warnings.warn(msg, RuntimeWarning)
+                self.correction = None
+        elif self.scatter_format == 'histogram':
+            msg = 'The legendre order will be ignored since the ' \
+                  'scatter format is set to histogram'
+            warnings.warn(msg)
 
         self._legendre_order = legendre_order
 
+    @histogram_bins.setter
+    def histogram_bins(self, histogram_bins):
+        cv.check_type('histogram_bins', histogram_bins, Integral)
+        cv.check_greater_than('histogram_bins', histogram_bins, 0)
+
+        self._histogram_bins = histogram_bins
+
     def load_from_statepoint(self, statepoint):
         """Extracts tallies in an OpenMC StatePoint with the data needed to
         compute multi-group cross sections.
@@ -3461,12 +3521,16 @@ class ScatterMatrixXS(MatrixMGXS):
             self._rxn_rate_tally = None
             self._loaded_sp = False
 
-        # Expand scores to match the format in the statepoint
-        # e.g., "scatter-P2" -> "scatter-0", "scatter-1", "scatter-2"
-        if self.correction != 'P0' or self.legendre_order != 0:
-            tally_key = '{}-P{}'.format(self.rxn_type, self.legendre_order)
-            self.tallies[tally_key].scores = \
-                [self.rxn_type + '-{}'.format(i) for i in range(self.legendre_order+1)]
+        if self.scatter_format == 'legendre':
+            # Expand scores to match the format in the statepoint
+            # e.g., "scatter-P2" -> "scatter-0", "scatter-1", "scatter-2"
+            if self.correction != 'P0' or self.legendre_order != 0:
+                tally_key = '{}-P{}'.format(self.rxn_type, self.legendre_order)
+                self.tallies[tally_key].scores = \
+                    [self.rxn_type + '-{}'.format(i)
+                     for i in range(self.legendre_order + 1)]
+        elif self.scatter_format == 'histogram':
+            self.tallies[self.rxn_type].scores = [self.rxn_type]
 
         super(ScatterMatrixXS, self).load_from_statepoint(statepoint)
 
@@ -3513,7 +3577,7 @@ class ScatterMatrixXS(MatrixMGXS):
         slice_xs._xs_tally = None
 
         # Slice the Legendre order if needed
-        if legendre_order != 'same':
+        if legendre_order != 'same' and self.scatter_format == 'legendre':
             cv.check_type('legendre_order', legendre_order, Integral)
             cv.check_less_than('legendre_order', legendre_order,
                                self.legendre_order, equality=True)
@@ -3522,7 +3586,8 @@ class ScatterMatrixXS(MatrixMGXS):
             # Slice the scattering tally
             tally_key = '{}-P{}'.format(self.rxn_type, self.legendre_order)
             expand_scores = \
-                [self.rxn_type + '-{}'.format(i) for i in range(self.legendre_order+1)]
+                [self.rxn_type + '-{}'.format(i)
+                 for i in range(self.legendre_order + 1)]
             slice_xs.tallies[tally_key] = \
                 slice_xs.tallies[tally_key].get_slice(scores=expand_scores)
 
@@ -3553,7 +3618,8 @@ class ScatterMatrixXS(MatrixMGXS):
         This method constructs a 5D NumPy array for the requested
         multi-group cross section data for one or more subdomains
         (1st dimension), energy groups in (2nd dimension), energy groups out
-        (3rd dimension), nuclides (4th dimension), and moments (5th dimension).
+        (3rd dimension), nuclides (4th dimension), and moments/histograms
+        (5th dimension).
 
         NOTE: The scattering moments are not multiplied by the :math:`(2l+1)/2`
         prefactor in the expansion of the scattering source into Legendre
@@ -3642,7 +3708,7 @@ class ScatterMatrixXS(MatrixMGXS):
                 filter_bins.append((self.energy_groups.get_group_bounds(group),))
 
         # Construct CrossScore for requested scattering moment
-        if moment != 'all':
+        if moment != 'all' and self.scatter_format == 'legendre':
             cv.check_type('moment', moment, Integral)
             cv.check_greater_than('moment', moment, 0, equality=True)
             cv.check_less_than(
@@ -3760,28 +3826,38 @@ class ScatterMatrixXS(MatrixMGXS):
         df = super(ScatterMatrixXS, self).get_pandas_dataframe(
             groups, nuclides, xs_type, distribcell_paths)
 
-        # Add a moment column to dataframe
-        if self.legendre_order > 0:
-            # Insert a column corresponding to the Legendre moments
-            moments = ['P{}'.format(i) for i in range(self.legendre_order+1)]
-            moments = np.tile(moments, int(df.shape[0] / len(moments)))
-            df['moment'] = moments
+        if self.scatter_format == 'legendre':
+            # Add a moment column to dataframe
+            if self.legendre_order > 0:
+                # Insert a column corresponding to the Legendre moments
+                moments = ['P{}'.format(i)
+                           for i in range(self.legendre_order + 1)]
+                moments = np.tile(moments, int(df.shape[0] / len(moments)))
+                df['moment'] = moments
 
-            # Place the moment column before the mean column
-            columns = df.columns.tolist()
-            mean_index = [i for i, s in enumerate(columns) if 'mean' in s][0]
-            if self.domain_type == 'mesh':
-                df = df[columns[:mean_index] + [('moment', '')] + columns[mean_index:-1]]
-            else:
-                df = df[columns[:mean_index] + ['moment'] + columns[mean_index:-1]]
+                # Place the moment column before the mean column
+                columns = df.columns.tolist()
+                mean_index \
+                    = [i for i, s in enumerate(columns) if 'mean' in s][0]
+                if self.domain_type == 'mesh':
+                    df = df[columns[:mean_index] + [('moment', '')] +
+                            columns[mean_index:-1]]
+                else:
+                    df = df[columns[:mean_index] + ['moment'] +
+                            columns[mean_index:-1]]
 
-        # Select rows corresponding to requested scattering moment
-        if moment != 'all':
-            cv.check_type('moment', moment, Integral)
-            cv.check_greater_than('moment', moment, 0, equality=True)
-            cv.check_less_than(
-                'moment', moment, self.legendre_order, equality=True)
-            df = df[df['moment'] == 'P{}'.format(moment)]
+            # Select rows corresponding to requested scattering moment
+            if moment != 'all':
+                cv.check_type('moment', moment, Integral)
+                cv.check_greater_than('moment', moment, 0, equality=True)
+                cv.check_less_than(
+                    'moment', moment, self.legendre_order, equality=True)
+                df = df[df['moment'] == 'P{}'.format(moment)]
+
+        elif self.scatter_format == 'histogram':
+            # Add a change-in-angle (mu) column to dataframe
+            ###TODO NOT SURE I NEED TO DO THIS
+            pass
 
         return df
 
@@ -3832,7 +3908,7 @@ class ScatterMatrixXS(MatrixMGXS):
 
         cv.check_value('xs_type', xs_type, ['macro', 'micro'])
 
-        if self.correction != 'P0':
+        if self.correction != 'P0' and self.scatter_format == 'legendre':
             rxn_type = '{0} (P{1})'.format(self.rxn_type, moment)
         else:
             rxn_type = self.rxn_type
diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py
index 35dff1223..5504b6e7c 100644
--- a/openmc/mgxs_library.py
+++ b/openmc/mgxs_library.py
@@ -1008,21 +1008,26 @@ class XSdata(object):
         check_type('temperature', temperature, Real)
         check_value('temperature', temperature, self.temperatures)
 
-        if self.scatter_format != 'legendre':
-            msg = 'Anisotropic scattering representations other than ' \
-                  'Legendre expansions have not yet been implemented in ' \
-                  'openmc.mgxs.'
-            raise ValueError(msg)
+        # Set the value of scatter_format based on the same value within
+        # scatter
+        self.scatter_format = scatter.scatter_format
 
         # If the user has not defined XSdata.order, then we will set
         # the order based on the data within scatter.
-        # Otherwise, we will check to see that XSdata.order to match
+        # Otherwise, we will check to see that XSdata.order matches
         # the order of scatter
-        if self.order is None:
-            self.order = scatter.legendre_order
-        else:
-            check_value('legendre_order', scatter.legendre_order,
-                        [self.order])
+        if self.scatter_format == 'legendre':
+            if self.order is None:
+                self.order = scatter.legendre_order
+            else:
+                check_value('legendre_order', scatter.legendre_order,
+                            [self.order])
+        elif self.scatter_format == 'histogram':
+            if self.order is None:
+                self.order = scatter.histogram_bins
+            else:
+                check_value('histogram_bins', scatter.histogram_bins,
+                            [self.order])
 
         i = np.where(self.temperatures == temperature)[0][0]
         if self.representation == 'isotropic':
@@ -1030,10 +1035,16 @@ class XSdata(object):
             self._scatter_matrix[i] = np.zeros((self.num_orders,
                                                 self.energy_groups.num_groups,
                                                 self.energy_groups.num_groups))
-            for moment in range(self.num_orders):
-                self._scatter_matrix[i][moment, :, :] = \
+            if self.scatter_format == 'legendre':
+                for moment in range(self.num_orders):
+                    self._scatter_matrix[i][moment, :, :] = \
+                        scatter.get_xs(nuclides=nuclide, xs_type=xs_type,
+                                       moment=moment, subdomains=subdomain)
+            else:
+                self._scatter_matrix[i][:, :, :] = \
                     scatter.get_xs(nuclides=nuclide, xs_type=xs_type,
-                                   moment=moment, subdomains=subdomain)
+                                   subdomains=subdomain)
+                import pdb; pdb.set_trace()
 
         elif self.representation == 'angle':
             msg = 'Angular-Dependent MGXS have not yet been implemented'
diff --git a/openmc/tallies.py b/openmc/tallies.py
index f9daa7a22..6922e77ac 100644
--- a/openmc/tallies.py
+++ b/openmc/tallies.py
@@ -1301,7 +1301,8 @@ class Tally(object):
 
                     # Create list of 2-tuples for energy boundary bins
                     elif isinstance(self_filter, (openmc.EnergyFilter,
-                        openmc.EnergyoutFilter)):
+                        openmc.EnergyoutFilter, openmc.MuFilter,
+                        openmc.PolarFilter, openmc.AzimuthalFilter)):
                         bins = []
                         for k in range(self_filter.num_bins):
                             bins.append((self_filter.bins[k], self_filter.bins[k+1]))

From f126d4387ada157e88126a3e013ec95e7dbde5d8 Mon Sep 17 00:00:00 2001
From: Adam Nelson <nelsonag@umich.edu>
Date: Tue, 18 Oct 2016 05:15:51 -0400
Subject: [PATCH 02/12] working on the writing of the MGXS lib

---
 openmc/filter.py       | 21 +++++++++++++++------
 openmc/mgxs/library.py | 36 ++++++++++++++----------------------
 openmc/mgxs/mgxs.py    | 19 +++++++++++++++----
 openmc/mgxs_library.py | 18 +++++++++++-------
 4 files changed, 55 insertions(+), 39 deletions(-)

diff --git a/openmc/filter.py b/openmc/filter.py
index 91ceb2931..fef7a25b0 100644
--- a/openmc/filter.py
+++ b/openmc/filter.py
@@ -868,15 +868,13 @@ class RealFilter(Filter):
             return np.allclose(self.bins, other.bins)
 
     def get_bin_index(self, filter_bin):
-        # Use lower energy bound to find index for RealFilters
-        deltas = np.abs(self.bins - filter_bin[1]) / filter_bin[1]
-        min_delta = np.min(deltas)
-        if min_delta < 1E-3:
-            return deltas.argmin() - 1
-        else:
+        i = np.where(self.bins == filter_bin[1])[0]
+        if len(i) == 0:
             msg = 'Unable to get the bin index for Filter since "{0}" ' \
                   'is not one of the bins'.format(filter_bin)
             raise ValueError(msg)
+        else:
+            return i[0] - 1
 
     def get_bin(self, bin_index):
         cv.check_type('bin_index', bin_index, Integral)
@@ -907,6 +905,17 @@ class EnergyFilter(RealFilter):
 
     """
 
+    def get_bin_index(self, filter_bin):
+        # Use lower energy bound to find index for RealFilters
+        deltas = np.abs(self.bins - filter_bin[1]) / filter_bin[1]
+        min_delta = np.min(deltas)
+        if min_delta < 1E-3:
+            return deltas.argmin() - 1
+        else:
+            msg = 'Unable to get the bin index for Filter since "{0}" ' \
+                  'is not one of the bins'.format(filter_bin)
+            raise ValueError(msg)
+
     def check_bins(self, bins):
         for edge in bins:
             if not isinstance(edge, Real):
diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py
index 96bfa6fab..46d404e4f 100644
--- a/openmc/mgxs/library.py
+++ b/openmc/mgxs/library.py
@@ -875,13 +875,13 @@ class Library(object):
         return pickle.load(open(full_filename, 'rb'))
 
     def get_xsdata(self, domain, xsdata_name, nuclide='total', xs_type='macro',
-                   order=None, subdomain=None):
+                   subdomain=None):
         """Generates an openmc.XSdata object describing a multi-group cross section
         dataset for writing to an openmc.MGXSLibrary object.
 
         Note that this method does not build an XSdata
         object with nested temperature tables.  The temperature of each
-        XSdata object will be left at the default value of 300K.
+        XSdata object will be left at the default value of 294K.
 
         Parameters
         ----------
@@ -896,10 +896,6 @@ class Library(object):
             Provide the macro or micro cross section in units of cm^-1 or
             barns. Defaults to 'macro'. If the Library object is not tallied by
             nuclide this will be set to 'macro' regardless.
-        order : int
-            Scattering order for this data entry.  Default is None,
-            which will set the XSdata object to use the order of the
-            Library.
         subdomain : iterable of int
             This parameter is not used unless using a mesh domain. In that
             case, the subdomain is an [i,j,k] index (1-based indexing) of the
@@ -929,10 +925,6 @@ class Library(object):
         cv.check_type('xsdata_name', xsdata_name, basestring)
         cv.check_type('nuclide', nuclide, basestring)
         cv.check_value('xs_type', xs_type, ['macro', 'micro'])
-        cv.check_type('order', order, (type(None), Integral))
-        if order is not None:
-            cv.check_greater_than('order', order, 0, equality=True)
-            cv.check_less_than('order', order, 10, equality=True)
         if subdomain is not None:
             cv.check_iterable_type('subdomain', subdomain, Integral,
                                    max_depth=3)
@@ -953,14 +945,6 @@ class Library(object):
             name += '_' + nuclide
         xsdata = openmc.XSdata(name, self.energy_groups)
 
-        if order is None:
-            # Set the order to the Library's order (the default behavior)
-            xsdata.order = self.legendre_order
-        else:
-            # Set the order of the xsdata object to the minimum of
-            # the provided order or the Library's order.
-            xsdata.order = min(order, self.legendre_order)
-
         # Right now only isotropic weighting is supported
         self.representation = 'isotropic'
 
@@ -1033,6 +1017,7 @@ class Library(object):
             using_multiplicity = False
 
         if using_multiplicity:
+            # import pdb; pdb.set_trace()
             nuscatt_mgxs = self.get_mgxs(domain, 'nu-scatter matrix')
             xsdata.set_scatter_matrix_mgxs(nuscatt_mgxs, xs_type=xs_type,
                                            nuclide=[nuclide],
@@ -1049,10 +1034,17 @@ class Library(object):
                 # accounted for approximately by using an adjusted
                 # absorption cross section.
                 if 'total' in self.mgxs_types or 'transport' in self.mgxs_types:
-                    for i in range(len(xsdata.temperatures)):
-                        xsdata._absorption[i] = \
-                            np.subtract(xsdata._total[i], np.sum(
-                                xsdata._scatter_matrix[i][0, :, :], axis=1))
+                    if xsdata.scatter_format == 'legendre':
+                        for i in range(len(xsdata.temperatures)):
+                            xsdata._absorption[i] = \
+                                np.subtract(xsdata._total[i], np.sum(
+                                    xsdata._scatter_matrix[i][0, :, :], axis=1))
+                    elif xsdata.scatter_format == 'histogram':
+                        for i in range(len(xsdata.temperatures)):
+                            xsdata._absorption[i] = \
+                                np.subtract(xsdata._total[i], np.sum(np.sum(
+                                    xsdata._scatter_matrix[i][:, :, :], axis=0),
+                                    axis=1))
 
         return xsdata
 
diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py
index cafda1ecf..c30cc87ac 100644
--- a/openmc/mgxs/mgxs.py
+++ b/openmc/mgxs/mgxs.py
@@ -3410,7 +3410,8 @@ class ScatterMatrixXS(MatrixMGXS):
             else:
                 filters = [[energy], [energy, energyout]]
         elif self.scatter_format == 'histogram':
-            bins = np.linspace(-1., 1., num=self.histogram_bins, endpoint=True)
+            bins = np.linspace(-1., 1., num=self.histogram_bins + 1,
+                               endpoint=True)
             filters = [[energy], [energy, energyout, openmc.MuFilter(bins)]]
 
         return filters
@@ -3758,9 +3759,19 @@ class ScatterMatrixXS(MatrixMGXS):
         else:
             num_out_groups = len(out_groups)
 
+        if self.scatter_format == 'histogram':
+            num_mu_bins = self.histogram_bins
+        else:
+            num_mu_bins = 1
+
         # Reshape tally data array with separate axes for domain and energy
-        num_subdomains = int(xs.shape[0] / (num_in_groups * num_out_groups))
-        new_shape = (num_subdomains, num_in_groups, num_out_groups)
+        num_subdomains = int(xs.shape[0] /
+                             (num_mu_bins * num_in_groups * num_out_groups))
+        if self.scatter_format == 'histogram':
+            new_shape = (num_subdomains, num_in_groups, num_out_groups,
+                         num_mu_bins)
+        else:
+            new_shape = (num_subdomains, num_in_groups, num_out_groups)
         new_shape += xs.shape[1:]
         xs = np.reshape(xs, new_shape)
 
@@ -3771,7 +3782,7 @@ class ScatterMatrixXS(MatrixMGXS):
         # Reverse data if user requested increasing energy groups since
         # tally data is stored in order of increasing energies
         if order_groups == 'increasing':
-            xs = xs[:, ::-1, ::-1, :]
+            xs = xs[:, ::-1, ::-1, ...]
 
         if squeeze:
             xs = np.squeeze(xs)
diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py
index 5504b6e7c..c78b9ff0c 100644
--- a/openmc/mgxs_library.py
+++ b/openmc/mgxs_library.py
@@ -1031,21 +1031,21 @@ class XSdata(object):
 
         i = np.where(self.temperatures == temperature)[0][0]
         if self.representation == 'isotropic':
-            # Get the scattering orders in the outermost dimension
-            self._scatter_matrix[i] = np.zeros((self.num_orders,
-                                                self.energy_groups.num_groups,
-                                                self.energy_groups.num_groups))
             if self.scatter_format == 'legendre':
+                # Get the scattering orders in the outermost dimension
+                self._scatter_matrix[i] = np.zeros((self.num_orders,
+                                                    self.energy_groups.num_groups,
+                                                    self.energy_groups.num_groups))
                 for moment in range(self.num_orders):
                     self._scatter_matrix[i][moment, :, :] = \
                         scatter.get_xs(nuclides=nuclide, xs_type=xs_type,
                                        moment=moment, subdomains=subdomain)
             else:
-                self._scatter_matrix[i][:, :, :] = \
+                self._scatter_matrix[i] = \
                     scatter.get_xs(nuclides=nuclide, xs_type=xs_type,
                                    subdomains=subdomain)
-                import pdb; pdb.set_trace()
-
+                # self._scatter_matrix[i] = \
+                #     np.swapaxes(self._scatter_matrix[i], 0, 2)
         elif self.representation == 'angle':
             msg = 'Angular-Dependent MGXS have not yet been implemented'
             raise ValueError(msg)
@@ -1124,6 +1124,10 @@ class XSdata(object):
                 scatt = scatter.get_xs(nuclides=nuclide,
                                        xs_type=xs_type, moment=0,
                                        subdomains=subdomain)
+                if scatter.scatter_format == 'histogram':
+                    scatt = np.sum(scatt, axis=4)
+                if nuscatter.scatter_format == 'histogram':
+                    nuscatt = np.sum(nuscatt, axis=4)
                 self._multiplicity_matrix[i] = np.divide(nuscatt, scatt)
         elif self.representation == 'angle':
             msg = 'Angular-Dependent MGXS have not yet been implemented'

From 45fb4201f643cb49045e3af36c8a771953d5e995 Mon Sep 17 00:00:00 2001
From: Adam Nelson <nelsonag@umich.edu>
Date: Sun, 23 Oct 2016 13:14:46 -0400
Subject: [PATCH 03/12] Fixed angular indexing

---
 openmc/mgxs/mgxs.py    | 5 +++++
 openmc/mgxs_library.py | 6 ++----
 2 files changed, 7 insertions(+), 4 deletions(-)

diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py
index fb81e4c31..470fd01d0 100644
--- a/openmc/mgxs/mgxs.py
+++ b/openmc/mgxs/mgxs.py
@@ -3776,6 +3776,11 @@ class ScatterMatrixXS(MatrixMGXS):
         if order_groups == 'increasing':
             xs = xs[:, ::-1, ::-1, ...]
 
+        # Place the histogram bins before the outgoing group index
+        if self.scatter_format == 'histogram':
+            xs = np.swapaxes(xs, 3, 1)
+            xs = np.swapaxes(xs, 3, 2)
+
         if squeeze:
             xs = np.squeeze(xs)
             xs = np.atleast_2d(xs)
diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py
index 8a7d771ac..6256da352 100644
--- a/openmc/mgxs_library.py
+++ b/openmc/mgxs_library.py
@@ -1043,8 +1043,6 @@ class XSdata(object):
                 self._scatter_matrix[i] = \
                     scatter.get_xs(nuclides=nuclide, xs_type=xs_type,
                                    subdomains=subdomain)
-                # self._scatter_matrix[i] = \
-                #     np.swapaxes(self._scatter_matrix[i], 0, 2)
         elif self.representation == 'angle':
             msg = 'Angular-Dependent MGXS have not yet been implemented'
             raise ValueError(msg)
@@ -1124,9 +1122,9 @@ class XSdata(object):
                                        xs_type=xs_type, moment=0,
                                        subdomains=subdomain)
                 if scatter.scatter_format == 'histogram':
-                    scatt = np.sum(scatt, axis=4)
+                    scatt = np.sum(scatt, axis=0)
                 if nuscatter.scatter_format == 'histogram':
-                    nuscatt = np.sum(nuscatt, axis=4)
+                    nuscatt = np.sum(nuscatt, axis=0)
                 self._multiplicity_matrix[i] = np.divide(nuscatt, scatt)
         elif self.representation == 'angle':
             msg = 'Angular-Dependent MGXS have not yet been implemented'

From 50a40795976ea5db3a7e34110b9a105115800b5c Mon Sep 17 00:00:00 2001
From: Adam Nelson <nelsonag@umich.edu>
Date: Sat, 12 Nov 2016 09:40:07 -0500
Subject: [PATCH 04/12] minor changes

---
 openmc/mgxs/mgxs.py    |  4 ++--
 openmc/mgxs_library.py | 15 +++++++++++----
 2 files changed, 13 insertions(+), 6 deletions(-)

diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py
index 2fb7459b0..324cbe56c 100644
--- a/openmc/mgxs/mgxs.py
+++ b/openmc/mgxs/mgxs.py
@@ -3776,10 +3776,10 @@ class ScatterMatrixXS(MatrixMGXS):
         if order_groups == 'increasing':
             xs = xs[:, ::-1, ::-1, ...]
 
-        # Place the histogram bins before the outgoing group index
+        # Place the histogram bins before the incoming groups
         if self.scatter_format == 'histogram':
-            xs = np.swapaxes(xs, 3, 1)
             xs = np.swapaxes(xs, 3, 2)
+            xs = np.swapaxes(xs, 2, 1)
 
         if squeeze:
             xs = np.squeeze(xs)
diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py
index 290b25424..e3af59632 100644
--- a/openmc/mgxs_library.py
+++ b/openmc/mgxs_library.py
@@ -1738,11 +1738,18 @@ class XSdata(object):
             if self.representation == 'isotropic':
 
                 g_out_bounds = np.zeros((G, 2), dtype=np.int)
-
                 for g_in in range(G):
-                    nz = np.nonzero(self._scatter_matrix[i][0, g_in, :])
-                    g_out_bounds[g_in, 0] = nz[0][0]
-                    g_out_bounds[g_in, 1] = nz[0][-1]
+                    if self.scatter_format == 'legendre':
+                        nz = np.nonzero(self._scatter_matrix[i][0, g_in, :])
+                    elif self.scatter_format == 'histogram':
+                        nz = np.nonzero(np.sum(
+                            self._scatter_matrix[i][:, g_in, :], axis=0))
+                    if len(nz[0]) == 0:
+                        g_out_bounds[g_in, 0] = 0
+                        g_out_bounds[g_in, 1] = 0
+                    else:
+                        g_out_bounds[g_in, 0] = nz[0][0]
+                        g_out_bounds[g_in, 1] = nz[0][-1]
 
                 # Now create the flattened scatter matrix array
                 matrix = self._scatter_matrix[i]

From 2455388315919dca1080fd060d8e2765eea5fd92 Mon Sep 17 00:00:00 2001
From: Adam Nelson <nelsonag@umich.edu>
Date: Sun, 13 Nov 2016 13:29:45 -0500
Subject: [PATCH 05/12] Fixed array ordering of scattering matrix

---
 openmc/mgxs/library.py     | 25 +++++++++++++++++-----
 openmc/mgxs/mgxs.py        | 17 ++++++++-------
 openmc/mgxs_library.py     | 43 ++++++++++++++++++++------------------
 scripts/openmc-update-mgxs | 40 +++++++++++++++++++++++------------
 src/mgxs_header.F90        |  6 +++---
 src/scattdata_header.F90   |  4 ++--
 6 files changed, 84 insertions(+), 51 deletions(-)

diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py
index 05cabf496..97bb6d96b 100644
--- a/openmc/mgxs/library.py
+++ b/openmc/mgxs/library.py
@@ -287,7 +287,7 @@ class Library(object):
     def by_nuclide(self, by_nuclide):
         cv.check_type('by_nuclide', by_nuclide, bool)
 
-        if by_nuclide == True and self.domain_type == 'mesh':
+        if by_nuclide and self.domain_type == 'mesh':
             raise ValueError('Unable to create MGXS library by nuclide with '
                              'mesh domain')
 
@@ -297,7 +297,7 @@ class Library(object):
     def domain_type(self, domain_type):
         cv.check_value('domain type', domain_type, openmc.mgxs.DOMAIN_TYPES)
 
-        if self.by_nuclide == True and domain_type == 'mesh':
+        if self.by_nuclide and domain_type == 'mesh':
             raise ValueError('Unable to create MGXS library by nuclide with '
                              'mesh domain')
 
@@ -372,7 +372,15 @@ class Library(object):
 
     @scatter_format.setter
     def scatter_format(self, scatter_format):
-        cv.check_value('scatter_format', scatter_format, openmc.mgxs.MU_TREATMENTS)
+        cv.check_value('scatter_format', scatter_format,
+                       openmc.mgxs.MU_TREATMENTS)
+
+        if scatter_format == 'histogram' and self.correction == 'P0':
+            msg = 'The P0 correction will be ignored since the ' \
+                  'scatter format is set to histogram'
+            warn(msg)
+            self.correction = None
+
         self._scatter_format = scatter_format
 
     @legendre_order.setter
@@ -405,6 +413,13 @@ class Library(object):
             msg = 'The histogram bins will be ignored since the ' \
                   'scatter format is set to legendre'
             warn(msg)
+        elif self.scatter_format == 'histogram':
+            if self.correction == 'P0':
+                msg = 'The P0 correction will be ignored since ' \
+                      'a histogram representation of the scattering '\
+                      'kernel is requested'
+                warn(msg, RuntimeWarning)
+                self.correction = None
 
         self._histogram_bins = histogram_bins
 
@@ -471,7 +486,7 @@ class Library(object):
                         mgxs.delayed_groups = None
                     else:
                         delayed_groups \
-                            = list(range(1,self.num_delayed_groups+1))
+                            = list(range(1, self.num_delayed_groups + 1))
                         mgxs.delayed_groups = delayed_groups
 
                 # If a tally trigger was specified, add it to the MGXS
@@ -1351,7 +1366,7 @@ class Library(object):
                          'scattering matrix is not provided.')
         # Total or transport can be present, but if using
         # self.correction=="P0", then we should use transport.
-        if (((self.correction is "P0") and
+        if (((self.correction == "P0") and
              ('nu-transport' not in self.mgxs_types))):
             error_flag = True
             warn('A "nu-transport" MGXS type is required since a "P0" '
diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py
index 324cbe56c..036feac74 100644
--- a/openmc/mgxs/mgxs.py
+++ b/openmc/mgxs/mgxs.py
@@ -3775,15 +3775,16 @@ class ScatterMatrixXS(MatrixMGXS):
         # tally data is stored in order of increasing energies
         if order_groups == 'increasing':
             xs = xs[:, ::-1, ::-1, ...]
-
-        # Place the histogram bins before the incoming groups
-        if self.scatter_format == 'histogram':
-            xs = np.swapaxes(xs, 3, 2)
-            xs = np.swapaxes(xs, 2, 1)
-
+        # import pdb; pdb.set_trace()
         if squeeze:
-            xs = np.squeeze(xs)
-            xs = np.atleast_2d(xs)
+            # We want to squeeze out everything but the in_groups, out_groups,
+            # and, if needed, num_mu_bins dimension. These must not be squeezed
+            # so 1-group problems have the correct shape.
+            if self.scatter_format == 'histogram':
+                axes = (5, 4, 0)
+            else:
+                axes = (4, 3, 0)
+            xs = np.squeeze(xs, axis=axes)
 
         return xs
 
diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py
index e3af59632..5256dd102 100644
--- a/openmc/mgxs_library.py
+++ b/openmc/mgxs_library.py
@@ -15,7 +15,7 @@ from openmc.checkvalue import check_type, check_value, check_greater_than, \
 # Supported incoming particle MGXS angular treatment representations
 _REPRESENTATIONS = ['isotropic', 'angle']
 _SCATTER_TYPES = ['tabular', 'legendre', 'histogram']
-_XS_SHAPES = ["[Order][G][G']", "[G]", "[G']", "[G][G']", "[DG]", "[G][DG]",
+_XS_SHAPES = ["[G][G'][Order]", "[G]", "[G']", "[G][G']", "[DG]", "[G][DG]",
               "[G'][DG]"]
 
 
@@ -134,7 +134,7 @@ class XSdata(object):
     Note that some cross sections can be input in more than one shape so they
     are listed multiple times:
 
-    [Order][G][G']: scatter_matrix
+    [G][G'][Order]: scatter_matrix
 
     [G]: total, absorption, fission, kappa_fission, nu_fission,
          prompt_nu_fission, inverse_velocity
@@ -298,9 +298,9 @@ class XSdata(object):
                                           self.num_delayed_groups)
             self._xs_shapes["[G'][DG]"] = (self.energy_groups.num_groups,
                                            self.num_delayed_groups)
-            self._xs_shapes["[Order][G][G']"] \
-                = (self.num_orders, self.energy_groups.num_groups,
-                   self.energy_groups.num_groups)
+            self._xs_shapes["[G][G'][Order]"] \
+                = (self.energy_groups.num_groups,
+                   self.energy_groups.num_groups, self.num_orders)
 
             # If representation is by angle prepend num polar and num azim
             if self.representation == 'angle':
@@ -718,7 +718,7 @@ class XSdata(object):
         """
 
         # Get the accepted shapes for this xs
-        shapes = [self.xs_shapes["[Order][G][G']"]]
+        shapes = [self.xs_shapes["[G][G'][Order]"]]
 
         # Convert to a numpy array so we can easily get the shape for checking
         scatter = np.asarray(scatter)
@@ -1534,11 +1534,10 @@ class XSdata(object):
         if self.representation == 'isotropic':
             if self.scatter_format == 'legendre':
                 # Get the scattering orders in the outermost dimension
-                self._scatter_matrix[i] = np.zeros((self.num_orders,
-                                                    self.energy_groups.num_groups,
-                                                    self.energy_groups.num_groups))
+                self._scatter_matrix[i] = \
+                    np.zeros(self.xs_shapes["[G][G'][Order]"])
                 for moment in range(self.num_orders):
-                    self._scatter_matrix[i][moment, :, :] = \
+                    self._scatter_matrix[i][:, :, moment] = \
                         scatter.get_xs(nuclides=nuclide, xs_type=xs_type,
                                        moment=moment, subdomains=subdomain)
             else:
@@ -1657,7 +1656,7 @@ class XSdata(object):
                 if self.num_polar is not None:
                     grp.attrs['num_polar'] = self.num_polar
 
-        grp.attrs['scatter_shape'] = np.string_("[Order][G][G']")
+        grp.attrs['scatter_shape'] = np.string_("[G][G'][Order]")
         if self.scatter_format is not None:
             grp.attrs['scatter_format'] = np.string_(self.scatter_format)
         if self.order is not None:
@@ -1736,14 +1735,13 @@ class XSdata(object):
             # Get the sparse scattering data to print to the library
             G = self.energy_groups.num_groups
             if self.representation == 'isotropic':
-
                 g_out_bounds = np.zeros((G, 2), dtype=np.int)
                 for g_in in range(G):
                     if self.scatter_format == 'legendre':
-                        nz = np.nonzero(self._scatter_matrix[i][0, g_in, :])
+                        nz = np.nonzero(self._scatter_matrix[i][g_in, :, 0])
                     elif self.scatter_format == 'histogram':
                         nz = np.nonzero(np.sum(
-                            self._scatter_matrix[i][:, g_in, :], axis=0))
+                            self._scatter_matrix[i][g_in, :, :], axis=1))
                     if len(nz[0]) == 0:
                         g_out_bounds[g_in, 0] = 0
                         g_out_bounds[g_in, 1] = 0
@@ -1757,8 +1755,8 @@ class XSdata(object):
                 for g_in in range(G):
                     for g_out in range(g_out_bounds[g_in, 0],
                                        g_out_bounds[g_in, 1] + 1):
-                        for l in range(len(matrix[:, g_in, g_out])):
-                            flat_scatt.append(matrix[l, g_in, g_out])
+                        for l in range(len(matrix[g_in, g_out, :])):
+                            flat_scatt.append(matrix[g_in, g_out, l])
 
                 # And write it.
                 scatt_grp = xs_grp.create_group('scatter_data')
@@ -1767,7 +1765,6 @@ class XSdata(object):
 
                 # Repeat for multiplicity
                 if self._multiplicity_matrix[i] is not None:
-
                     # Now create the flattened scatter matrix array
                     matrix = self._multiplicity_matrix[i][:, :]
                     flat_mult = []
@@ -1793,7 +1790,13 @@ class XSdata(object):
                 for p in range(Np):
                     for a in range(Na):
                         for g_in in range(G):
-                            matrix = self._scatter_matrix[i][p, a, 0, g_in, :]
+                            if self.scatter_format == 'legendre':
+                                matrix = \
+                                    self._scatter_matrix[i][p, a, g_in, :, 0]
+                            elif self.scatter_format == 'histogram':
+                                matrix = \
+                                    np.sum(self._scatter_matrix[i][p, a, g_in, :, :],
+                                           axis=1)
                             nz = np.nonzero(matrix)
                             g_out_bounds[p, a, g_in, 0] = nz[0][0]
                             g_out_bounds[p, a, g_in, 1] = nz[0][-1]
@@ -1806,8 +1809,8 @@ class XSdata(object):
                         for g_in in range(G):
                             for g_out in range(g_out_bounds[p, a, g_in, 0],
                                                g_out_bounds[p, a, g_in, 1] + 1):
-                                for l in range(len(matrix[:, g_in, g_out])):
-                                    flat_scatt.append(matrix[l, g_in, g_out])
+                                for l in range(len(matrix[g_in, g_out, :])):
+                                    flat_scatt.append(matrix[g_in, g_out, l])
 
                 # And write it.
                 scatt_grp = xs_grp.create_group('scatter_data')
diff --git a/scripts/openmc-update-mgxs b/scripts/openmc-update-mgxs
index cad48d0e3..54a553f77 100755
--- a/scripts/openmc-update-mgxs
+++ b/scripts/openmc-update-mgxs
@@ -23,12 +23,10 @@ optional arguments:
 
 from __future__ import print_function
 import os
-from shutil import move
 import warnings
 import xml.etree.ElementTree as ET
 
 import argparse
-import h5py
 import numpy as np
 
 import openmc.mgxs_library
@@ -51,7 +49,7 @@ def parse_args():
 
     if args['output'] == '':
         filename = args['input'].name
-        extension = filenameos.path.splitext()
+        extension = os.path.splitext(filename)
         if extension == '.xml':
             filename = filename[:filename.rfind('.')] + '.h5'
         args['output'] = filename
@@ -84,6 +82,8 @@ if __name__ == '__main__':
     temp = tree.find('group_structure').text.strip()
     temp = np.array(temp.split())
     group_structure = temp.astype(np.float)
+    # Convert from MeV to eV
+    group_structure *= 1.e6
     energy_groups = openmc.mgxs.EnergyGroups(group_structure)
     temp = tree.find('inverse-velocity')
     if temp is not None:
@@ -103,7 +103,7 @@ if __name__ == '__main__':
         temperature = get_data(xsdata_elem, 'kT')
         if temperature is not None:
             temperature = \
-                float(temperature) / openmc.data.K_BOLTZMANN
+                float(temperature) / openmc.data.K_BOLTZMANN * 1.E6
         else:
             temperature = 294.
         temperatures = [temperature]
@@ -163,7 +163,7 @@ if __name__ == '__main__':
         if temp is not None:
             temp = np.array(temp.split(), dtype=float)
             total = temp.astype(np.float)
-            total.shape = xsd[i].vector_shape
+            total.shape = xsd[i].xs_shapes['[G]']
             xsd[i].set_total(total, temperature)
 
         if inverse_velocity is not None:
@@ -172,41 +172,55 @@ if __name__ == '__main__':
         temp = get_data(xsdata_elem, 'absorption')
         temp = np.array(temp.split())
         absorption = temp.astype(np.float)
-        absorption.shape = xsd[i].vector_shape
+        absorption.shape = xsd[i].xs_shapes['[G]']
         xsd[i].set_absorption(absorption, temperature)
 
         temp = get_data(xsdata_elem, 'scatter')
         temp = np.array(temp.split())
         scatter = temp.astype(np.float)
-        scatter.shape = xsd[i].pn_matrix_shape
+        # This is now a flattened-array of something that started with a
+        # shape of [Order][G][G']; we need to unflatten and then switch the
+        # ordering
+        in_shape = (order_dim, energy_groups.num_groups,
+                    energy_groups.num_groups)
+        if representation == 'angle':
+            in_shape = (n_pol, n_azi) + in_shape
+            scatter.shape = in_shape
+            scatter = np.swapaxes(scatter, 2, 3)
+            scatter = np.swapaxes(scatter, 3, 4)
+        else:
+            scatter.shape = in_shape
+            scatter = np.swapaxes(scatter, 0, 1)
+            scatter = np.swapaxes(scatter, 1, 2)
+
         xsd[i].set_scatter_matrix(scatter, temperature)
 
         temp = get_data(xsdata_elem, 'multiplicity')
         if temp is not None:
             temp = np.array(temp.split())
             multiplicity = temp.astype(np.float)
-            multiplicity.shape = xsd[i].matrix_shape
+            multiplicity.shape = xsd[i].xs_shapes["[G][G']"]
             xsd[i].set_multiplicity_matrix(multiplicity, temperature)
 
         temp = get_data(xsdata_elem, 'fission')
         if temp is not None:
             temp = np.array(temp.split())
             fission = temp.astype(np.float)
-            fission.shape = xsd[i].vector_shape
+            fission.shape = xsd[i].xs_shapes['[G]']
             xsd[i].set_fission(fission, temperature)
 
         temp = get_data(xsdata_elem, 'kappa_fission')
         if temp is not None:
             temp = np.array(temp.split())
             kappa_fission = temp.astype(np.float)
-            kappa_fission.shape = xsd[i].vector_shape
+            kappa_fission.shape = xsd[i].xs_shapes['[G]']
             xsd[i].set_kappa_fission(kappa_fission, temperature)
 
         temp = get_data(xsdata_elem, 'chi')
         if temp is not None:
             temp = np.array(temp.split())
             chi = temp.astype(np.float)
-            chi.shape = xsd[i].vector_shape
+            chi.shape = xsd[i].xs_shapes['[G]']
             xsd[i].set_chi(chi, temperature)
         else:
             chi = None
@@ -216,9 +230,9 @@ if __name__ == '__main__':
             temp = np.array(temp.split())
             nu_fission = temp.astype(np.float)
             if chi is not None:
-                nu_fission.shape = xsd[i].vector_shape
+                nu_fission.shape = xsd[i].xs_shapes['[G]']
             else:
-                nu_fission.shape = xsd[i].matrix_shape
+                nu_fission.shape = xsd[i].xs_shapes["[G][G']"]
             xsd[i].set_nu_fission(nu_fission, temperature)
 
     # Build library as we go, but first we have enough to initialize it
diff --git a/src/mgxs_header.F90 b/src/mgxs_header.F90
index 17104d508..4fb1f7272 100644
--- a/src/mgxs_header.F90
+++ b/src/mgxs_header.F90
@@ -355,7 +355,7 @@ module mgxs_header
       if (attribute_exists(xs_id, "scatter_shape")) then
         call read_attribute(temp_str, xs_id, "scatter_shape")
         temp_str = trim(temp_str)
-        if (to_lower(temp_str) /= "[order][g][g']") then
+        if (to_lower(temp_str) /= "[g][g'][order]") then
           call fatal_error("Invalid scatter_shape option!")
         end if
       end if
@@ -1794,8 +1794,8 @@ module mgxs_header
             order_dim = order + 1
           end if
 
-          ! Convert temp_1d to a jagged array ((gin) % data(l, gout)) for passing
-          ! to ScattData
+          ! Convert temp_1d to a jagged array ((gin) % data(l, gout)) for
+          ! passing to ScattData
           allocate(input_scatt(energy_groups, this % n_azi, this % n_pol))
 
           index = 1
diff --git a/src/scattdata_header.F90 b/src/scattdata_header.F90
index 6004cbd19..684be88b8 100644
--- a/src/scattdata_header.F90
+++ b/src/scattdata_header.F90
@@ -274,7 +274,7 @@ contains
       ! Get this by summing the un-normalized P0 coefficient in matrix
       ! over all outgoing groups
       do gin = 1, groups
-        this % scattxs(gin) = sum(matrix(gin) % data(1, :), dim=1)
+        this % scattxs(gin) = sum(matrix(gin) % data(:, :))
       end do
 
       allocate(energy(groups))
@@ -282,7 +282,7 @@ contains
       ! while also normalizing matrix itself (making CDF of f(mu=1)=1)
       do gin = 1, groups
         allocate(energy(gin) % data(gmin(gin):gmax(gin)))
-        do gout = 1, groups
+        do gout = gmin(gin), gmax(gin)
           norm = sum(matrix(gin) % data(:, gout))
           energy(gin) % data(gout) = norm
           if (norm /= ZERO) then

From 7c4ed737afe8b9ece0b35523bd7c7030962548a6 Mon Sep 17 00:00:00 2001
From: Adam Nelson <nelsonag@umich.edu>
Date: Sun, 13 Nov 2016 14:38:57 -0500
Subject: [PATCH 06/12] Fixing failing tests

---
 openmc/mgxs/mgxs.py                 |   9 ++++-
 tests/1d_mgxs.h5                    | Bin 128136 -> 134280 bytes
 tests/test_tallies/inputs_true.dat  |   2 +-
 tests/test_tallies/results_true.dat |   2 +-
 tests/test_tallies/test_tallies.py  |  51 ++++++++++------------------
 5 files changed, 28 insertions(+), 36 deletions(-)

diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py
index 036feac74..1337b1fa7 100644
--- a/openmc/mgxs/mgxs.py
+++ b/openmc/mgxs/mgxs.py
@@ -3784,7 +3784,14 @@ class ScatterMatrixXS(MatrixMGXS):
                 axes = (5, 4, 0)
             else:
                 axes = (4, 3, 0)
-            xs = np.squeeze(xs, axis=axes)
+            # Squeeze will return a ValueError if the axis has a size greater
+            # than 1, so try each axis in axes one at a time, catching the
+            # ValueError as needed.
+            for axis in axes:
+                try:
+                    xs = np.squeeze(xs, axis=axis)
+                except ValueError:
+                    pass
 
         return xs
 
diff --git a/tests/1d_mgxs.h5 b/tests/1d_mgxs.h5
index e3b85c24163fe4020435db5c5d674575f2b65be0..8f28e02cbcfb2d51e94d81bafdd97ba1220313d3 100644
GIT binary patch
delta 1872
zcmeBp$=)%OV}d3V57S2N^&FEQ>6%Udp}}a#>mD2Jt{xlhUzCzs6ua4vGeiZ4-1K|-
zj82;)^o%5MDA@eZZ3$2=9-G0F{Wz^Q2P9mQ!66SZ{6_67aU2T31{6%Yh12qV>!0Z0
za2V9(cRwlOP_TL5e;Z947H`kfWc=fbOA;6+XDo0@%B^I~@($)<@L*tJU|?Wm;A0SB
z;9w|DPAn-&Es8JBNGwPNic2ywG6(=^5M~5o2A~=cfUp^Wx)>Nb{y=yTQefgiQK}lt
z1h$as#xm2!SWsXBgNxdMz)Blq*=S=dJ8g`GBt5E!EhlYU%S9VwxoKl855(Aw2cI);
zyx_-J@9ymH1J2N(oDT%Z86L!sMk$&Z7$CV6QXpWHhh#V0@{oK3m#@dL0FndX3NYj$
zNg201Bu(O$ha@hTd_ATQA;|=$08<_k-?-%=5r<nI5(7|qSu8;e3RVy(;Df}{4PZ=R
H!}*K=*xXlq

delta 457
zcmeB}$<gtWeS#*_7sid+>p3P1aq9E?7p0^Y#YVdWk$Nlx!{i0J7L%_jOHI7Mwz+`Q
zNQHo!$qHHun<Mm$BnYV6tl-`xG5MCRJ}-)Wljp@-L2a8{z$vjgAmNe>epMj*Zq&XK
z$FB-xR>8De1YPxD{UIF!E{FK=?k7b8>NY?4FQtjsRof?MGTw5<t7KaN<0=chDz-bU
IWc2U`08!MZ3jhEB

diff --git a/tests/test_tallies/inputs_true.dat b/tests/test_tallies/inputs_true.dat
index 69472c8d6..a8b40554b 100644
--- a/tests/test_tallies/inputs_true.dat
+++ b/tests/test_tallies/inputs_true.dat
@@ -1 +1 @@
-b1f616bc0e342e3886b4b27eca22c4d1e55e445979254d5ba5f9ea01048f53d0e9a01c575745ae130eb34108080d2256ab9e4cb926ea0050891a44ae37ecdd88
\ No newline at end of file
+ac118c5796593df0f53491920875bbe77921882253e7f24942e365e5812d726fb3d3bab86da02b3aad93f6f8f20a090125ffa51f9fdf9aef84009702f9cc363d
\ No newline at end of file
diff --git a/tests/test_tallies/results_true.dat b/tests/test_tallies/results_true.dat
index 639e8baf3..1d6f6bfaf 100644
--- a/tests/test_tallies/results_true.dat
+++ b/tests/test_tallies/results_true.dat
@@ -1 +1 @@
-ceb1420d48e097185dd0798a495676c8e968c37c6e75fd6232137a4df8f6fae33d8232b95c68f6479815a014d20ccfc7611f96b22fcb7da159b42346e00cdf4e
\ No newline at end of file
+a8172f7492fcc69d21f53c20523ffd0ad2882a753a3e1d263d135153cc6b3ee8d8e38f8823a4a2eb3288d730579892875d2d19f85f083dd659f91414988115ea
\ No newline at end of file
diff --git a/tests/test_tallies/test_tallies.py b/tests/test_tallies/test_tallies.py
index 9b48729ef..9e525d2cf 100644
--- a/tests/test_tallies/test_tallies.py
+++ b/tests/test_tallies/test_tallies.py
@@ -3,6 +3,7 @@
 import os
 import sys
 sys.path.insert(0, os.pardir)
+import numpy as np
 from testing_harness import PyAPITestHarness
 from openmc.filter import *
 from openmc import Mesh, Tally, Tallies
@@ -21,33 +22,27 @@ class TalliesTestHarness(PyAPITestHarness):
         self._input_set.settings.source = Source(space=Box(
             [-160, -160, -183], [160, 160, 183]))
 
-        azimuthal_bins = (-3.1416, -1.8850, -0.6283, 0.6283, 1.8850, 3.1416)
-        azimuthal_filter1 = AzimuthalFilter(azimuthal_bins)
+        azimuthal_bins = (-3.14159, -1.8850, -0.6283, 0.6283, 1.8850, 3.14159)
+        azimuthal_filter = AzimuthalFilter(azimuthal_bins)
         azimuthal_tally1 = Tally()
-        azimuthal_tally1.filters = [azimuthal_filter1]
+        azimuthal_tally1.filters = [azimuthal_filter]
         azimuthal_tally1.scores = ['flux']
         azimuthal_tally1.estimator = 'tracklength'
 
         azimuthal_tally2 = Tally()
-        azimuthal_tally2.filters = [azimuthal_filter1]
+        azimuthal_tally2.filters = [azimuthal_filter]
         azimuthal_tally2.scores = ['flux']
         azimuthal_tally2.estimator = 'analog'
 
-        azimuthal_filter2 = AzimuthalFilter(5)
-        azimuthal_tally3 = Tally()
-        azimuthal_tally3.filters = [azimuthal_filter2]
-        azimuthal_tally3.scores = ['flux']
-        azimuthal_tally3.estimator = 'tracklength'
-
         mesh_2x2 = Mesh(mesh_id=1)
         mesh_2x2.lower_left  = [-182.07, -182.07]
         mesh_2x2.upper_right = [182.07,  182.07]
         mesh_2x2.dimension = [2, 2]
         mesh_filter = MeshFilter(mesh_2x2)
-        azimuthal_tally4 = Tally()
-        azimuthal_tally4.filters = [azimuthal_filter2, mesh_filter]
-        azimuthal_tally4.scores = ['flux']
-        azimuthal_tally4.estimator = 'tracklength'
+        azimuthal_tally3 = Tally()
+        azimuthal_tally3.filters = [azimuthal_filter, mesh_filter]
+        azimuthal_tally3.scores = ['flux']
+        azimuthal_tally3.estimator = 'tracklength'
 
         cellborn_tally = Tally()
         cellborn_tally.filters = [CellbornFilter((10, 21, 22, 23))]
@@ -76,20 +71,17 @@ class TalliesTestHarness(PyAPITestHarness):
         material_tally.filters = [MaterialFilter((1, 2, 3, 4))]
         material_tally.scores = ['total']
 
+        mu_bins = (-1.0, -0.5, 0.0, 0.5, 1.0)
+        mu_filter = MuFilter(mu_bins)
         mu_tally1 = Tally()
-        mu_tally1.filters = [MuFilter((-1.0, -0.5, 0.0, 0.5, 1.0))]
+        mu_tally1.filters = [mu_filter]
         mu_tally1.scores = ['scatter', 'nu-scatter']
 
-        mu_filter = MuFilter(5)
         mu_tally2 = Tally()
-        mu_tally2.filters = [mu_filter]
+        mu_tally2.filters = [mu_filter, mesh_filter]
         mu_tally2.scores = ['scatter', 'nu-scatter']
 
-        mu_tally3 = Tally()
-        mu_tally3.filters = [mu_filter, mesh_filter]
-        mu_tally3.scores = ['scatter', 'nu-scatter']
-
-        polar_bins = (0.0, 0.6283, 1.2566, 1.8850, 2.5132, 3.1416)
+        polar_bins = (0.0, 0.6283, 1.2566, 1.8850, 2.5132, 3.14159)
         polar_filter = PolarFilter(polar_bins)
         polar_tally1 = Tally()
         polar_tally1.filters = [polar_filter]
@@ -101,17 +93,11 @@ class TalliesTestHarness(PyAPITestHarness):
         polar_tally2.scores = ['flux']
         polar_tally2.estimator = 'analog'
 
-        polar_filter2 = PolarFilter((5,))
         polar_tally3 = Tally()
-        polar_tally3.filters = [polar_filter2]
+        polar_tally3.filters = [polar_filter, mesh_filter]
         polar_tally3.scores = ['flux']
         polar_tally3.estimator = 'tracklength'
 
-        polar_tally4 = Tally()
-        polar_tally4.filters = [polar_filter2, mesh_filter]
-        polar_tally4.scores = ['flux']
-        polar_tally4.estimator = 'tracklength'
-
         universe_tally = Tally()
         universe_tally.filters = [UniverseFilter((1, 2, 3, 4, 6, 8))]
         universe_tally.scores = ['total']
@@ -173,10 +159,9 @@ class TalliesTestHarness(PyAPITestHarness):
         self._input_set.tallies = Tallies()
         self._input_set.tallies += (
             [azimuthal_tally1, azimuthal_tally2, azimuthal_tally3,
-             azimuthal_tally4, cellborn_tally, dg_tally, energy_tally,
-             energyout_tally, transfer_tally, material_tally, mu_tally1,
-             mu_tally2, mu_tally3, polar_tally1, polar_tally2, polar_tally3,
-             polar_tally4, universe_tally])
+             cellborn_tally, dg_tally, energy_tally, energyout_tally,
+             transfer_tally, material_tally, mu_tally1, mu_tally2,
+             polar_tally1, polar_tally2, polar_tally3, universe_tally])
         self._input_set.tallies += score_tallies
         self._input_set.tallies += flux_tallies
         self._input_set.tallies += (scatter_tally1, scatter_tally2)

From 1a921e7d08fc41b72bf1dd65cd17e922222b78b1 Mon Sep 17 00:00:00 2001
From: Adam Nelson <nelsonag@umich.edu>
Date: Sun, 13 Nov 2016 14:56:57 -0500
Subject: [PATCH 07/12] Fixed final failing test

---
 tests/test_mg_tallies/results_true.dat | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_mg_tallies/results_true.dat b/tests/test_mg_tallies/results_true.dat
index a765c0ba0..8ab726bbb 100644
--- a/tests/test_mg_tallies/results_true.dat
+++ b/tests/test_mg_tallies/results_true.dat
@@ -1 +1 @@
-3d2ce1b8bdd558fe9f8560e8bb91455e1fedd80a47349344399e22f5b0fac1391f1721e7f12a42bb0c458fd6c87ebf9fa38fe9539fcf69292b21cebf8e5f8989
\ No newline at end of file
+864328b2c4f3c4bfa9c80c756acbedac6fbdf3d502c03c2f2a3e7461b774729c47a55341c9515eff246e1bd445485f0d0a05a5712663ee499046afdbe0c77aef
\ No newline at end of file

From df7ae534138080dedfa3e351cfe068b64f54dd2a Mon Sep 17 00:00:00 2001
From: Adam Nelson <nelsonag@umich.edu>
Date: Sun, 13 Nov 2016 15:27:06 -0500
Subject: [PATCH 08/12] Cleaned up mgxs_library and made sure the mgxs-part-iv
 notebook was up to date

---
 .../pythonapi/examples/mgxs-part-iv.ipynb     | 217 ++++++++++--------
 openmc/mgxs_library.py                        | 140 +++++------
 2 files changed, 172 insertions(+), 185 deletions(-)

diff --git a/docs/source/pythonapi/examples/mgxs-part-iv.ipynb b/docs/source/pythonapi/examples/mgxs-part-iv.ipynb
index cd011d862..312bb8ee6 100644
--- a/docs/source/pythonapi/examples/mgxs-part-iv.ipynb
+++ b/docs/source/pythonapi/examples/mgxs-part-iv.ipynb
@@ -429,7 +429,7 @@
    "outputs": [
     {
      "data": {
-      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAACBjSFJN\nAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEX////pgJFyEhJNv8RV\nUZDeAAAAAWJLR0QAiAUdSAAAAAd0SU1FB+AKHxEkM8uGp70AAAWFSURBVGje7Zs7cttADIZ9CSvX\ncrP0iCxUqbBc8Ag6xR6BhV2EvYvwFD4CCx1ABT1jMdgndpegRQnOrCbjpPlGESISC4A/gd27e8H5\n83CX3b4+iKJrRHkS4vkghMPBonRYWGwtfgD2YN+dRDUOoh6lACw0Noi9w2fESuEoAR/uVuMolX03\n9oXGT7F3eFL2iEfhUX1f4cPdL/ishs+68ai+udE4xPhexbjX2FfjGNoPj/DPNX4Tsd+EODr8FvsV\ndf1Hd9P2VvCi4+s/aXvrf+upAD+1/9GV1mkOH5X9vV6THtfvACslcaUCbESL61drBPtdI8SrFMWr\nELsXCkuFDYW75gbiP7d9Cf7bAYI/aCwUShrBvh30+lWQkzVgZ/HD4OixNCgcQpJ3BxU/Ln91elKo\nM5VEE38QtJ+Yv6cQ9xjKNYayyl8TypP8DfJnQ2H/b/N3ye9P83cT33SQv/sQh9gV7zZ/0dNj5HQa\nC5vVzv9+/WFN2w8KVaZ2BwL1+pv4g0x1QRfjq0dB4Q3kT277oP6VNL6gKxNU9a8zK+WLbi/Wwpdi\nhbboKqyxFOulHMj6v4W/AXbmUeAxrv9J/CqEBXaRKsXaodD4nsYvkT/G6H1D4SR/iPy1Roj9JsQ5\ne18/7EUHv1+Fvx/Xj5V9Ugb5K8TW4TZEEdcvoz/up0VTe9qsVIppKVX6a7D6y9ZvwEKjrtQxPtv6\nfXII9vCxKOGaIeAIfEF8IvAG8ie3vRK9rRQl+PPpSctbhfpTUCpviH+kxsZgpT91+snoX1l49KK3\niUQvICRy5aUw6l8leoVwoo3Uv1rKreF/UFLY6d9QP4L9Wf2r7EP9GOSfcsjZ56f60kz+XmVPXv+R\nuP49ff0T/53Rv6n/7m2lvXT9Wqd/VUz8hvh5M/ED6ILmt4mfHYZSaePnTWpsf/SvqV9O6dLYYClL\nEetnoH/LBLFoBvrX189uTv8++kot5vTvQD4/9jP690g9P/4z/bvo/XVG/xYoZZx+8fr3MxAtsf7t\nUOkG2JqsTtCIpgCt/qX1226KqZS7gfzJbe+c9jLrtIZ8lXD+s4umlW6AKIVrlML2/cXjgPFjlJqI\nRC+Fj0bVJe+vSh56pSdR6YkQ1ygF10Wqf0FeLta/iKn9Mv1L24ti2e+7W4n1b3T/W+L+t9H9T/Sv\nVboUmqJJon1/hZq8LnzRDlDrX1u0xRT1+6vEpomMmyYkqi95vIH8yW1PN+122KkLcNLKi/WTF01z\n/cNASrWE/l3ev6T17zX909z9X27/euK/Rf3zWP+Waf9eEv37KkWJ+rfDl6ZglNDa+cEBhwYDvkoN\nP/rX69814NaI3imq0l7OYDy/qSdDGwr7r+Y3VbzoKZr6XX2lfxfOb87qXzr+b1j/Xlp/nP6dn98M\ncdH7cn7zjPObKsYWS3Eb9w8n85smHtqQuPuZ30T2dlIT6F9xFl+n8xslegL9a4c2KRr9W4rp/GYq\numiM9Nec/j2v/yj9u1h//hv9e93vc++f63/u+rPjL3f+5Lbn1j9m/eXWf+7zh/v8+2b9e/Hzn6s/\nuPqHrb8g71n6L3f+5Lbnvn8w33+4718/+5d47//c/gO7/5E7/nPbc/tv3P4fs//I7X9y+6/fqH+v\n6j9z+9/c/ju3/8+eP+TOn9z23PkXc/7Gnf9x5483q38Xzn+582fu/Js9fy8kb/6fO39y23P3n3S8\n/S/c/Tfc/T83uX/pgv1XE/9duP+Lu/+Mvf8td/znti8kb/8ld/9nx9t/Sjw/Ltr/yt1/+337f6/b\nf0zoB3nJ/ucVc/81d/83e/957vzJbc89/8A8f8E9/5HE78XnT/4H/cs5f8Q9/8Q9f8U+/5U7f3Lb\nc88fdrzzjyvm+cuf/Uu887/c88fs88954/8vO4SjPC+2QRIAAAAldEVYdGRhdGU6Y3JlYXRlADIw\nMTYtMTAtMzFUMTI6MzY6NTEtMDU6MDD2oqCYAAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE2LTEwLTMx\nVDEyOjM2OjUxLTA1OjAwh/8YJAAAAABJRU5ErkJggg==\n",
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAACBjSFJN\nAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEX////pgJFyEhJNv8RV\nUZDeAAAAAWJLR0QAiAUdSAAAAAd0SU1FB+ALDQ8YOYloD4IAAAWFSURBVGje7Zs7cttADIZ9CSvX\ncrP0iCxUqbBc8Ag6xR6BhV2EvYvwFD4CCx1ABT1jMdgndpegRQnOrCbjpPlGESISC4A/gd27e8H5\n83CX3b4+iKJrRHkS4vkghMPBonRYWGwtfgD2YN+dRDUOoh6lACw0Noi9w2fESuEoAR/uVuMolX03\n9oXGT7F3eFL2iEfhUX1f4cPdL/ishs+68ai+udE4xPhexbjX2FfjGNoPj/DPNX4Tsd+EODr8FvsV\ndf1Hd9P2VvCi4+s/aXvrf+upAD+1/9GV1mkOH5X9vV6THtfvACslcaUCbESL61drBPtdI8SrFMWr\nELsXCkuFDYW75gbiP7d9Cf7bAYI/aCwUShrBvh30+lWQkzVgZ/HD4OixNCgcQpJ3BxU/Ln91elKo\nM5VEE38QtJ+Yv6cQ9xjKNYayyl8TypP8DfJnQ2H/b/N3ye9P83cT33SQv/sQh9gV7zZ/0dNj5HQa\nC5vVzv9+/WFN2w8KVaZ2BwL1+pv4g0x1QRfjq0dB4Q3kT277oP6VNL6gKxNU9a8zK+WLbi/Wwpdi\nhbboKqyxFOulHMj6v4W/AXbmUeAxrv9J/CqEBXaRKsXaodD4nsYvkT/G6H1D4SR/iPy1Roj9JsQ5\ne18/7EUHv1+Fvx/Xj5V9Ugb5K8TW4TZEEdcvoz/up0VTe9qsVIppKVX6a7D6y9ZvwEKjrtQxPtv6\nfXII9vCxKOGaIeAIfEF8IvAG8ie3vRK9rRQl+PPpSctbhfpTUCpviH+kxsZgpT91+snoX1l49KK3\niUQvICRy5aUw6l8leoVwoo3Uv1rKreF/UFLY6d9QP4L9Wf2r7EP9GOSfcsjZ56f60kz+XmVPXv+R\nuP49ff0T/53Rv6n/7m2lvXT9Wqd/VUz8hvh5M/ED6ILmt4mfHYZSaePnTWpsf/SvqV9O6dLYYClL\nEetnoH/LBLFoBvrX189uTv8++kot5vTvQD4/9jP690g9P/4z/bvo/XVG/xYoZZx+8fr3MxAtsf7t\nUOkG2JqsTtCIpgCt/qX1226KqZS7gfzJbe+c9jLrtIZ8lXD+s4umlW6AKIVrlML2/cXjgPFjlJqI\nRC+Fj0bVJe+vSh56pSdR6YkQ1ygF10Wqf0FeLta/iKn9Mv1L24ti2e+7W4n1b3T/W+L+t9H9T/Sv\nVboUmqJJon1/hZq8LnzRDlDrX1u0xRT1+6vEpomMmyYkqi95vIH8yW1PN+122KkLcNLKi/WTF01z\n/cNASrWE/l3ev6T17zX909z9X27/euK/Rf3zWP+Waf9eEv37KkWJ+rfDl6ZglNDa+cEBhwYDvkoN\nP/rX69814NaI3imq0l7OYDy/qSdDGwr7r+Y3VbzoKZr6XX2lfxfOb87qXzr+b1j/Xlp/nP6dn98M\ncdH7cn7zjPObKsYWS3Eb9w8n85smHtqQuPuZ30T2dlIT6F9xFl+n8xslegL9a4c2KRr9W4rp/GYq\numiM9Nec/j2v/yj9u1h//hv9e93vc++f63/u+rPjL3f+5Lbn1j9m/eXWf+7zh/v8+2b9e/Hzn6s/\nuPqHrb8g71n6L3f+5Lbnvn8w33+4718/+5d47//c/gO7/5E7/nPbc/tv3P4fs//I7X9y+6/fqH+v\n6j9z+9/c/ju3/8+eP+TOn9z23PkXc/7Gnf9x5483q38Xzn+582fu/Js9fy8kb/6fO39y23P3n3S8\n/S/c/Tfc/T83uX/pgv1XE/9duP+Lu/+Mvf8td/znti8kb/8ld/9nx9t/Sjw/Ltr/yt1/+337f6/b\nf0zoB3nJ/ucVc/81d/83e/957vzJbc89/8A8f8E9/5HE78XnT/4H/cs5f8Q9/8Q9f8U+/5U7f3Lb\nc88fdrzzjyvm+cuf/Uu887/c88fs88954/8vO4SjPC+2QRIAAAAldEVYdGRhdGU6Y3JlYXRlADIw\nMTYtMTEtMTNUMTU6MjQ6NTctMDU6MDCpBKn7AAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE2LTExLTEz\nVDE1OjI0OjU3LTA1OjAw2FkRRwAAAABJRU5ErkJggg==\n",
       "text/plain": [
        "<IPython.core.display.Image object>"
       ]
@@ -575,7 +575,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "/home/romano/openmc/openmc/mgxs/library.py:370: RuntimeWarning: The P0 correction will be ignored since the scattering order 0 is greater than zero\n",
+      "/home/nelsonag/git/openmc/openmc/mgxs/library.py:398: RuntimeWarning: The P0 correction will be ignored since the scattering order 0 is greater than zero\n",
       "  warn(msg, RuntimeWarning)\n"
      ]
     }
@@ -731,9 +731,9 @@
       "         Copyright | 2011-2016 Massachusetts Institute of Technology\n",
       "           License | http://openmc.readthedocs.io/en/latest/license.html\n",
       "           Version | 0.8.0\n",
-      "          Git SHA1 | da5563eddb5f2c2d6b2c9839d518de40962b78f2\n",
-      "         Date/Time | 2016-10-31 12:36:52\n",
-      "    OpenMP Threads | 4\n",
+      "          Git SHA1 | 1a921e7d08fc41b72bf1dd65cd17e922222b78b1\n",
+      "         Date/Time | 2016-11-13 15:24:57\n",
+      "    OpenMP Threads | 8\n",
       "\n",
       " ===========================================================================\n",
       " ========================>     INITIALIZATION     <=========================\n",
@@ -743,12 +743,12 @@
       " Reading geometry XML file...\n",
       " Reading materials XML file...\n",
       " Reading cross sections XML file...\n",
-      " Reading U235 from /home/romano/openmc/scripts/nndc_hdf5/U235.h5\n",
-      " Reading U238 from /home/romano/openmc/scripts/nndc_hdf5/U238.h5\n",
-      " Reading O16 from /home/romano/openmc/scripts/nndc_hdf5/O16.h5\n",
-      " Reading Zr90 from /home/romano/openmc/scripts/nndc_hdf5/Zr90.h5\n",
-      " Reading H1 from /home/romano/openmc/scripts/nndc_hdf5/H1.h5\n",
-      " Reading B10 from /home/romano/openmc/scripts/nndc_hdf5/B10.h5\n",
+      " Reading U235 from /opt/xsdata/nndc/U235.h5\n",
+      " Reading U238 from /opt/xsdata/nndc/U238.h5\n",
+      " Reading O16 from /opt/xsdata/nndc/O16.h5\n",
+      " Reading Zr90 from /opt/xsdata/nndc/Zr90.h5\n",
+      " Reading H1 from /opt/xsdata/nndc/H1.h5\n",
+      " Reading B10 from /opt/xsdata/nndc/B10.h5\n",
       " Maximum neutron transport energy: 2.00000E+07 eV for U235\n",
       " Reading tallies XML file...\n",
       " Building neighboring cells lists for each surface...\n",
@@ -819,20 +819,20 @@
       "\n",
       " =======================>     TIMING STATISTICS     <=======================\n",
       "\n",
-      " Total time for initialization     =  5.3479E-01 seconds\n",
-      "   Reading cross sections          =  3.8403E-01 seconds\n",
-      " Total time in simulation          =  3.9455E+01 seconds\n",
-      "   Time in transport only          =  3.9331E+01 seconds\n",
-      "   Time in inactive batches        =  4.6776E+00 seconds\n",
-      "   Time in active batches          =  3.4778E+01 seconds\n",
-      "   Time synchronizing fission bank =  1.0514E-02 seconds\n",
-      "     Sampling source sites         =  7.2826E-03 seconds\n",
-      "     SEND/RECV source sites        =  3.1267E-03 seconds\n",
-      "   Time accumulating tallies       =  3.7644E-04 seconds\n",
-      " Total time for finalization       =  7.7600E-06 seconds\n",
-      " Total time elapsed                =  4.0029E+01 seconds\n",
-      " Calculation Rate (inactive)       =  10689.2 neutrons/second\n",
-      " Calculation Rate (active)         =  5750.82 neutrons/second\n",
+      " Total time for initialization     =  2.4361E-01 seconds\n",
+      "   Reading cross sections          =  1.6584E-01 seconds\n",
+      " Total time in simulation          =  7.3766E+00 seconds\n",
+      "   Time in transport only          =  7.3482E+00 seconds\n",
+      "   Time in inactive batches        =  8.5442E-01 seconds\n",
+      "   Time in active batches          =  6.5222E+00 seconds\n",
+      "   Time synchronizing fission bank =  4.9435E-03 seconds\n",
+      "     Sampling source sites         =  3.3815E-03 seconds\n",
+      "     SEND/RECV source sites        =  1.5249E-03 seconds\n",
+      "   Time accumulating tallies       =  9.5694E-05 seconds\n",
+      " Total time for finalization       =  2.9260E-06 seconds\n",
+      " Total time elapsed                =  7.6377E+00 seconds\n",
+      " Calculation Rate (inactive)       =  58519.0 neutrons/second\n",
+      " Calculation Rate (active)         =  30664.7 neutrons/second\n",
       "\n",
       " ============================>     RESULTS     <============================\n",
       "\n",
@@ -970,7 +970,20 @@
    "metadata": {
     "collapsed": false
    },
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/nelsonag/git/openmc/openmc/tallies.py:1944: RuntimeWarning: invalid value encountered in true_divide\n",
+      "  self_rel_err = data['self']['std. dev.'] / data['self']['mean']\n",
+      "/home/nelsonag/git/openmc/openmc/tallies.py:1945: RuntimeWarning: invalid value encountered in true_divide\n",
+      "  other_rel_err = data['other']['std. dev.'] / data['other']['mean']\n",
+      "/home/nelsonag/git/openmc/openmc/tallies.py:1946: RuntimeWarning: invalid value encountered in true_divide\n",
+      "  new_tally._mean = data['self']['mean'] / data['other']['mean']\n"
+     ]
+    }
+   ],
    "source": [
     "# Create a MGXS File which can then be written to disk\n",
     "mgxs_file = mgxs_lib.create_mg_library(xs_type='macro', xsdata_names=['fuel', 'zircaloy', 'water'])\n",
@@ -1107,9 +1120,9 @@
       "         Copyright | 2011-2016 Massachusetts Institute of Technology\n",
       "           License | http://openmc.readthedocs.io/en/latest/license.html\n",
       "           Version | 0.8.0\n",
-      "          Git SHA1 | da5563eddb5f2c2d6b2c9839d518de40962b78f2\n",
-      "         Date/Time | 2016-10-31 12:37:32\n",
-      "    OpenMP Threads | 4\n",
+      "          Git SHA1 | 1a921e7d08fc41b72bf1dd65cd17e922222b78b1\n",
+      "         Date/Time | 2016-11-13 15:25:05\n",
+      "    OpenMP Threads | 8\n",
       "\n",
       " ===========================================================================\n",
       " ========================>     INITIALIZATION     <=========================\n",
@@ -1133,56 +1146,56 @@
       "\n",
       "  Bat./Gen.      k            Average k         \n",
       "  =========   ========   ====================   \n",
-      "        1/1    1.00711                       \n",
-      "        2/1    1.01538                       \n",
-      "        3/1    1.01664                       \n",
-      "        4/1    1.03592                       \n",
-      "        5/1    1.00771                       \n",
-      "        6/1    1.00555                       \n",
-      "        7/1    1.02573                       \n",
-      "        8/1    1.04322                       \n",
-      "        9/1    1.02270                       \n",
-      "       10/1    1.02354                       \n",
-      "       11/1    1.02023                       \n",
-      "       12/1    1.03047    1.02535 +/- 0.00512\n",
-      "       13/1    1.04476    1.03182 +/- 0.00711\n",
-      "       14/1    1.02223    1.02942 +/- 0.00557\n",
-      "       15/1    1.02082    1.02770 +/- 0.00465\n",
-      "       16/1    1.01472    1.02554 +/- 0.00437\n",
-      "       17/1    1.02104    1.02489 +/- 0.00375\n",
-      "       18/1    1.04471    1.02737 +/- 0.00408\n",
-      "       19/1    1.02806    1.02745 +/- 0.00360\n",
-      "       20/1    1.02044    1.02675 +/- 0.00330\n",
-      "       21/1    1.02592    1.02667 +/- 0.00298\n",
-      "       22/1    1.02242    1.02632 +/- 0.00275\n",
-      "       23/1    0.99969    1.02427 +/- 0.00325\n",
-      "       24/1    1.02213    1.02412 +/- 0.00301\n",
-      "       25/1    1.02080    1.02390 +/- 0.00281\n",
-      "       26/1    1.01033    1.02305 +/- 0.00277\n",
-      "       27/1    1.02881    1.02339 +/- 0.00262\n",
-      "       28/1    1.01649    1.02300 +/- 0.00250\n",
-      "       29/1    1.03817    1.02380 +/- 0.00250\n",
-      "       30/1    1.00958    1.02309 +/- 0.00247\n",
-      "       31/1    1.01811    1.02285 +/- 0.00236\n",
-      "       32/1    1.02709    1.02305 +/- 0.00226\n",
-      "       33/1    1.01823    1.02284 +/- 0.00217\n",
-      "       34/1    1.01208    1.02239 +/- 0.00213\n",
-      "       35/1    1.01380    1.02204 +/- 0.00207\n",
-      "       36/1    1.02358    1.02210 +/- 0.00199\n",
-      "       37/1    1.03653    1.02264 +/- 0.00199\n",
-      "       38/1    1.03117    1.02294 +/- 0.00194\n",
-      "       39/1    1.00915    1.02247 +/- 0.00193\n",
-      "       40/1    1.03107    1.02275 +/- 0.00189\n",
-      "       41/1    1.02316    1.02277 +/- 0.00182\n",
-      "       42/1    1.02677    1.02289 +/- 0.00177\n",
-      "       43/1    0.99361    1.02200 +/- 0.00193\n",
-      "       44/1    1.04841    1.02278 +/- 0.00203\n",
-      "       45/1    0.99768    1.02206 +/- 0.00210\n",
-      "       46/1    1.02694    1.02220 +/- 0.00204\n",
-      "       47/1    1.03540    1.02256 +/- 0.00202\n",
-      "       48/1    1.03539    1.02289 +/- 0.00199\n",
-      "       49/1    1.02498    1.02295 +/- 0.00194\n",
-      "       50/1    1.00692    1.02255 +/- 0.00193\n",
+      "        1/1    0.98369                       \n",
+      "        2/1    1.01520                       \n",
+      "        3/1    1.03642                       \n",
+      "        4/1    1.02658                       \n",
+      "        5/1    1.03102                       \n",
+      "        6/1    1.05382                       \n",
+      "        7/1    1.01978                       \n",
+      "        8/1    1.01753                       \n",
+      "        9/1    1.02420                       \n",
+      "       10/1    0.99889                       \n",
+      "       11/1    1.04874                       \n",
+      "       12/1    1.01382    1.03128 +/- 0.01746\n",
+      "       13/1    1.03987    1.03414 +/- 0.01048\n",
+      "       14/1    1.02282    1.03131 +/- 0.00793\n",
+      "       15/1    1.03282    1.03162 +/- 0.00615\n",
+      "       16/1    0.99669    1.02579 +/- 0.00769\n",
+      "       17/1    1.00052    1.02218 +/- 0.00743\n",
+      "       18/1    1.01124    1.02082 +/- 0.00658\n",
+      "       19/1    1.00629    1.01920 +/- 0.00602\n",
+      "       20/1    1.05322    1.02260 +/- 0.00637\n",
+      "       21/1    1.00763    1.02124 +/- 0.00592\n",
+      "       22/1    1.01841    1.02101 +/- 0.00541\n",
+      "       23/1    1.03430    1.02203 +/- 0.00508\n",
+      "       24/1    1.03064    1.02264 +/- 0.00474\n",
+      "       25/1    1.03272    1.02331 +/- 0.00447\n",
+      "       26/1    1.01226    1.02262 +/- 0.00424\n",
+      "       27/1    1.00883    1.02181 +/- 0.00406\n",
+      "       28/1    1.02712    1.02211 +/- 0.00384\n",
+      "       29/1    1.03146    1.02260 +/- 0.00367\n",
+      "       30/1    1.02964    1.02295 +/- 0.00350\n",
+      "       31/1    0.99832    1.02178 +/- 0.00353\n",
+      "       32/1    1.03420    1.02234 +/- 0.00341\n",
+      "       33/1    1.01860    1.02218 +/- 0.00326\n",
+      "       34/1    1.03328    1.02264 +/- 0.00316\n",
+      "       35/1    1.01865    1.02248 +/- 0.00303\n",
+      "       36/1    1.02643    1.02264 +/- 0.00292\n",
+      "       37/1    1.01070    1.02219 +/- 0.00284\n",
+      "       38/1    1.01871    1.02207 +/- 0.00274\n",
+      "       39/1    0.98827    1.02090 +/- 0.00289\n",
+      "       40/1    1.01740    1.02079 +/- 0.00279\n",
+      "       41/1    1.02920    1.02106 +/- 0.00272\n",
+      "       42/1    1.02496    1.02118 +/- 0.00263\n",
+      "       43/1    1.04288    1.02184 +/- 0.00264\n",
+      "       44/1    1.03749    1.02230 +/- 0.00260\n",
+      "       45/1    1.04338    1.02290 +/- 0.00259\n",
+      "       46/1    1.03146    1.02314 +/- 0.00253\n",
+      "       47/1    1.04668    1.02377 +/- 0.00254\n",
+      "       48/1    1.02707    1.02386 +/- 0.00248\n",
+      "       49/1    1.02589    1.02391 +/- 0.00241\n",
+      "       50/1    1.02100    1.02384 +/- 0.00235\n",
       " Creating state point statepoint.50.h5...\n",
       "\n",
       " ===========================================================================\n",
@@ -1192,27 +1205,27 @@
       "\n",
       " =======================>     TIMING STATISTICS     <=======================\n",
       "\n",
-      " Total time for initialization     =  3.6445E-02 seconds\n",
-      "   Reading cross sections          =  4.9377E-03 seconds\n",
-      " Total time in simulation          =  3.0983E+01 seconds\n",
-      "   Time in transport only          =  3.0902E+01 seconds\n",
-      "   Time in inactive batches        =  3.2106E+00 seconds\n",
-      "   Time in active batches          =  2.7772E+01 seconds\n",
-      "   Time synchronizing fission bank =  9.7451E-03 seconds\n",
-      "     Sampling source sites         =  6.9236E-03 seconds\n",
-      "     SEND/RECV source sites        =  2.6796E-03 seconds\n",
-      "   Time accumulating tallies       =  3.2976E-04 seconds\n",
-      " Total time for finalization       =  7.4870E-06 seconds\n",
-      " Total time elapsed                =  3.1057E+01 seconds\n",
-      " Calculation Rate (inactive)       =  15573.4 neutrons/second\n",
-      " Calculation Rate (active)         =  7201.41 neutrons/second\n",
+      " Total time for initialization     =  2.5229E-02 seconds\n",
+      "   Reading cross sections          =  1.1001E-02 seconds\n",
+      " Total time in simulation          =  7.3074E+00 seconds\n",
+      "   Time in transport only          =  7.2846E+00 seconds\n",
+      "   Time in inactive batches        =  6.2995E-01 seconds\n",
+      "   Time in active batches          =  6.6774E+00 seconds\n",
+      "   Time synchronizing fission bank =  4.8069E-03 seconds\n",
+      "     Sampling source sites         =  3.3693E-03 seconds\n",
+      "     SEND/RECV source sites        =  1.3618E-03 seconds\n",
+      "   Time accumulating tallies       =  8.0542E-05 seconds\n",
+      " Total time for finalization       =  2.9260E-06 seconds\n",
+      " Total time elapsed                =  7.3508E+00 seconds\n",
+      " Calculation Rate (inactive)       =  79370.9 neutrons/second\n",
+      " Calculation Rate (active)         =  29951.7 neutrons/second\n",
       "\n",
       " ============================>     RESULTS     <============================\n",
       "\n",
-      " k-effective (Collision)     =  1.02346 +/-  0.00203\n",
-      " k-effective (Track-length)  =  1.02255 +/-  0.00193\n",
-      " k-effective (Absorption)    =  1.02775 +/-  0.00139\n",
-      " Combined k-effective        =  1.02594 +/-  0.00120\n",
+      " k-effective (Collision)     =  1.02315 +/-  0.00205\n",
+      " k-effective (Track-length)  =  1.02384 +/-  0.00235\n",
+      " k-effective (Absorption)    =  1.02372 +/-  0.00194\n",
+      " Combined k-effective        =  1.02369 +/-  0.00173\n",
       " Leakage Fraction            =  0.00000 +/-  0.00000\n",
       "\n"
      ]
@@ -1294,8 +1307,8 @@
      "output_type": "stream",
      "text": [
       "Continuous-Energy keff = 1.024739\n",
-      "Multi-Group keff = 1.025941\n",
-      "bias [pcm]: -120.2\n"
+      "Multi-Group keff = 1.023689\n",
+      "bias [pcm]: 105.0\n"
      ]
     }
    ],
@@ -1392,7 +1405,7 @@
     {
      "data": {
       "text/plain": [
-       "<matplotlib.text.Text at 0x7f4c89d6c4a8>"
+       "<matplotlib.text.Text at 0x7f4ad8f44eb8>"
       ]
      },
      "execution_count": 40,
@@ -1401,9 +1414,9 @@
     },
     {
      "data": {
-      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAADDCAYAAACS2+oqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXd4FdXW/78rCEoLHaQrzS4IylXBggVFQbkgCAe4the8\nYntf/FlBUa/1esFyFSsiXj2iUgQbguBVIArSu7Q0AgQIJSShJuv3x0zOmTlzzqwJSTgn4/o8T57M\nXrNm7X32rFmzZ8/svYmZoSiKolR8kuJdAEVRFKVs0ICuKIriEzSgK4qi+AQN6IqiKD5BA7qiKIpP\n0ICuKIriExIuoBPRaiK6PN7l+DNDRN8R0ZBSHP82EY0syzL9GSGiIiJq5bLft9eK+uBxwsziH4AA\ngN8BHACQBeBbAF28HCvYnQDg2dLaieef+RsOA8g1/w4AWBbvcnko92gARyxlzgXw/+JdrhKUeQ+A\n+QAuLsHxPwG48wSUMw3AIQB1I+TLABQBaOHRTiGAVhY/K9G1AqAygKcArDfPcaZ57V4b73MZ5Xyq\nD5bBn9hCJ6IRAMYCeA5AQwAtALwFoJd07J+Il5k52fyrycwXlHUGRFSprG0CmGQpczIz/6sc8ihr\nJjFzMoD6AP4L4Mv4FicqDCAVwMBiARGdC6Cquc8rVMpyTIFxnQ4GUAfA6QBeB3BD1MzKx8ck1AfL\nEuFukgzjztnHRacKgNdgtNy3AngVQGVz3xUwWgUjAGSbOreb+4bCuNMdgnG3m27KUwFcZbkbfg5g\noqmzCkBHS95FMFswZtrWijHz2AhgN4CvADQ25S3NY5Oi3TkBtIZxovYB2AngM5ffH7PlZMnnbwDS\nTVtPWPYTgMcAbAKwC8AkALUjjr3TPPa/pvxvMFqAuwCMKq4vAI0A5AOoY7Hf0cyzUoyWxsdSK8Kt\nLsxznQ1gP4AVAM4uyXmwnMO7AWyA0eJ5U2gdfWxJnwWjFVvPTNcG8LVZzhxzu4m57zkAxwAUmL70\nhik/E8AsU38dgH4W+zcAWGPqZwIY4bEVlgrgCQCLLLJXADxulrdFtNYagNsAzIv0b3i4VqKU4RrT\nHxp7KOsj5vk7CKMb9iyzbHthXHO9ovmGS5nvB7DZPA//9Ho+1QdL74NSC/0SACebFRCLUQA6Azgf\nQHtze5Rl/6kAagJoAuB/ALxFRLWY+X0An8I44cnMfHMM+70ABAHUMivnLcu+mK0dIroKwAsAbgHQ\nGEAGjIApHgvgHwB+YObaAJoB+LeLrhe6AGgL4yJ7iojOMOUPALgJwGUw6mcvgHERx14O44RfR0Rn\nwfj9A2H8plrmcWDmbBgXQX/LsYNhOH9hKcoetS6IqDuArgDaMHMtM9+cyIM9nAcAuBFAJxj+09+0\n7QoRVYERTHJg1BtgBKMPATSH8SRZANNfmHkUgHkA7jP97QEiqgbjQvoERmtrAIBxRHSmae8DAEPZ\naI2dC2CuVC4LvwGoSURnEFESgFvNfKRWt8MvS3CtWLkawEJm3u5BdwCAHjCCURKAGQBmAmgAw0c/\nJaK2JShzbxiNiY4AbiaiOz2UwQ31QY8+KAX0egB2M3ORi04AwDPMnMPMOQCeAWB9mXEEwD+YuZCZ\nvweQB+CMKHZiMZ+Zf2DjdvUfGDeOYtwujgCA8cy8gpmPwmgdXUJELTzkeRRASyJqysxHmDlF0H+Y\niPYQ0V7z/wTLPgbwtGlnJYxWRHtz390ARjLzdrOMzwK4xQwAxceOZuaDzHwYhkPOYOZfmfkYjP5R\nKx/DrHvTxkAYdRaLWyPKfWoJ6uIojBv12UREzPyHeVOJxMt5eJGZDzBzJoybUgepzDAulLsA3FLs\nn8y8h5mnMfNhZs4H8CKMG2IsegJIZeaP2WAFjG6Kfub+IwDOIaKazLyfmZe72IrGf2Bc8NfCaHlt\nK+HxpaE+gB3FCSKqY57nfUR0MEL3dWbeZvrYxQCqM/PLzHyMmX8C8A0s3UceeMmsr60wnt7djlUf\nLEMflAJ6DoD6lgATjSYw7njFpJuykI2IG0IBgBpCvlZ2WLYLAJwilMdarvTihFm5OQCaejj2YRh1\ns4iIVhHRHQBARI8T0QEiyiUia0v6FWauy8x1zP93RNizOpn197cEMM105D0A1sJw0kYW/a0RvynT\n8psOwt4imQ7gLCJqCaA7gH3MvNjld34eUe4dUXSi1oV5ob8Jo/WRTUTvEFG08+rlPMSqn5hlhvE+\nZzWAC4t3EFFVInqXiNKIaB+AnwHUJqJYN/6WAC4urn8i2gvj4i+u/74wWm7pRPQTEV3sUq5ofGLa\nux3GzbbcMP2y2DebwajjxsX7mXkvM9eB0QqtEnF4TB8zSYe36yaavch4EIn6YBn6oBQYf4XxBUdv\nF50ss1DWAnptiZTkBVE0CgBUs6Std/dt1nIRUXUYTxxbYfQtItaxzLyTmYcxc1MAf4fxCNSKmV/k\n8Mub4aUsO2DcCHuYjlzs1NUjHpOtdbQdxiNn8W+qav6m4nIfBvAFjFb6YLi3zj0Rqy7MfW8y84UA\nzobx1PVwFBNu56E05doD4wnnaSIqdv6HYHRtXWQ+nhe3jIovpkh/y4TxbsJa/8nMfJ+ZxxJm7g2j\n62E6jLotSRkzYPRR9wAwNYpKPmL7r8OckFdNi29uBTAHwEVEFC2YRgYXq+1tMLoLrLSAcZ17LbP1\n+BYo5ZOJ+qB3H3QN6MycC+MlwFtEdLN59zmJiHoQ0Uum2iQAo4ioPhHVB/AkvAeSbBgvfUqC1RmX\nAQgQURIRXQ/jJWwxnwG4g4jOJ6KTYfSh/cbMmcy8G4aDDjaPvRPGixcjA6JbiKj47r0PxksTt24n\nr+WN5F0ALxQ/+hFRAyK6yeXYyQB6EdHFRFQZwNNRbP4HRouwF8ogoMeqCyK6kIg6E9FJMF6mHUL0\nOop5HkpbNmbeAKOv91FTVNMsSy4R1YWzfiL97RsA7YhosOnXlc3fdaa5HSCiZDbeQRyA8fKrpNwJ\n48VlZDcHACwH0Me8rtrAeHyPRYmuFWaeDaPr4CvzPFU2z9UlcL85LARQQESPmHVyJYxugc9KUOaH\niag2ETUH8CCc/dUlQn3Quw+KXRfMPBbGVyqjYLy5zQAwHOEXpc8BWAyguH94MYDn3UxatsfD6B/a\nQ0RTo+yXjv9fGC8V98Lop5tmKfccGDeXqTCC9+kwXjgUMxTG2/3dMN5UL7DsuwjAQiLKNX/nA8yc\n5lKmR8xH3VzzsXdnjPJGpl+HcdedRUT7AaTAeKkc9VhmXgvjC4LPYbQ6cmGck8MWnRQYTr20FA5r\nzTdWXSQDeB/GVwGpMOrxFYch+Ty41Y8X/gVgqNmYeA1G63E3jLr8LkL3dQD9iCiHiF5j5jwYXVMD\nYNTnNgAvIdwlMQRAqvnoPAzGo7AXQr+BmVOZeWm0fTC+0DgKo1txAowumqh2cHzXyl9hBIxPYFwj\nW2BcJ9YXfpE+dhRGY+AGGPX4JoAhzLzRY5kBw6eXAFgK40OGD4VyRkN90KBEPkjMpe31UOKF+ei4\nD8Zb/nSLfA6AT5n5eC4kRTluiKgIhj9uiXdZ/owk3NB/xR0i6mk+7lYHMAbAyohgfhGAC2C04hVF\n+ROhAb3icTOMx7KtMPr9Q4+ORPQRjG9aHzTf5CvKiUYf+eOIdrkoiqL4BG2hK4qi+ISTysOo+Qnh\nazBuGOOZ+eUoOvpooJQrzFzaya0cqG8riUAs3y7zLhcyRnFugDGXxDYY0+4OYOb1EXps/2S0L4wR\nr1a8TKHygKzykvwbv3ukm6hzQ+Anu2BeX+Aye5mfD44Q7Yz8bqyoM+jGD1z3z+PLRBsZ30WZYeH5\nvsDIcJn73SAPYPxy/G2iTngcngteXG16FKX0vkDLcJm7rvpRNDOfupd5QC+Zb/9ikYyC8XVvMXPk\nzP4xWtZ5cqas89X1okqbm1Y6ZNv7jkDjKWE/3bKztUMnkvYN5ZkRll3aRdTBJ7KjdG79i0O2oe9T\naDfl2VC6behLy9h8Om6oXJ7qsgp+lsv8wYRBtvS4vvMwfIr9Ov4IkYPM7VyEeniVOsX07fLocukM\nYCMzp5vftE6C8SJPUSo66ttKQlMeAb0p7HNBbEXJ5oFQlERFfVtJaMqlD907fS3bO2HMkmvFbV6p\nYiKPicJy+XHop2C0SdoiSIvIi4scshXB9RBZLpc5df9C1/35vN9DPkucsqIi4L/h/DP2/SrbWVhZ\n1tktq3jqctkXRYmLgH3hMu8MrnaoFKxNR8G6DIc8flhnkN4DYLYlvUY+fIUHv8YKWeXnPaLKgbwo\n9VZUhAPB8CBHzm3k1IlgT7KHfrfdHnRmyI6yu1GU66ywCLuD4e64k6LO8xXB7x76U06WVbBFLvPC\nYJotXVTIDtnOKLPjWn17h2NeNTvlEdCzYEzIU0wzhCf2icDa/xyEc1SrY2rjKHgYjd1BruxugfdF\nnVe+iZLXaXZZ+4B8E/qitlzm028scN2/1UMfek7tGLMUXxnOv8UNx0Q7Cw96qOOy6kPfHEPJUmcN\nA9760MuBEvi2tc98NowZdIvx8GDc3kOdT64r61wh96HXjNKHDgA1A+GFjXZ56EOv66EPPf1ND33o\nN8mOUj9KHzoA1A9cE9o+3UMfeso+D/XspQ/9iFzmvwS+jSI7zZZeh6tcbRT3oceiPLpcfgfQhoha\nkjEB/AAYE+YrSkVHfVtJaMq8hc7MhUR0H4wRi8Wfdq0r63wU5USjvq0kOuXSh87MM1GyVYkUpUKg\nvq0kMvF9KTra8inlKgLOi/i08p37ZRvZHjpmf5M/R77xUXm5yJHBkbb06uAanBuwv+CqDnkKlcLW\n8uLqO1DLdf8wvCfayNgUJe5kw1iS2uQ6+kG00/wueRbesVNGiTpYIZ+rpis3OWQFwWxUC4Tl8/+4\nxqGTcJzTNby9LwOobUnnd3XqR+LhnXfR4zeIOtxOtnM013l9TDrIGJAb9vdxDYeJdh7e6Zi51kH+\nbNn3z62+StTJ5WSH7CBXtck/Xfo/op3qt+0SdfJX1Bd16g6O8SrFwsB8+7TwSYcZt+b/ZpN9VN39\nO3QJHfqvKIriEzSgK4qi+AQN6IqiKD5BA7qiKIpP0ICuKIriEzSgK4qi+AQN6IqiKD5BA7qiKIpP\niOvAohajwzOm5Qe3oXrAPoNaRgd5QF6Dm+UZ9r5Af1FnL9cRdXrje1s6iCACEZODPYeHRTtV68kz\n4B1i94FFbfGiaKOwp3MQR7CIEeg5JJQmLhTt8HZ5MEha39NEnam1B4k6WRktncKc+thrkTdt5xx8\n5LAjapQzoyyDdRYQ0MWS/tzD8dd5GDA3T1b55YzOos4V9JtDdlLVIKokh337zkJ5ysEHh8uD3Wiy\n7G+bPpX9bfog5+Rrv2AbLrdUSu+O3zt0ItkJedKxyZf2FXVexBOiTvVdh+yC3Em4c9cAu+wr95lN\nqzV2z0Nb6IqiKD5BA7qiKIpP0ICuKIriEzSgK4qi+AQN6IqiKD5BA7qiKIpP0ICuKIriE+L6HXpG\nUqoltRM5g1MjNFqJNnZcebqoQ3Plb1+TPhNV0CBg/+b9EHLwf7DLKtE9op2H6o8RdXKPvey6f8w8\neXHnGh13OmTHakzBsHrh72rz3pK/+e0+/CtRZzZ6iTqcJuf11NWPOWSr663BuS3CYxTeo7tFO3Hn\nQst35NvYnvbwHXrhJrmu8HSRqHJprtxm46POvDiPwTnh8QrJsRbvtrBn8imiziJ0E3VOGnSpqNN7\nl3NhloJcoPeu8OIY++vK5dlS6SxRZzg+FHW2Y6So83zW83bBnkpAlv278+8C7vVTHxdhlst+baEr\niqL4BA3oiqIoPkEDuqIoik/QgK4oiuITNKAriqL4BA3oiqIoPkEDuqIoik/QgK4oiuIT4jqwCP97\nXXh7fQ5w5nW23SljOoomkt6TBzzw6/IgjcLzRRUsoAts6dm0B9fSRJtsE7cV7dw5Th7FdGR4Fdf9\n53Vb5bofABbiYofs22oHcGOtV0Lpl+99QLTTh6aKOvfyRlGn8139RJ3PIhYMAYA8fItVuDGU3vVz\nC9FOvKlz2rbQ9pEGe1HFkt57XRPx+KQlHvx6quzXleUxMUCqMy/aD9DrYXnGmvqimRbDd4s6zd/O\nFHVyUE/UWdbgTIcsLXk/ljUILwzTBNscOpH8Zc8KUYcXy/V8YfdrRR1HtK3klA084h4brqIqAP4V\nc7+20BVFUXyCBnRFURSfoAFdURTFJ2hAVxRF8Qka0BVFUXyCBnRFURSfoAFdURTFJ2hAVxRF8Qnl\nMrCIiNIA7AdQBOAoM3eOqjjTcj/ZnwSk2e8vXZouFfMqfIhEnV/wF1Gnc95iUacr23UyOIiubB8I\ns5kGiHZeuec+UechvOmusES+Fz/WabRDtharsBfnhdIvwakTCe+R85pZ50pRpwfmijpvYIFDVogC\n1MK+sKCpaKbc8OzbbjYulHX4bnk1IqyXz8vYr+UVtEbQW05hMAgEwr7dor2c1zfLrxJ1evKPog7e\n87DKUg/ndb82h9Ehc0coTc3llcqQ4aFN+5uscnP3maLOJ5372NIpmzJxaefJNtngDPflrPJPcY93\n5TVStAjAlcy8t5zsK0q8UN9WEpby6nKhcrStKPFEfVtJWMrLMRnAbCL6nYiGllMeihIP1LeVhKW8\nuly6MPN2ImoAw/nXMfP8cspLUU4k6ttKwlIuAZ2Zt5v/dxHRNACdATidfmtfy0HOFxi8TH4xFAzK\nL0XXI0fUST0kz25X+ZSgLZ2SkuLQ+ZXSRTt5RXtEnSAF3RVSRRNY+4dzRsasFPtsd0EI+QDgfDmv\nldWzRZ29HvLagw0OWX5KxO/IjlLHm9YCm9eJ9kuLV9/O62dpvBfa/ZjT64j5BP/wUBh5MkEsXeKs\nT0deUXzN4dv7HCoOlgV3iDq5kl8DwCJZhQ87r9eUxYDxAGVA9TzkJV+uwEo5NiAo55UC+7W3ISVK\nHMiZ5JRtXAdsWg8AWC70qZR5QCeiagCSmDmPiKoD6A7gmajKzaaEt/cHgVr2L0boAjmgBwJevnJ5\nXdTpnCdHyFNqOKd2DQTssqM0Q7Szu0ieijRAzrxsLBks2ljZ6byo8rMDYXkgynS1kfBeOa+6dRqJ\nOj085DUmylcuAFA3EJ6eNH1zF9EO2pZ9b2JJfLvGl++Hto98Ng1VBv41lC5YLk+fG+jkoUDr5fOy\n44x2cl4xfM3m2y/LeSUHThV1ekp+DQB5cl7RvnIBGIHeYTk195DXcjkvHJVjDAIe8sJkh+TSQHNb\nelyG+1dyHU4hzDo19tTa5dFCbwRgGhGxaf9TZp5VDvkoyolGfVtJaMo8oDNzKoAOZW1XUeKN+raS\n6MR3xaL1Cy2JTcD2hbbdXFseEFTpc7lbptWA8aLO/BpdRZ2Td9pXLuFcBu8cYpPVbiSvXLIoSR6L\nMpO7ue5f20keMPLvXOcApmMHp2BObvjdRf/kc0U76+r0EXXOIrn/ekDRRFFn2Yu3OYUr0pGeZulm\ncT65JhznVFoT2t6ZtBUNLem/dJogHp+DF0SdD8+8V9S5EEtEHVwbZUWeHQxMsPj2G7KZ6w7IA8eK\njsmr/6wfdrqoc85EZxdp0hogKTnc312UIue16t02ok77xzaJOvyRnNfAiLFHnM4YOMM+eHLMJPcX\nCK1R23W/fk+rKIriEzSgK4qi+AQN6IqiKD5BA7qiKIpP0ICuKIriEzSgK4qi+AQN6IqiKD5BA7qi\nKIpPiOvAor8Xhuft2BjcgLaBGrb9E3LbijYI8sQ5L+AJUSeL5GVwejT8wJbemzwLYxp2t8newP2i\nnXxUE3XexTDX/dN/GSja6H75dIdse9U/0Dg5PJfUEaos2ilkedBEx4XywKJJte4QdYY+8b5DNje4\nC1cFxoXS1942R7SD5rJKefJ1Ua/Q9pdciH5Fr4bSnSrJg336s/vKNQDw6EPyaJ9hY2WdDrOWOWT5\nwULsDYTPe52XDot2TpkhrxD0xhh5xuF7l3wo6rx/2yCHbGHlVOQHwoOShi77VLRz/qUeBg19Japg\nV8Pqos6w2961pbcGF2BywD4v0bJH3OcpaiBMzaMtdEVRFJ+gAV1RFMUnaEBXFEXxCRrQFUVRfIIG\ndEVRFJ+gAV1RFMUnaEBXFEXxCRrQFUVRfEJcBxZ9k9QztF2QxPjDkgaAe5LfEW2MpcfljP77jajC\ng+WFYJdmXmJLB5GKAOwyXGNfgSQaXSbKK81QE/dBGpwv34uP7nf+pkkFjAH7Xwulq9SSB4O8jM9E\nnYHt5GWEKteW87rqXecgph2LGFcdsAwAeflk0U68OT9pZWg7n77GP5LCA40e5ZfE4y/CSlFn7Rh5\nZZ8anCfq1KGDDll1CqKOZUHn9MflRcD/yfKgunvxgahDL8gDi+56IOiQVV3HCPz8ayjNr8mrmeF6\n+Tqa2qCHqNMH34o6/cm+8tcCykQX2maTTb9bGDBY1X23ttAVRVF8ggZ0RVEUn6ABXVEUxSdoQFcU\nRfEJGtAVRVF8ggZ0RVEUn6ABXVEUxSdoQFcURfEJcR1Y1JkXhbYzeTOaW9IA8No8edDQhsuniDqD\nr7hJ1EnLlAdpjMNGW7oAO/BYhKzp7LminXvwtqhzc+EprvsH9fhStPEQxjhkG6rtxoJa9UPpbqPk\n1Yim3i6vCvV5G7mO+9SV8/oi568OWUqNTHDAsgTRPfKAs3izKzdcx4UHa6LAkn6+1kjx+Fu4lqiz\nH7LPnjkuXdThWc7zwlsZ/MWQsGB6fYdOJN3wk5zXl7IPTJtyvajTJ3WmQ5ZUH0hqEfZVPlPOiyaJ\nKujzr+9FnaK75LwGHcy3C/Z8jnFZt9pl9whGOrnv1ha6oiiKT9CAriiK4hM0oCuKovgEDeiKoig+\nQQO6oiiKT9CAriiK4hM0oCuKovgEDeiKoig+4bgHFhHReAA9AWQz8/mmrA6AzwG0BJAGoD8z749l\nY+qsQeHEKsLvswK2/fd2f0UsR3UqEHXO49Wizu17Joo6B/fUtaWD2YzApodtspFtR4l2crieqFMt\n/4jr/ptrTRdtLGbnKIQt2IBaaBdKdzv2m2hnTpsuos6Ah2eIOs/E9IQwd9HPDtk+OowraUtYsExe\nHQkXDJd1YlAWvt06OVzefVV3orYlfRrSxDK8SfKqVvfjTVFnwfALRJ0u1y5zyOhrgHqFB+mswnmi\nnV9xqajTce5zok6z/ltFnWdOe9ghW1l/LTaednYoPXqfHD/4PVEFK8e1EXWaQS7zJ7DHt5S6mbi0\nqf06/uesR11tnI5awMux95emhT4BwHURsscA/MjMZwCYC8DD+nCKknCobysVkuMO6Mw8H8DeCPHN\nAIqbuhMB9D5e+4oSL9S3lYpKWfehN2TmbABg5h0AGpaxfUWJF+rbSsJT3i9F5VmdFKVior6tJBxl\nPdtiNhE1YuZsIjoVwE5X7X/0DW8XFTl2b9jtfFkTycnk/vIQAL7lPFGnME+evTCYb7+GU5YCkdf1\nmlNXiXaS+YCoM+mg+/5F1VJFG/lczSFLS9luSwfXimawOrhL1Nm5TrazykMI/Cp42CFbnHLULsgI\nOg/cshbY4qEQx0+JfDuj7yOhbY7w7SzaLWZWGfL5nYKjos4O3iPqpO9wygzfDrOscRSlCDZ6+Pgg\nuFFUwZbgPlFnZZHTcTNTsux5HZLzYg/lyQjK12tdFIo6K5FpS29IcZ6bvfjBITu0NhWH16UBAHJR\n2TWP0gZ0Mv+KmQHgdhjvYW8D4P4pxpOWqW9/CgLd7G+B23XPgoSXr1xu5OWizlN7+ok6gT13R0gY\ngV5kk6xpK38N0IDlADkgd47r/kO15KlT98WYgvWCQPgrl8DK2aKdHwMNRJ1rlm0QdTbKs5Cid+Bk\nUf7AikBUHRsXlPrhs1S+3WLKP0Pb+4IzUTsQnhK2KaWJmZ8L5804kr4epqvdwHVFnS4bo0+xG+gV\n3k5ud6popwqfK+oE5n0l6vweqC3qFBWdHVV+fiAsDzz0rWiH24oqWBmoKeo0g/wJVw00d8guDdhl\n8x3v4u1cglp4Nyn6bwdK0eVCREEAKQDaEVEGEd0B4CUA1xLRHwCuNtOKUqFQ31YqKsfdQmfmWM2k\na47XpqIkAurbSkUlrisWVe6UG9ouSjuIJEsaAHoLPTYAcDXPlzNq9oyoMna7/GhKhfZ+MmoUBLWx\nX/sv9JQfergviTp0u/OdgpWhy+V8Zra/wiHbjWycB8t7h5fc8wGAzoXRu0FsVJJ/0+hCuZ9xPTm7\nkg4RIY+qh9JvtB8q2nlA1ChfnqJnQ9vzKQtdKbwaVxLLdd4X34g6n1aSHxICX8vdjbghSnlODQJt\nw77dc5jsb5nvOrsUHLwt//Yn+GtRZ9bfnV+NBjcxAr98F0rzDjmvmUny77pu9WZRh36RfbspXWxL\n16G9aEr2lY46wP29Yeso3TZWdOi/oiiKT9CAriiK4hM0oCuKovgEDeiKoig+QQO6oiiKT9CAriiK\n4hM0oCuKovgEDeiKoig+Ia4Di56u93Roe0WN9Whfb41t/zv0d9HGq0VTRJ1v/yKXZfgt8sxR/IF9\nEAAvZHDBEJts99c1RDsBfCrqzBpeyXX/u+P+JtpoHjEZEAAUIQmFlvs4z3XPBwCOXCXPK0KXyxNF\nFb0u53X7g3Mdst38IyZzeJDmoqSLRDvAeA865ceTHB5YlMvf4wfuEUqvn9ZRPL7gNrmuBv1TVAFd\nKPt1UaYzL85hcGbYt4e/M1a08/aIEaLOF2PlycIOkGxn7TunOWRZwTysDYSvv8mV5Doc7Rx75yB/\npjxobiDkyf02Y5wtvR8zMRXX22S7Ud/VxkGc4rpfW+iKoig+QQO6oiiKT9CAriiK4hM0oCuKovgE\nDeiKoig+QQO6oiiKT9CAriiK4hM0oCuKoviEuA4sWoJOoe1MHMYxSxoApk2UFwM+don8E3iKvHIJ\nnpbvbZOfvtGWXlQtC1UCTW2y23MniHam1LpF1Ekd18h1/0e4Q7Tx+7bODhnv/Rxjt90aSi+7qr1o\n57xn5RVb+Cm5jmvly4tj582JsiD1mixsmWMdAbJJtBNvNnzXIZxYvhY7aofT/f76sXj815f1EnUO\n1Y++CLjFnlywAAAJu0lEQVSVZ/CoqPPUiFccMloPUGp4UNKNY+UVlOBh7e42JJ+79lgh6swi52LK\nq2gjKlN41efR7d+RCzReHjRU/eRjos42kldOqwr7gvYHcdgh645ZrjbOQ1PXIUzaQlcURfEJGtAV\nRVF8ggZ0RVEUn6ABXVEUxSdoQFcURfEJGtAVRVF8ggZ0RVEUn6ABXVEUxSfEdWBRCi4JbR/ELmRa\n0gDAteWP/u9vKy/b8u+N8solr46+R9QZgbds6aMIol/EaIpbsuV75Lo68u9qVVjour82TRdt7Gzs\nXP1kSu1j6Nv43lC6HvJFO9lP1RZ16uyX67h9rR9EnQVnXOYUrjsKnHHYImgj2ok3bW5YGdo+sC8D\nNS3pxXSheHzL+mmiztv8mqjzZJo8iAljowwKCwaBQNi3W6G1bKeR7NcdeY2oczedI+qswnkOWXXk\noS4sKyItlQe78V75en2RHxJ1lr43RtRxjIdbl4G0ZV3tdt7pCje6X+2ehbbQFUVRfIIGdEVRFJ+g\nAV1RFMUnaEBXFEXxCRrQFUVRfIIGdEVRFJ+gAV1RFMUnaEBXFEXxCcc9sIiIxgPoCSCbmc83ZaMB\nDAWw01R7gplnxrKx441W4cSShti/u5VdoQNEptJfRZ0DbWuKOiMgDwzgAfbBM5zO4BlDbLKcSdVF\nO2dfLg/mWYEzXfcP9FA5k2igQ7aYNuMohQeJDE+VBwS90WqkqNO41nZRZ8HH14g69w5xrp6zoc4y\ntGu6M5R+y8vSOKWgLHz7bKwNbW9FFppZ0jMWDhDL8HrnYaLOfSd9KOrcc/tEUWfF+DMcsgzkYgWe\nCaWbYLdop9J2eWWfvzcbK+q88+BqUSfrjboOWT4dwXW0PpSey11EOzl1eoo6Ty79l6jz5TDZTr/1\n39oF3xDQM2IwlttyRAAghLLStNAnAHCuAwWMZeaO5l9Mh1eUBEZ9W6mQHHdAZ+b5APZG2SWP/1WU\nBEZ9W6molEcf+n1EtJyIPiAieRVbRak4qG8rCU1ZT841DsCzzMxE9ByAsQDuiqn9Yd/wdlGUiXTk\nbjsczNop6qSy3Cf3PXJFndXpbEun7AYAuywvKPcj1swWVZARdC9PGqWLNg7zKQ7ZlhR75sFd7NCJ\nZE1DD32avE/Uwa9BUWVDpWUO2faUtAhJ1ShHbjT/yo0S+faivq+GtjnStzfLk0Yt3hQ5k5MTZrk+\ng5tFlai+tjzloC1dG0fl8mz2cH43LpXt/CHbmRo84pD9vsB+7WXxLtFOHg7J5UmTyzNvfZaog20R\ndpYtcOpEe712ZC1wdB0AYPnP7lmUaUBnttXg+wC+dj3gzinh7SVBoFPEyy4PL0WrXrFF1Dmdq4k6\nPTBP1Dl/RuTNgxFoaX8KzwnIVVr//cOizopAsuv+5dRStJHHNaLKLwyEX4oGUn8R7axtda6o05jl\nl6IfH5NfZrYLRL8w2gUuCG3PHiK/CAeae9DxTkl9u/OU/wttbw0uQLNA+AXd0kXyS9ELO/9X1Plk\niFyfgdZDRJ1YvtbDIm+CAtHOvYs8nN+L5NbMXA92+gTuiyGvEtpeyw1EOzmoJ+q8uVQuz2UdPxN1\n3lgfxU7PCNlbThVYvrPocBkw65PYPX+l7XIhWPoViehUy74+AOSmnaIkJurbSoWjNJ8tBgFcCaAe\nEWUAGA2gGxF1AFAEIA3A3WVQRkU5oahvKxWV4w7ozBztOWRCKcqiKAmB+rZSUYnrikX4xtIXtIOA\n7RF9Q8fkF3b3X/6mqDPyUXnQUOeXF4k6Z79jf+lWNBk4dou9jI1Gyy9X2/60QtTpRO4vj1qx/Lbr\n+aufdwqzg5g4PhyvPp9zq2hn/uprRR3IC9HgnL8tFnXeWvGwU5gZxOyV1hj7jFMnwZjxlKWffHUR\nlq63pAfJXz8+OPE9OZNOsp2kO+UXsPSF8zrj34IYeVK4zm/qP0m0w/Pl8vznHLlPn9vJPcFNpkf5\nqnRxEPdVD5eZDsnx49b+H8nlOSiXp98N34g6jo88cgB8ESGL/torTLTvASzo0H9FURSfoAFdURTF\nJ2hAVxRF8Qka0BVFUXxC4gT0vLWyToKx7o94l+A4yK949YwtFbDMVnZVwPJvXRfvEpSczApWzwfL\nvryJE9DzK54Dra+QAb3i1TO2VMAyW9ldAcufVQHLXNFuQofKvryJE9AVRVGUUhHX79A7tg1vb94M\ntG4bodBQtnEqmsn5NJXtNIA8NwpV6hgh2Ayq1Nom6thYzquF9DEpgNOEOSaa4KDrfsBev8VsTrXX\ncztpxnwABc45vpzUkVVaQ55T5+QoVbO5EtDaKvdQyUvlOaDKlY5NwtubqwKtLWmc7MGAPMUIhDVQ\nDOT1VkBRzt3mKkBri7yVhxPc0cP1Wi1JDjkF8hQsQJS5LjdXBlpb5CRfZjjdQ0V39FCHaONBp749\nuTkbaB15DoVpntq0AGa57Cdm+eP78oCI4pOx8qeBmeMyf7n6tlLexPLtuAV0RVEUpWzRPnRFURSf\noAFdURTFJyREQCei64loPRFtIKJH410eLxBRGhGtIKJlRCTP7BUHiGg8EWUT0UqLrA4RzSKiP4jo\nh0RaSi1GeUcT0VYiWmr+XR/PMpYE9evyoaL5NXDifDvuAZ2IkgC8CWOV9XMADCQiL+/v400RgCuZ\n+QJm7hzvwsQg2ur1jwH4kZnPADAXwOMnvFSxiVZeABjLzB3Nv5knulDHg/p1uVLR/Bo4Qb4d94AO\noDOAjcyczsxHAUwCcHOcy+QFQmLUX0xirF5/M4CJ5vZEAL1PaKFciFFewLJyUAVC/bqcqGh+DZw4\n306EE9cUQKYlvdWUJToMYDYR/U5EQ+NdmBLQkJmzAYCZd8DT1/5x5z4iWk5EHyTao7QL6tcnloro\n10AZ+3YiBPSKShdm7gjgBgD3ElHXeBfoOEn071bHAWjFzB0A7AAwNs7l8Tvq1yeOMvftRAjoWQBa\nWNLNTFlCw2wsc2+uBj8NxiN2RSCbiBoBoYWPd8a5PK4w8y4OD5Z4H8BF8SxPCVC/PrFUKL8Gyse3\nEyGg/w6gDRG1JKIqAAYAmBHnMrlCRNWIqIa5XR1AdyTuKvC21eth1O3t5vZtAKaf6AIJ2MprXpzF\n9EHi1nMk6tflS0Xza+AE+HZ81xQFwMyFRHQfjCkKkgCMZ+ZEnzatEYBp5hDvkwB8ysxuUyzEhRir\n178E4EsiuhNAOoD+8SuhnRjl7UZEHWB8fZEG4O64FbAEqF+XHxXNr4ET59s69F9RFMUnJEKXi6Io\nilIGaEBXFEXxCRrQFUVRfIIGdEVRFJ+gAV1RFMUnaEBXFEXxCRrQFUVRfIIGdEVRFJ/w/wGvJaRK\nz6tq0gAAAABJRU5ErkJggg==\n",
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAEPCAYAAADf8cexAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XecFdX5x/HPAyIICliwBCuCiTFqBKOJBey9xZgEsP5i\n7LEQo8aKJTFRoxgTTdQYFQuIorFEwRoVuy4WFAsIYqMoShWF5fz+OHPl7uzd3XmW3Xt3me/79bov\n2LnPnXPmzpxznznTLISAiIiI5FebSldAREREKkvJgIiISM4pGRAREck5JQMiIiI5p2RAREQk55QM\niIiI5JySARERkZxTMiAiIpJzSgZERERybplJBszsf2b2RKXrIS2LmY0xs4ebYb63mtl7TT1fKS8z\nO9/MFmeMPcLMFpvZus1dr2WBme2cfF/bNPF8N0zmO7Ap55t3jUoGzKyHmV1rZhPN7Cszm5V0uieZ\nWYemrmRRuRub2eA6GmMAMjXqZYWZTU4aRanXg5WuX3Mys371LPvtRaEheTW1imxvRR1h4VVtZp+b\n2QNmttVSzHfbpG2t2JT1XRpmdnjRcpb8QTGzD5P372tkMbW2DzM708z2zxJbH4sOM7OHzWyGmX1j\nZtPMbLSZHWVmyzeyzhWX7HyVanvVZrZRUWhz3e++IvfRT3YCipd3gZm9nbSdRq9PMzvbzPZtyrp6\nLef9gJntBdwJLACGAuOA5YHtgEuB7wPHNmEdi30fGAw8AUxJvbdrM5XZkgVgLPAXwFLvfVL+6lTE\nlcDLqWmTi/6/I83TcRxB7e+8nG4BRgNtge8CJwCPm9mWIYS3GzG/7YDzgOuBuU1Wy6bxFTAQeLZ4\nopn1A7oT+6KmdBaxj7s3NX0oMCyE8E1DM0h2iv4D7AY8A1wGTANWAfoBVwNbAUc1XbXLKgAfAr+n\njr4nhPCYma2Q5ftyFRzCxOaYb9bigXnA0cTl7gocQPxdWh/4v0bO9xxim75/6avYOK5kwMzWB4YD\nk4CdQgjTi97+h5mdC+zdZLUrUQXq6NhDCIuasdyW7OMQwrBKVwJiBxhCaOqOuSFjQgh31/Vmc20X\nIYTq5pivwyshhG9HQMzsOWJHcixwSiPmV8nEpiEPAj83s5NCCMWjMQOJieBq5ahEiE91y/oDdCVx\nB+WkEMLfU+8NMbMNaWAHxszaAm1CCAvdlS2PWQ31Pc31g12hRKBgYWq5rzGzF4BDzezUEMLMSlVs\nqYQQMr+AfwDVwNYZ49sC5wITiNn7JOAPwPKpuMnAfcC2wAvEPYGJwKFFMYcTh2Wrk38L/++bvP8/\n4PGi+H5JzM+Bs4lZ7FfAo8CGJcr/d4n615hnMq0bcAMwNZnfq8BhqZhC2X1T09dLph9WNG0N4Mak\nfguIWfV/gHUzfL+TgPsyxN0EzAG+k8x7DjCduLdiqVgj/piMS5ZvKvBPoGsd62w34KWk7icl73UA\nrgJmALOTMr+TLPt5ScyOyd/7l6jvwOS9Orezou/4wAaWfQzwcGraycCbxAx/JvAi8POi9zsn9Z+c\nLNc04l74pkUxtwLvpea7IjCkaF2OB04p0SYWA1cABybf8wLgDWCXDOtyw+TzJ6Wmd06m35+avjlw\nM/B+sj4/Je79r1wUcxGl29Z3Uu3vZWA+8DlwW/H7ScxGwN0saRtTkrhOnn4mVWZ18j0tAnYveq9d\nUo9TSLUDfO1vMFBd9Hf6e1hM0jcQR4MW00DbBNYmJg0POJa1ULffJtvnBGAhsFkz9Ts3EfuBDZJt\ney7wMXBuxvo+AbzeQMzOSbnbeLYRYA9iu/0yqePbwIUl2sDAVHm7EEdh5gFfJOVslIr5Q/LZ9Ykj\nPV8msdcD7TMs9y3AzBLTr0i2mz6p6Wckdfqc2HZeAg4o0R+kt7nrimK6J+trKkv6isNL1KHefq2h\nl/cwwT7A+yGEFzLG3wAcBowgDmVvTRyC2xj4WVFcAHoRh+ZuSBb8V8CNZvZyCGE88BSxgz6RuEIL\nQ6Hji+ZRyu+JX/RlQBfiyrkV+Emq/FLSxxI7EBOEDYG/EX8sfg7cZGZdQgh/yzDPtLuJ38dVwAfA\n6sQ9hnWpfSiklHZmtmqJ6fPCkr30QDw/ZDTwPHAqseH8ltjpXFv0ueuI6+zfwF+JncWJwA/NbNuw\nZI84AN8Dbk8+fx3wTvLezcBBxMb2ArGT+i9F30kI4QkzmwIcTO3h2IOBCRm3s5VKLP/MkLQOaq/D\n44g/2MOSf1cANiNum3cmYdcD+xLX8dvEPc/tiOvpjaL5hqL5WrKM2yaffx3YE7jCzNYKIZyRquMO\nxG3nGmJHfAow0szWDSHMyrDcaRsk/36Rmr47cVsq/JD8ADgmWZbtkpgRQE/gF8BviB0kxA4FMxtM\nPIRwe7JsqxM7nq3MbIsQwlwzaw88TNzOriQmUGsTv8fOxA6qsSYTt9sBxG0YYK9kvsOTuqQ19tDQ\nIcTv6gXiNg1xx6Qwzyzz3ZPYyd/WiPJ/BbQntqmvgZnN1O8U+oRRwHPAacQf4QvMrG0I4fwM82hb\nou0tCCEUr+viNtLgNmJmmxL7g1eIQ+dfE38b6j0J0cx2Bx4A3iXugHYibhfPJNvoR0X1CcBIYt93\nBrAl8Xufmny2MepqfyclZd1KPJw+kNjO9wwhPBxCqDazQ4g7hGOI2x5J3TCzNYk/6t8QfyM+J277\nN5pZpxDCNUlcln6tfo7MdSVixnJ3xvjNkvh/pqZfSvxx7lc0bVIyrTiDXI2YOV5aNO1nFI0GlMhU\nS40MjAPaFk0/MZnH91PllxoZSM/z5OSz/VOZ3TPALJLsNim7Vj1JZejE5GQx8Nus6yE1v0nUzCaL\n9+pOL4q7MZl2VurzrwAvFv29XfL5X6bidk2m90+VXU1qbxbYIon9S2r6v5P484qm/ZGYLa+UWu/f\n0MAeStH6TWfU1RTtuQFPUzQyQBxKr2pg3rOBKxqIuQV4N7VtLgZ+l4obSdzDW7doe1mcLHdxPQvf\n29ENlFvYKzoTWJX4w7w9ca+9Gtg3FV9rb4eYbNUY4SN2ijVGA5LpPYh75aempm+aLNfvkr/7JPXa\nt776O7fvwshAb+B4YpLSPnnvDuDRom0xPTLQYPtLptUYGUimzaF0f1CoT0MjA5cncZumprdL1lnh\ntUqJun1RPD15r0n7nWRaoU8Ykoq9n9jvrtLAMj5B6X7n30UxO1PUr2fZRog7KtUU9Qn1tIGBRdPe\nII5sFPclP0zmdX3RtMIo2DWped4LfJJhm7wlWUeFddgDOD0p5+US8e1Tfy9H3Ht/KDX9K4pGA4qm\n30TcKeySmj4C+AxoV7Te6u3XGnp5ribonPw7J2P8XsQMbEhq+uXEoej0uQVvhRC+PUEohPAZcU+z\nh6OOpfw71Dy++3RSfmPmuycwNYQwvDAhmfdVxCHifs75fUX84dvBzLo2oj4Q95h2Ju7pF167EjPE\ntGtTfz9Nze/hIGKH+5iZrVp4EU9SnEsc2i82KYTwaGraHsT1/o/U9L9R+7j0UOIhhYOKpvXHt1d1\nAbWXfWo98V8C65rZFvXEzAJ+nGTlWe1JXJdXp6ZfQVyePVLTR4UQvh35CSGMJe49Z90u/0A8DDMV\neJK4Z39yCKHGCUghhK8L/zez9sn6fIG4LnpnKKcwgjcytU18Sjz0UNgmCqMJezbTFUUjgI7APskV\nD/vQuD3vcij0lekTMfcirrPCa3KJz94Vah9zbup+p1h6e/07cQ92lwyfnUTNvmdX4s5eXbJsI4WY\nn2YoHwAzWxvYBLghhPDt71MI4VXgcWr/1gRK94VrZNx2u7BkHU4A/kwcuT4wHZhqf12JJxyOIUPb\nS0Ybf0pMVJZLtb+HgZWJCQ9k69fq5TlMMDv5d6WM8YVsdELxxBDCNDP7Mnm/WKkh8S+IC7w0Piwx\nTxo53/WAUteWjyd2rullqlcI4RszO4N4CGWamT1PHOoaGkKYBmBmnYlDPgXfhBCKh6I+CyE8kaG4\nBSGEz1PT0t9vL+LGOp3aAnEvtNikEnGF9Z5+b0I6MITwjpm9RNxTvTGZPBB4PoTwfqmFKGFcCOHx\njLEAfyIO0b9i8T4BDwO3hRCeL4o5jTiS8ZGZvUw8gW1oCGFyPfNdD/gohPBVavr4oveLpbdLiA06\n63b5D+IhphWIHfJviElHDUnHcT7xEEC3orcCsVNrSE/isG6p9RFI+oUQz/D+K3FY9HAze4p4Tsmt\nxR10Y4UQPjOzR4nbR6ekTnct7XyXRj1ts7C86cs0x7DkR/Z0Sg99Ty4xrUn7nSKLqb1e33XMc17G\nvgfIvI3czpJDxJcRz/G6mzgiXdchkEJd3y3x3nhgJzNrF2qeiJn+vSn+Xfi0gUWZC+xP3AbXJo6q\nrU7cuavBzPYjHhrfnHj4pyDLCZBrEn9vjydeLZRW3Cdn6dfqlTkZCCHMMbNPiMODWRT2AutagWl1\nnZ29tGc5Z5lvXXVsSxwi9dalvvnVDAzhr8k10gcQj+9eCJxpZjuGEF4jHrc/vOgj/wN2yliPYlnO\nfm9DPI43kNLLOiP1d62Nvx51fSdDgSvN7DvEjvXHxI2/WYQQ3jKz7xL3LPcg7vmeYGbnhhD+mMQM\nN7MniVn5rsTk4Awz27/ESEiBdztd2u393aIk6L9xJ4K/mNn/ku2mYCRxePYS4nkM84jD1Q+S7T4j\nbYhtID2yUVC8JzbIzG4gdpS7EfcyzzCzH4cQ6hutyapwzsJaxGHWupKMzO1vKdXVNt8mrscfsOQc\nE5Jk/HEAMzu0jnmWalNN3u/Uo1mvKmloGwkhfGVm2xFHnPYmbncDiD9udW2Djanz0rS/RcVJUJKk\njicm6AcVTd8RuIe4zo8ljuItJF5OWnzOXF0K7fNm4jkHpbwG2fq1hnhPIHwAOMrMtg4Nn9w1mbgw\nvVhyYhlmtjpx7/MDZ9nQfDea+IJYp7T1WHLyEMRlKpUMbZz8W1imL1hyDWqx9UsVHkKYRDycUrjk\n6DXisbPDiJ34Lam6NpeJxL3MZ4uHt5w+IK73Daj53W1UOpxhxKH0AcRh4G+IQ8LNJoQwPyljhJm1\nIw7DnWtmfy4cUgohfEo8ue8aM+tGXCdnEfdUSpkMbJdc/1zcoae3jeZyEXBk8u9+8O2oQF/gzBDC\nJYVAM/teic/X1bYmEn9M3m9gZCTOJIRxxPN0/ph06k8Rr8m+MPOS1O0e4vDu1sAv64lztb8SsvYz\ndbXNh4g/NgdT+nCd12Saod8httMe1By1K7TTZtteG9pGkhGAx5PXqckl6+ebWd8QwlMlZjk5+fe7\nJd77HjAtNOPlmSGEj5MRj7PMrHcIoSp560Bi8r1H8aFqMzum1GxKTJuafL5NltHPLP1afbx3ILyU\neOLTv5If9Ros3h3tpOTPB1lymVqxU4kL/l9n2RC/mFIb+9KaSDxG/G1ylNwNap1U3IPAmmb2y6K4\ntsSTEucQj91CbEjVxI642PHUPLt2heQM22KTknm1BwghvB1CeLzoNbaRy5jFCGKCeF76DTNra2ZZ\nhpVHE9dReu/+REps8Mnx0YeAQ4md56gSx0ybjJmtkip/IXFPrg3xyoy2ZrZSKmYGcegwva6KPUg8\n1ppe7kHEbeGhpax6vZLh6euBvc1sk2RyoQNIt/NB1F4XhTPA021rZBI7uFS5he/TzDqbWbqcccln\n6/veMgvxLPVjiYc96rs5S6b2V495ZOhj6mqbIYQPiYeZ9jSzUsO74Ot7m7TfSflNib+/AR5z1C+T\nLNtIun0mCiNdJbejEK8UGAccUdx2zWxz4kjNA0tZ9Sz+Srzy4fdF0wonN387MmNmPYhXT6TV2uaS\nH/B7gF+Y2cbpD5jZakX/r7dfy7IArpGBEML7Fu8HPRwYb2bFdyDchni5y41J7OtmdjNwtJmtTNxg\ntybu7d4dQniyVBkNeJX4BZ+RnIzxNfBYcrLh0vgXcXhntJmNIJ6tegi1j3NfR7ws6yYz25Ill/j8\nhHjy1jyAEMJsM7sTOCkZvp1I3ADSN0fZiHiy3gjgLeJw7IHE40BZ9yi6m9nBJabPDSGkL9mrVwjh\nKTO7Fvi9mf2QODS3MKnnQcRjfXXe4CeZR5WZjQROSTbW54knOPUqhJT42FDi8d9AvJyoOT1u8ZLG\n54iHRDYhdpb3hhAWJHvTk5L19waxke5GPFHnpDrmCbHRPgVcYmY9WXJp4d7AZckPRHO7kvgDcQbx\nzPEvzexZ4mGnFYj3sNiDeKlhejj0lWTan5JlXwj8J4TwnsVLCy9MRq3uIx4z7UE8jPI34olsuxJH\ntu4kHt9uRxxCX0gD20wDatQzhHBLXYFFMVnbX11eAXYxs0HE72xSCOFFV63jTtD6wFVm1p+YvExP\n6rBtUp/xdX66pqbudwq+BvZI+unniSc47gn8scT5RY1VvP7q20ZGJjEXmNmPicnzB8Tj5scn/69x\nB8qU3xF/9J8zs38Tz9c4kXh5bFOMStUrOaflZuDXZtYzhDAhqc9JxN+VYcTDW8cTR8o3Sc3iFWA3\nMzuFuOMxMYTwMvHckr7Ai2Z2PXGbWYV4OeT2xO8HGujXsi6E+0X8sfwncWP7inji01PAcSSXOiRx\nbYide+GmQ5OJw5jtUvN7P6l0qctXHktN+xVxQ/qGmjcdqhHLkstsDkx9fr1kevqGHacQTyqZT0xc\ntqij/NWIycM0ltz849ASdV+VuKc9h3gJyNXEYb1vy05W6lXES01mEzfcZ9N1rmc9TErmV+r1flHc\njcS7haU/P5h4/Cs9/Ujita1zk3X7KnAxsEZD6yx5r/imQ7OIP/Q9iVnyaSXi2yXf0RekbkhVz7KX\nXL8l4p4GRhf9fUyyfqcn6/pd4iWOhcuzlicO/45Nln0WsaEemZrvLcA7qWmdiIc8PiJu728TO+vi\nmLZJvS8vUdcpwLUZ2l41cGId7w8ldvKFSxm7EzvamSy5WdBayTzOTH32XOKJjQupfdOhA4ltfHby\nepOYfPRI3u9BbBfvEROo6cAjlLgM2NHPHJ7Uo3cDcbW2RTK0v7raADH5fSLZ/r+9XI6MlxYWzceI\nOz+PJG3ha2K/8TDw6+JtnSX90qA65tVk/U5RnzCbmLCMSuI/wXfTodcaiElfWtjgNkLck7+HJTeJ\nm5Js0xuUaAPpmw7tTGzvc4l9yUigVyrmouSznVPTj0xv83Us0y3A53W815PYdq5LzfcdYl8zjriT\neRHxZNPiz36PeM5JYZsrnkc34rkVk4n9ysfEEdjDi2Lq7deyvCyZkUizSkYaqoCDQ+oWpsmQ5yfE\nDv3oStRPJE/M7EbgZyGEzg0GSy4sM48wlpajxHkQEEdeqol7l2k/Je75DG3OeomISGnupxaKZHC6\nmfUhDnstIh6L3J04BP5xIcjiI3c3Jx5KqgohjKlAXUVEck/JgDSH54gnC51DPJFnCvHY7MWpuOOI\nVxCMpfGP/hSRxtExYvmWzhkQERHJuRY/MpBc6rU7S86kFJHG6UA8e3x0aLpLx5qc2rxIk8rU7lt8\nMkDsFFrqA0lEWqODibf2banU5kWaXr3tviLJQHJXrt8Rb5jwGvGa6ZfqCJ8c/7mVeClmsUHUfigi\nxHsieQxwxgMn+g6vXHW4/4q5k866rvQbYwfBFrWX+/iLr3CXcc2Y37ri99z+Hlf82OB/iNbUMeuX\nfuP6QXBU7eXeZTv/DcYe/U+pm4DVo6FHl5TiPQL3ZB0f+HQQrFV7uX84oqE7gqeqM/4TXjvkj1D6\nYTjNztHuJwPQ61bomLrx2qRBsEGpNk+8itvD+3ifxmwDx/g2gm6H1H1vqlmDLqbLkLNqTZ/x9/SN\nUuu35omTXfEAR3KD7wPObf/At+q+QeegITBkUIkiFvvK+N+mP/Z9AFjk/Ilsj/8u7q+zWcnpjw56\nhF2G7Fprer+SF2TV7YvxfRh0yHBooN2XPRlIbql5OfFe1C8Sf9FHm9lGofSdBJNhwu9R+6mPXUtM\ng/iYb48sT3JN6e7b2nv2zvqwxyKr1FGvdl1Lvte9d607RDdsqm/ZV+ld1XBQkeVD+kZbGUwtdYtx\noFNX6Fm7viv39v4KAGMbsc69vMnACnV8oG1XWKF2fVfs7btrc1jS3Ms+9O5s97F+HTeGFVPL3bZL\n7WnfFuKs1AoNh9Qs2xkPsLZvI1i+d92X/bfpuhLL9y7Rntbs6Sqjfe+OrniA9VjV9wHnuWi96/lh\n77Ii9C7xNA1vMvBhPd9tXRZmu5Pvtzq4nt0WTWetktPbd2nPmr1rv/ddfOvvsyUPm6233VfiPgOD\niJeYDQ0hvE281/h84p0FRWTZpHYv0oKVNRlInqTUh6KHYIR4OcOjxPtsi8gyRu1epOUr98jAasTB\ntmmp6dNY8sAFEVm2qN2LtHAt5WoCo8EjrIOo/VTRdZupOi3cev0rXYPK6JfT5e7iX+4Zw55gxrAn\nakwLs+Y3VY2aSv3tftKgeI5Asfbes/6WHSv036fSVaiIAbtVugaV8f0B/vOtHh72BY8M+7LGtG9m\nzc702XInA58R70+/Rmr66tTea0gZQqNO9FsWrdeIqx+WBf1yutxd/cvdbcCOdBuwY41poWoSz/Q5\npqlq5dG4dr/BkLpPFsyhjgNymgzsXukaVMYmA37g/sxuA1ZmtwEr15j2WVU/9u1zVYOfLethghDC\nQuLjYHcuTLP44O2dqf9Z1SLSSqndi7R8lThMcAVws5m9wpJLjDoCN1WgLiJSHmr3Ii1Y2ZOBEMII\nM1sNuJA4bPgqsHsIYUa56yIi5aF2L9KyVeQEwhDCNcA1lShbRCpD7V6k5WopVxM07GiDtTLeYuyf\nJ/rmPa0RT2583ne7s73PeNxdxNm3n+2K78Q8dxnVG/puqzaVLg0HFTmaOm6pXI8pE+q4A2EddrfR\n7jLWObLu276WcsXIc9xl8Jpvu+r++gRX/Jh3dnHFd57iu3tkxR0HbOSIP805/16+9dP3jUZsZ8G3\nnd124K/dZeC7ASHLB/8tc399yK2u+Ktv891L6q99jnLFA5y85fWu+P3/8rC7jDd36uGK/w8HuMs4\nJVzpiu/24lxXfNVn2a5KqMQdCEVERKQFUTIgIiKSc0oGREREck7JgIiISM4pGRAREck5JQMiIiI5\np2RAREQk55QMiIiI5JySARERkZxTMiAiIpJzSgZERERyTsmAiIhIzrWaBxWtefQklu/dMVPslB/6\nHnTTbf8p7vqM4Beu+C/Cyu4yDuAhV/wf3E9qgRVWnemKXxB8DyrqxZ9c8QDV+/genmSh2l1G+NRX\nxuSfre8u4+6uB7viP56yniu++0a+Bxu1n/sRs12fqLA7A6zseJjQK8/45v+7bV3hz3/+Y9/8gSce\n3ssVf8vd/gcVbYLvAVTjJv3IXQa3LnaFn3Cfbz/zin2Pc8UDhFd8dbIN/fu+m2z3vit+ztDH3GU8\nbLu54rfc6mVX/MyqlTLFaWRAREQk55QMiIiI5JySARERkZxTMiAiIpJzSgZERERyTsmAiIhIzikZ\nEBERyTklAyIiIjmnZEBERCTnlAyIiIjknJIBERGRnGs1zyaYuuUnQLZnE0AP37x32MBdH3vcdz/8\nNsPcRdBtoO+ZCW3Nf3/vU1e73BU/e9ElrvjLn17kigdYsfd0V/zcq33PGQDY7fj/uOIfYV93GWGy\nr17n7fx7V/x1dowrvq21c8VX3G4GPSx7/Hu+Zw3YUF91FvzS91wOgGcOcTxbAWg3YHN3GftzgCt+\n9XU/cJdxsvO5J+33PcEVP93WcMUDPBl8z4qYPmEfdxk/t/tc8euwmruMlcIcV/xGX/h+F+YuyDZ/\njQyIiIjknJIBERGRnFMyICIiknNKBkRERHJOyYCIiEjOKRkQERHJOSUDIiIiOadkQEREJOeUDIiI\niOSckgEREZGcUzIgIiKSc0oGREREcq7VPKiIgT+B1XtnCn328mxxBW2u8z1MBCD81fcQmurN3EXw\njG3hip8QernL+NU1vicofXP88q74TXd8wxUP8AK+B5BccsJJ7jIOtLtd8SeE99xlbHXkz13xwxjo\nip/x5Lqu+M4TP3PFV9zMAJ0cbfMuXzvuvuEEV/yT4VBXPEC/3/niw0nj3GXs8udHXfEbdXzHXcY5\np/seaFZ16cau+N53jnfFAyzYyxdf3da/7+vt51c/1v9bclGXc13xG67s2247d1gPuL3BOI0MiIiI\n5JySARERkZxTMiAiIpJzSgZERERyTsmAiIhIzikZEBERyTklAyIiIjmnZEBERCTnlAyIiIjknJIB\nERGRnFMyICIiknOt59kEz7aBDtlyl227V7lmXX2quavzFFu74rea+7K7jO2C7zMTrb+7jMuO+40r\n/lT+7ivgFX+++fs+g13xf8YXDxBm+uo1auUd3GXsyeOu+Kt4xldAd184s5zxlbY20DN7eJfu01yz\nn4LvWR53WTdXPACb+voWO6TaXcSsZ333z991m3vdZXDpYlf4c/ZrV/waP/OtO4Du5nzWxh6N2Pft\n5Ft/7U73fU8AR/EDV/xP7R5X/F42JVNcWUcGzGywmS1Ovd4qZx1EpLzU7kVavkqMDIwDdgYKKdei\nCtRBRMpL7V6kBatEMrAohDCjAuWKSOWo3Yu0YJU4gbCXmX1sZhPN7FYzW6cCdRCR8lK7F2nByp0M\nPA8cAewOHAtsADxlZp3KXA8RKR+1e5EWrqyHCUIIo4v+HGdmLwIfAL8Abqz3w9MGQdsuNad1HhBf\nIlLT/cPggeE1Js2f82VFqtLodn/tIOjUtea0HfrDjmrzIqXMH/YA84c9UGPaqFlfZfpsRS8tDCHM\nMrN3yXIB0RpDoEPv5q+UyLJg3wHxVaTjW1XM3n/LClVoiczt/pgh0EttXiSrjgP2oeOAfWpM26Nq\nCtf12bXBz1b0pkNmtiKwIfBpJeshIuWjdi/S8pT7PgOXmVlfM1vPzLYB7iFeYjSsnPUQkfJRuxdp\n+cp9mGBt4HZgVWAGMAb4cQjh8zLXQ0TKR+1epIUr9wmEOvNHJGfU7kVaPj2oSEREJOdaz4OKJr8B\nLMwUGrr6HiLU9g7/wyV69L/BFT9mxe3cZbSf7nsASdc1Gj5jNO3FNlu54keFHV3xb/U5zhUP8LfZ\nvocn/aJeMDN0AAATnklEQVSz70EfAONXPtAVv7GNd5fRf/HNrvixfzrcV8BdvvBWl/r/G1gxe3jo\n53uoTNs7giv+wV/OccUD2EW+MsInvjYP0GZrXxmP3rCfu4ywju+7fXu3P7vi131/uiseYNEDvu+q\n+g53EbS7x/fdLj7Hv/5+tLrz4VRdGw4pNteybbetrXsQERGRJqZkQEREJOeUDIiIiOSckgEREZGc\nUzIgIiKSc0oGREREck7JgIiISM4pGRAREck5JQMiIiI5p2RAREQk55QMiIiI5FyreTbBz156jW69\nsz3x9MbZvVzzNnz3nwa4mLNc8R9bd3cZe67+L1f8VZzoLmMeHV3x13K0K/7ep/wPrNut772u+G+s\nnbuM6uC7h3jvF/zPJhje5f9c8Ueddb0rftfDH3PF8+Y3sIfvIxW1rUH37PfEn335Gq7Z73jWf13x\nT9HXFQ+wzZvPuuJXOnORu4wt+z3piv97P38/sfWzr7viZ9jqrvgnN/Q9IwXAhvr67SEnn+Au47TZ\nV7viJ/zB38/vx3BXfHu+ccX3qsrWhjQyICIiknNKBkRERHJOyYCIiEjOKRkQERHJOSUDIiIiOadk\nQEREJOeUDIiIiOSckgEREZGcUzIgIiKSc0oGREREck7JgIiISM4pGRAREcm5VvOgoqfbbM/ybTbJ\nFHtc53+65n2Fnemv0P8ecIWHQ7I/cKWg6sOf+D6wS5W7jG1v9j28w75T7YoP8/z55sJZvu9q+S6+\nOgFcwjBX/ICN7nKX0a6rr147Xet7eBKXtPfFL7+8L77SngywYvaH0fx1jO8hWidynSt+OL4HSQGs\n+D3ntjlhsbuMwbaLK35qWNNdBtu86goftqmv3V/8pr9/3K7a912ddpq/Lxpzha9e2500xV3GHhzh\nij/206Gu+Kq5A7kgQ5xGBkRERHJOyYCIiEjOKRkQERHJOSUDIiIiOadkQEREJOeUDIiIiOSckgER\nEZGcUzIgIiKSc0oGREREck7JgIiISM4pGRAREcm5VvNsgk3Cm6wc5mWKvfJp37MG3u070l2fQ/rt\n54qf/OEG7jKu4T1XfPdHHneXcRz/cMXvX93BFX/wnne64gFO5XJX/I7nOO/pD9x9RPZ73gPc0dO3\nvgEOXMVXrxGf/9RXwHG+Z3DQeaovvtLeWAQszBx+ynJbuGa/2iLf9z3gjvtc8QAM8oVfEk52F3HW\nlAdd8U+t63zmCcDlvm25zd98sz/7PV97BAi3+up01mWD3WX86cQsd/Vf4qq2/r7opBd9yz7/B77n\nJSxcIVucRgZERERyTsmAiIhIzikZEBERyTklAyIiIjmnZEBERCTnlAyIiIjknJIBERGRnFMyICIi\nknNKBkRERHJOyYCIiEjOKRkQERHJuVbzbIInntsTPuudKfaE3S5zzbuTzXfXZ9MwzhV/xMyb3WV8\nNXMVV/zZvc5xl/F5WNUV33HeN674/bvc64oHeDn0ccXvuOh5dxmP9dzWFd//NP996S+Y5Ys/0p70\nfWDsXb74d6qg/4W+z1TQ1i89Q+fe2Z+n8Ej341zzf8m+dMUP6ODflu86YW9X/BnXOW/qD6x59DRX\n/NP0dZex7Y/GuuLv7burK37/Xo+44gHsYN89/f+0wPecAYDbjvHFn+SsE8D1fQ5xxZ/DH1zx+7V9\nH7i9wbgmHRkws+3N7D4z+9jMFptZrae7mNmFZvaJmc03s0fMrGdT1kFEykvtXqT1a+rDBJ2AV4ET\ngFopkpmdAfwGOAbYCpgHjDaz5Zu4HiJSPmr3Iq1ckx4mCCGMAkYBmFmp5yyeDFwUQrg/iTkMmAYc\nAIxoyrqISHmo3Yu0fmU7gdDMNgDWBB4rTAshzAZeABrxgG0RaenU7kVah3JeTbAmcQgxfbbLtOQ9\nEVn2qN2LtAIt4dJCo8RxRhFZpqndi7Qg5by0cCqxA1iDmnsJqwMNX7dy7SDo1LXmtB36w44Dmq6G\nIsuKh4bBqOE1Js2f57uUrok0ut2/89sbWK5LpxrT1uy/PWsN8F8aJ5IHC4bdy9fD768x7ZEv52X6\nbNmSgRDCJDObCuwMvA5gZp2BrYGrG5zBMUOgV7b7DIjk3p4D4qtIx3eqmN1/y7JWY2na/XevOJLO\nvTds/kqKLCM6DNifDgP2rzFt16r3uWHLnRr8bJMmA2bWCehJ3BMA6GFmmwMzQwgfAlcC55jZBGAy\ncBHwEeC/m4eItAhq9yKtX1OPDGwJPEE8FhiAy5PpNwO/CiFcamYdgWuBrsDTwJ4hBN9t7USkJVG7\nF2nlmvo+A0/SwEmJIYTzgfObslwRqRy1e5HWryVcTSAiIiIV1GoeVLTcxnOxzWdnij3AeShy5zDG\nX6G1fQ+9uOJT30OHAKy62hV/8T7+3C78rNQN4+pmRyx2xR/1qr9Oozbv5/vAn311Atiqur3vA219\n3xPAYOf6e9s2cMVftflRrvhp1fP4o+sTlTU+bMxyYdPM8T0+ftM1/ys40xW/cMezXfEAM1jN94Gj\n/dvy4ff72tiT+27tLoO+vnqNa+ur037H+NsXT/jal23l74sGXuKs1+98dQI4yrn+xu63hSt+fVsx\nU5xGBkRERHJOyYCIiEjOKRkQERHJOSUDIiIiOadkQEREJOeUDIiIiOSckgEREZGcUzIgIiKSc0oG\nREREck7JgIiISM4pGRAREcm5VvNsgqO7Xkf3VVfPFPtPO9Y17yGLR7rr81/n7b2PPyi4ywj/auuK\n/+z+bPegLjaQ21zxDx/vq9O11xzmigdYhw9d8eFxX50Avtmpoyve+i50l7H4r756HXHy4674F9v8\nyBXfufOrwHDXZypp9vRV4KNsbR7g1+v+y1fAx5u5wtv91zd7gOP2udkV/9Y6Q91lbDzDF7/l1y+5\ny7igo29bHryLb/72biP6x4t9dQrvu4ugzURfvRbf5O+LLj7it674NZjmiu+UcZ9fIwMiIiI5p2RA\nREQk55QMiIiI5JySARERkZxTMiAiIpJzSgZERERyTsmAiIhIzikZEBERyTklAyIiIjmnZEBERCTn\nlAyIiIjknJIBERGRnGs1Dyoaz8ZMZf1MsffcPNA170U/8X8NYeRi3wfO9+ddd52/tyv+iNk3ussY\n2eUgV/yka9Zwxd/E/7niAV76ZCtX/NidNneXsemFE13x4Tzn+ga6zPM9QWbuY92cJUxwxn/mjK+w\nG5eDbu0yh1/668Gu2S/cIvu8AS4/+ixXPIDd52v33/+d/4E9jHOG9/2Bu4jB4193xc9xNsnO2/ji\nAbjV913N/9rcRXT8R7Ur3ob7+/nVbborfg4rueLbsWqmOI0MiIiI5JySARERkZxTMiAiIpJzSgZE\nRERyTsmAiIhIzikZEBERyTklAyIiIjmnZEBERCTnlAyIiIjknJIBERGRnFMyICIiknOt5tkEr7Mp\n7ch2T+3Q1XcP6hN7Xequz9/ea+uKHzL4OHcZv+VqV/xB0/y53fiVfd9Vj2rfvbq72r2ueIDpa63m\nil+Vee4ypp3X1RW/8izf+gbYvMtoV/wz393eWUJPZ/xsZ3yFvWHQ0bF9Puab/ZDbfc8auPJ5/3MD\nzjv+dFf84P3+7C4D5/3w7+Gn7iK23uhVV/zbX23mir8z/NwVD3CpneuK73hhI/Z9L3a2+738RRzz\n2lDfB+b6wgd+XpUpTiMDIiIiOadkQEREJOeUDIiIiOSckgEREZGcUzIgIiKSc0oGREREck7JgIiI\nSM4pGRAREck5JQMiIiI5p2RAREQk55QMiIiI5JySARERkZxrNQ8q+nzY2jCmR7bgH/rmfbf5H9wx\np9dKrvjfcrm7jNDf95CMz4d3cpfx/b6+h/y8xvdc8QO8KwMYbgNc8cdP8j9E6KoeZ7vi1+ryqbuM\nZ4bu4oo/4dDLXPFXM9AVD5854ytr/dvH06F39vgP563jmv+8Cd1c8eEDVzgAF83wPUxnfLeb3WWM\n+O9i3wcO9j9wae4iX/918eJ3XPGXXj3YFQ/Ah+f74vv4i7i7/x6u+J+/f7+7jAN73OaK34ZnXfEr\nVq3P7RnimnRkwMy2N7P7zOxjM1tsZvul3r8xmV78erAp6yAi5aV2L9L6NfVhgk7Aq8AJQF3p50PA\nGsCaycu3GygiLY3avUgr16SHCUIIo4BRAGZW14PIvw4hzGjKckWkctTuRVq/SpxAuIOZTTOzt83s\nGjNbpQJ1EJHyUrsXacHKfQLhQ8BIYBKwIfAn4EEz+0kIwX9Wi4i0Bmr3Ii1cWZOBEMKIoj/fNLM3\ngInADsAT9X74nkGwQtea03r3hz469ChS23+Ae2tMmT9/TkVq0th2P23Q5bTtWvOqnc79d6fzAN8Z\n3iJ5MXbYe4wd/l6Nact92THTZyt6aWEIYZKZfQb0pKFk4KdDYB3HdUYiuXZA8lqiY8c3mT278j+k\nWdv9GkNOpUPvjctXMZFWbosBvdhiQK8a01asWp9jtzyjwc9W9KZDZrY2sCrgv4hbRFoltXuRlqdJ\nRwbMrBMx2y+cUdzDzDYHZiavwcRjh1OTuEuAd4HRTVkPESkftXuR1q+pDxNsSRz2C8mrcNuqm4Hj\ngc2Aw4CuwCfEzuC8EMLCJq6HiJSP2r1IK9fU9xl4kvoPPVT+gKWINCm1e5HWr9U8m4CnDTrXdT+T\nlEW+q5VO7Pt3d3XOPsN3r+6tLnnRXcb3//lew0FF1hg8211Grydec8X3sSpXfI8w0RUP8Med/+iK\nv+OxX7rLGDNuV98H3nQXwSaHveyKv/q105wlXOCMn+mMr6zJ534PVt0i+wd+4Cygi/Oqxsu93zf0\nuNT33JMRAw93l8Hwj33xJ3d3F3HNc6e64m/Z/FBX/OzP13TFAzDLt/6Cv3vk+/aWK75njzfcZdx5\n1WGu+LYHLHDFD5z9aqY4PbVQREQk55QMiIiI5JySARERkZxTMiAiIpJzSgZERERyrnUnA1OHVboG\nFXHHyErXoEKm5XN981BOl7uUSTn+Lj7I57IPG1fpGlTGsFfK+wwvJQOtUH6TgeGVrkFljMrpcpcy\nOcffxQf5XPa8JgPDfVdxL7XWnQyIiIjIUlMyICIiknNKBkRERHKuNdyOuAMA88bXfmfRLJhd4sDK\nR74TLz6umu6v1TTfAZ0Pqz5zFzF2Tunps2bD2BJ3EQ6f+g8yLah61xU/0z5wxa8QprriAZhTx3Is\n+rLke3OrfLdtBmDiqr74Sf4ivqp62/eBd+vIzed+CeNLfSe+JwBXV88q/LeD64PlF+s3q8T3t/BL\n+LyO7cN5V15mNRxSQ/A/cXlBVYl+qz4zF9X93sIvYWapZXf2X9On+eIB3vGFVwffra+r6vlqZ31d\nx/vOLjVM8cUDTK762hW/IDjXN1D1Yenfq1lf1fHeuLGu+c/86NuVV2+7txDKe8ail5kNBG6rdD1E\nliEHhxBur3Ql6qI2L9Is6m33rSEZWBXYHZgM+J7QICLFOgDrA6NDCJ9XuC51UpsXaVKZ2n2LTwZE\nRESkeekEQhERkZxTMiAiIpJzSgZERERyTsmAiIhIzikZEBERyblWmQyY2QlmNsnMvjKz583sR5Wu\nU3Mzs8Fmtjj1eqvS9WpqZra9md1nZh8ny7hfiZgLzewTM5tvZo+YWc9K1LUpNbTcZnZjifX/YKXq\nWwl5a/dq8zVilrk2Dy2r3be6ZMDMfglcDgwGtgBeA0ab2WoVrVh5jAPWANZMXttVtjrNohPwKnAC\nUOu6VzM7A/gNcAywFTCPuP6XL2clm0G9y514iJrrf0B5qlZ5OW73avPLbpuHFtTuW8PtiNMGAdeG\nEIYCmNmxwN7Ar4BLK1mxMlgUQphR6Uo0pxDCKGAUgJlZiZCTgYtCCPcnMYcB04ADgBHlqmdTy7Dc\nAF8v6+u/Hnlt92rzy2ibh5bV7lvVyICZtQP6AI8VpoV416RHgZ9Uql5l1CsZTppoZrea2TqVrlA5\nmdkGxMy4eP3PBl4gH+t/BzObZmZvm9k1ZrZKpStUDjlv92rz+W7zUKZ236qSAWA1oC0xKyw2jbjB\nLMueB44g3qb1WGAD4Ckz61TJSpXZmsShtDyu/4eAw4CdgNOBfsCD9exNLEvy2u7V5vPd5qGM7b41\nHiYoxaj7eMsyIYQwuujPcWb2IvAB8AvgxsrUqsXIw/ovHg5908zeACYCOwBPVKRSlbdMr3e1+Xot\n0+u+oJztvrWNDHwGVBNPpii2OrUzx2VaCGEW8C6wTJxVm9FUYieg9R/CJGJ7yMP6V7tHbT41PVfr\nvqA5232rSgZCCAuBV4CdC9OS4ZKdgWcrVa9KMLMVgQ3xPtS+FUsawlRqrv/OwNbkb/2vDaxKDta/\n2n2kNh/ltc1D87b71niY4ArgZjN7BXiReJZxR+CmSlaquZnZZcD9xGHC7sAFwCJgWCXr1dSS46E9\niXsDAD3MbHNgZgjhQ+BK4Bwzm0B8xO1FwEfAvRWobpOpb7mT12BgJLFj7AlcQtxLHF17bsuk3LV7\ntfllu81DC2v3IYRW9wKOJ24UXwHPAVtWuk5lWOZhxAbwFTAFuB3YoNL1aobl7AcsJg4LF7/+XRRz\nPvAJMD9pFD0rXe/mXG7i88hHJR3CAuB94B9At0rXu8zfUa7avdr8st3mG1r2crd7SyokIiIiOdWq\nzhkQERGRpqdkQEREJOeUDIiIiOSckgEREZGcUzIgIiKSc0oGREREck7JgIiISM4pGRAREck5JQMi\nIiI5p2RAREQk55QMiIiI5Nz/A+/fWb9QS6P+AAAAAElFTkSuQmCC\n",
       "text/plain": [
-       "<matplotlib.figure.Figure at 0x7f4c8c1169b0>"
+       "<matplotlib.figure.Figure at 0x7f4ad8fc5940>"
       ]
      },
      "metadata": {},
diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py
index 5256dd102..1f13e18ee 100644
--- a/openmc/mgxs_library.py
+++ b/openmc/mgxs_library.py
@@ -1735,108 +1735,82 @@ class XSdata(object):
             # Get the sparse scattering data to print to the library
             G = self.energy_groups.num_groups
             if self.representation == 'isotropic':
-                g_out_bounds = np.zeros((G, 2), dtype=np.int)
-                for g_in in range(G):
-                    if self.scatter_format == 'legendre':
-                        nz = np.nonzero(self._scatter_matrix[i][g_in, :, 0])
-                    elif self.scatter_format == 'histogram':
-                        nz = np.nonzero(np.sum(
-                            self._scatter_matrix[i][g_in, :, :], axis=1))
-                    if len(nz[0]) == 0:
-                        g_out_bounds[g_in, 0] = 0
-                        g_out_bounds[g_in, 1] = 0
-                    else:
-                        g_out_bounds[g_in, 0] = nz[0][0]
-                        g_out_bounds[g_in, 1] = nz[0][-1]
-
-                # Now create the flattened scatter matrix array
-                matrix = self._scatter_matrix[i]
-                flat_scatt = []
-                for g_in in range(G):
-                    for g_out in range(g_out_bounds[g_in, 0],
-                                       g_out_bounds[g_in, 1] + 1):
-                        for l in range(len(matrix[g_in, g_out, :])):
-                            flat_scatt.append(matrix[g_in, g_out, l])
-
-                # And write it.
-                scatt_grp = xs_grp.create_group('scatter_data')
-                scatt_grp.create_dataset("scatter_matrix",
-                                         data=np.array(flat_scatt))
-
-                # Repeat for multiplicity
-                if self._multiplicity_matrix[i] is not None:
-                    # Now create the flattened scatter matrix array
-                    matrix = self._multiplicity_matrix[i][:, :]
-                    flat_mult = []
-                    for g_in in range(G):
-                        for g_out in range(g_out_bounds[g_in, 0],
-                                           g_out_bounds[g_in, 1] + 1):
-                            flat_mult.append(matrix[g_in, g_out])
-
-                    scatt_grp.create_dataset("multiplicity matrix",
-                                             data=np.array(flat_mult))
-
-                # And finally, adjust g_out_bounds for 1-based group counting
-                # and write it.
-                g_out_bounds[:, :] += 1
-                scatt_grp.create_dataset("g_min", data=g_out_bounds[:, 0])
-                scatt_grp.create_dataset("g_max", data=g_out_bounds[:, 1])
-
+                Np = 1
+                Na = 1
             elif self.representation == 'angle':
                 Np = self.num_polar
                 Na = self.num_azimuthal
-                g_out_bounds = np.zeros((Np, Na, G, 2), dtype=np.int)
 
-                for p in range(Np):
-                    for a in range(Na):
-                        for g_in in range(G):
-                            if self.scatter_format == 'legendre':
+            g_out_bounds = np.zeros((Np, Na, G, 2), dtype=np.int)
+            for p in range(Np):
+                for a in range(Na):
+                    for g_in in range(G):
+                        if self.scatter_format == 'legendre':
+                            if self.representation == 'isotropic':
+                                matrix = \
+                                    self._scatter_matrix[i][g_in, :, 0]
+                            elif self.representation == 'angle':
                                 matrix = \
                                     self._scatter_matrix[i][p, a, g_in, :, 0]
-                            elif self.scatter_format == 'histogram':
+                        elif self.scatter_format == 'histogram':
+                            if self.representation == 'isotropic':
+                                matrix = \
+                                    np.sum(self._scatter_matrix[i][g_in, :, :],
+                                           axis=1)
+                            elif self.representation == 'angle':
                                 matrix = \
                                     np.sum(self._scatter_matrix[i][p, a, g_in, :, :],
                                            axis=1)
-                            nz = np.nonzero(matrix)
-                            g_out_bounds[p, a, g_in, 0] = nz[0][0]
-                            g_out_bounds[p, a, g_in, 1] = nz[0][-1]
+                        nz = np.nonzero(matrix)
+                        g_out_bounds[p, a, g_in, 0] = nz[0][0]
+                        g_out_bounds[p, a, g_in, 1] = nz[0][-1]
+
+            # Now create the flattened scatter matrix array
+            flat_scatt = []
+            for p in range(Np):
+                for a in range(Na):
+                    if self.representation == 'isotropic':
+                        matrix = self._scatter_matrix[i][:, :, :]
+                    elif self.representation == 'angle':
+                        matrix = self._scatter_matrix[i][p, a, :, :, :]
+                    for g_in in range(G):
+                        for g_out in range(g_out_bounds[p, a, g_in, 0],
+                                           g_out_bounds[p, a, g_in, 1] + 1):
+                            for l in range(len(matrix[g_in, g_out, :])):
+                                flat_scatt.append(matrix[g_in, g_out, l])
+
+            # And write it.
+            scatt_grp = xs_grp.create_group('scatter_data')
+            scatt_grp.create_dataset("scatter_matrix",
+                                     data=np.array(flat_scatt))
+
+            # Repeat for multiplicity
+            if self._multiplicity_matrix[i] is not None:
 
                 # Now create the flattened scatter matrix array
-                flat_scatt = []
+                flat_mult = []
                 for p in range(Np):
                     for a in range(Na):
-                        matrix = self._scatter_matrix[i][p, a, :, :, :]
+                        if self.representation == 'isotropic':
+                            matrix = self._multiplicity_matrix[i][:, :]
+                        elif self.representation == 'angle':
+                            matrix = self._multiplicity_matrix[i][p, a, :, :]
                         for g_in in range(G):
                             for g_out in range(g_out_bounds[p, a, g_in, 0],
                                                g_out_bounds[p, a, g_in, 1] + 1):
-                                for l in range(len(matrix[g_in, g_out, :])):
-                                    flat_scatt.append(matrix[g_in, g_out, l])
+                                flat_mult.append(matrix[g_in, g_out])
 
                 # And write it.
-                scatt_grp = xs_grp.create_group('scatter_data')
-                scatt_grp.create_dataset("scatter_matrix",
-                                         data=np.array(flat_scatt))
+                scatt_grp.create_dataset("multiplicity_matrix",
+                                         data=np.array(flat_mult))
 
-                # Repeat for multiplicity
-                if self._multiplicity_matrix[i] is not None:
-
-                    # Now create the flattened scatter matrix array
-                    flat_mult = []
-                    for p in range(Np):
-                        for a in range(Na):
-                            matrix = self._multiplicity_matrix[i][p, a, :, :]
-                            for g_in in range(G):
-                                for g_out in range(g_out_bounds[p, a, g_in, 0],
-                                                   g_out_bounds[p, a, g_in, 1] + 1):
-                                    flat_mult.append(matrix[g_in, g_out])
-
-                    # And write it.
-                    scatt_grp.create_dataset("multiplicity_matrix",
-                                             data=np.array(flat_mult))
-
-                # And finally, adjust g_out_bounds for 1-based group counting
-                # and write it.
-                g_out_bounds[:, :, :, :] += 1
+            # And finally, adjust g_out_bounds for 1-based group counting
+            # and write it.
+            g_out_bounds[:, :, :, :] += 1
+            if self.representation == 'isotropic':
+                scatt_grp.create_dataset("g_min", data=g_out_bounds[0, 0, :, 0])
+                scatt_grp.create_dataset("g_max", data=g_out_bounds[0, 0, :, 1])
+            elif self.representation == 'angle':
                 scatt_grp.create_dataset("g_min", data=g_out_bounds[:, :, :, 0])
                 scatt_grp.create_dataset("g_max", data=g_out_bounds[:, :, :, 1])
 

From 5af55a74d49520665da4e852522249ab93e00bfe Mon Sep 17 00:00:00 2001
From: Adam Nelson <nelsonag@umich.edu>
Date: Sun, 13 Nov 2016 15:38:19 -0500
Subject: [PATCH 09/12] Removed errant pdb import

---
 openmc/mgxs/mgxs.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py
index 1337b1fa7..06527d5c2 100644
--- a/openmc/mgxs/mgxs.py
+++ b/openmc/mgxs/mgxs.py
@@ -3775,7 +3775,7 @@ class ScatterMatrixXS(MatrixMGXS):
         # tally data is stored in order of increasing energies
         if order_groups == 'increasing':
             xs = xs[:, ::-1, ::-1, ...]
-        # import pdb; pdb.set_trace()
+
         if squeeze:
             # We want to squeeze out everything but the in_groups, out_groups,
             # and, if needed, num_mu_bins dimension. These must not be squeezed

From a29e719c501afd198c42de22d4b15a862518defc Mon Sep 17 00:00:00 2001
From: Adam Nelson <nelsonag@umich.edu>
Date: Mon, 14 Nov 2016 19:07:16 -0500
Subject: [PATCH 10/12] Removed mistake from a poor merge

---
 openmc/mgxs/library.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py
index 97bb6d96b..513884d96 100644
--- a/openmc/mgxs/library.py
+++ b/openmc/mgxs/library.py
@@ -946,6 +946,7 @@ class Library(object):
         if nuclide != 'total':
             name += '_' + nuclide
         xsdata = openmc.XSdata(name, self.energy_groups)
+        xsdata.num_delayed_groups = self.num_delayed_groups
 
         # Right now only isotropic weighting is supported
         self.representation = 'isotropic'

From ecd3112e21d0cb6908e116e7e224f9cab166a133 Mon Sep 17 00:00:00 2001
From: Adam Nelson <nelsonag@umich.edu>
Date: Mon, 14 Nov 2016 20:47:38 -0500
Subject: [PATCH 11/12] Updating example problems for the new format

---
 .../python/pincell_multigroup/build-xml.py    |  10 ++++++++--
 examples/xml/pincell_multigroup/mgxs.h5       | Bin 16696 -> 16696 bytes
 2 files changed, 8 insertions(+), 2 deletions(-)

diff --git a/examples/python/pincell_multigroup/build-xml.py b/examples/python/pincell_multigroup/build-xml.py
index ee74233c1..cc57b61e2 100644
--- a/examples/python/pincell_multigroup/build-xml.py
+++ b/examples/python/pincell_multigroup/build-xml.py
@@ -1,3 +1,5 @@
+import numpy as np
+
 import openmc
 import openmc.mgxs
 
@@ -26,7 +28,7 @@ uo2_xsdata.set_total(
      0.5644058])
 uo2_xsdata.set_absorption([8.0248E-03, 3.7174E-03, 2.6769E-02, 9.6236E-02,
                            3.0020E-02, 1.1126E-01, 2.8278E-01])
-uo2_xsdata.set_scatter_matrix(
+scatter_matrix = np.array(
     [[[0.1275370, 0.0423780, 0.0000094, 0.0000000, 0.0000000, 0.0000000, 0.0000000],
       [0.0000000, 0.3244560, 0.0016314, 0.0000000, 0.0000000, 0.0000000, 0.0000000],
       [0.0000000, 0.0000000, 0.4509400, 0.0026792, 0.0000000, 0.0000000, 0.0000000],
@@ -34,6 +36,8 @@ uo2_xsdata.set_scatter_matrix(
       [0.0000000, 0.0000000, 0.0000000, 0.0001253, 0.2714010, 0.0102550, 0.0000000],
       [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0012968, 0.2658020, 0.0168090],
       [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0085458, 0.2730800]]])
+scatter_matrix = np.swapaxes(np.swapaxes(scatter_matrix, 0, 1), 1, 2)
+uo2_xsdata.set_scatter_matrix(scatter_matrix)
 uo2_xsdata.set_fission([7.21206E-03, 8.19301E-04, 6.45320E-03,
                         1.85648E-02, 1.78084E-02, 8.30348E-02,
                         2.16004E-01])
@@ -50,7 +54,7 @@ h2o_xsdata.set_total([0.15920605, 0.412969593, 0.59030986, 0.58435,
 h2o_xsdata.set_absorption([6.0105E-04, 1.5793E-05, 3.3716E-04,
                            1.9406E-03, 5.7416E-03, 1.5001E-02,
                            3.7239E-02])
-h2o_xsdata.set_scatter_matrix(
+scatter_matrix = np.array(
     [[[0.0444777, 0.1134000, 0.0007235, 0.0000037, 0.0000001, 0.0000000, 0.0000000],
       [0.0000000, 0.2823340, 0.1299400, 0.0006234, 0.0000480, 0.0000074, 0.0000010],
       [0.0000000, 0.0000000, 0.3452560, 0.2245700, 0.0169990, 0.0026443, 0.0005034],
@@ -58,6 +62,8 @@ h2o_xsdata.set_scatter_matrix(
       [0.0000000, 0.0000000, 0.0000000, 0.0000714, 0.1391380, 0.5118200, 0.0612290],
       [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0022157, 0.6999130, 0.5373200],
       [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.1324400, 2.4807000]]])
+scatter_matrix = np.swapaxes(np.swapaxes(scatter_matrix, 0, 1), 1, 2)
+h2o_xsdata.set_scatter_matrix(scatter_matrix)
 
 mg_cross_sections_file = openmc.MGXSLibrary(groups)
 mg_cross_sections_file.add_xsdatas([uo2_xsdata, h2o_xsdata])
diff --git a/examples/xml/pincell_multigroup/mgxs.h5 b/examples/xml/pincell_multigroup/mgxs.h5
index 1d5561b0006a8afb89663ca33573eedaabad7a91..ae3b17e90f2dafb1a350d9bd1392d61254f2a37c 100644
GIT binary patch
delta 194
zcmdnd#JHo0af2r(k9%yiyLxQ2e^E+m(PUoT$&(i-N^IiFVPW2wp|v@mPn8$Mm^@##
zeY2d*IYzMPds$ux<G&n~!LI;i$Sc-Dbk0}xWr8SHhbh)zVnMNOGq3J@h~oc-s;p3J
fIyY}HKLatn-tsIn*a&}Hs9E#vA!bc}uigm&6H`Mw

delta 194
zcmdnd#JHo0af2r(kAG20YEf*oI}oW)=GC1%d4ZzDCaxS7=GeK?oAdcpc|nZH^Htk7
z%gLN$1dG0x<%KZ*%Rw3Z3Q&f;Vl71Hd{tj2h+=h^Vhtu16x%lQ>b{34{%@$t3bm$l
d^9J)X5Yy`|&oYCJ@VAAUHQyd$*5vo<od7nSK#2eV


From 6ba1f046311aef9ef9b67985b494fbd00fe0caff Mon Sep 17 00:00:00 2001
From: Adam Nelson <nelsonag@umich.edu>
Date: Sun, 27 Nov 2016 06:34:47 -0500
Subject: [PATCH 12/12] Resolving @paulromano comments

---
 .../python/pincell_multigroup/build-xml.py    |  4 +--
 openmc/filter.py                              | 11 +++++---
 openmc/mgxs/library.py                        |  2 +-
 openmc/mgxs/mgxs.py                           | 25 +++++++++++++------
 openmc/mgxs_library.py                        | 10 ++++++--
 tests/test_tallies/test_tallies.py            |  3 ++-
 6 files changed, 38 insertions(+), 17 deletions(-)

diff --git a/examples/python/pincell_multigroup/build-xml.py b/examples/python/pincell_multigroup/build-xml.py
index cc57b61e2..9cc23300d 100644
--- a/examples/python/pincell_multigroup/build-xml.py
+++ b/examples/python/pincell_multigroup/build-xml.py
@@ -36,7 +36,7 @@ scatter_matrix = np.array(
       [0.0000000, 0.0000000, 0.0000000, 0.0001253, 0.2714010, 0.0102550, 0.0000000],
       [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0012968, 0.2658020, 0.0168090],
       [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0085458, 0.2730800]]])
-scatter_matrix = np.swapaxes(np.swapaxes(scatter_matrix, 0, 1), 1, 2)
+scatter_matrix = np.rollaxis(scatter_matrix, 0, 3)
 uo2_xsdata.set_scatter_matrix(scatter_matrix)
 uo2_xsdata.set_fission([7.21206E-03, 8.19301E-04, 6.45320E-03,
                         1.85648E-02, 1.78084E-02, 8.30348E-02,
@@ -62,7 +62,7 @@ scatter_matrix = np.array(
       [0.0000000, 0.0000000, 0.0000000, 0.0000714, 0.1391380, 0.5118200, 0.0612290],
       [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0022157, 0.6999130, 0.5373200],
       [0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.1324400, 2.4807000]]])
-scatter_matrix = np.swapaxes(np.swapaxes(scatter_matrix, 0, 1), 1, 2)
+scatter_matrix = np.rollaxis(scatter_matrix, 0, 3)
 h2o_xsdata.set_scatter_matrix(scatter_matrix)
 
 mg_cross_sections_file = openmc.MGXSLibrary(groups)
diff --git a/openmc/filter.py b/openmc/filter.py
index 61818f89f..a7a8b934a 100644
--- a/openmc/filter.py
+++ b/openmc/filter.py
@@ -100,10 +100,13 @@ class Filter(object):
     @classmethod
     def _recursive_subclasses(cls):
         """Return all subclasses and their subclasses, etc."""
-        subs = cls.__subclasses__()
-        subsubs = [grand for s in subs for grand in s.__subclasses__()]
-        subsubsubs = [grand for s in subsubs for grand in s.__subclasses__()]
-        return subs + subsubs + subsubsubs
+        all_subclasses = []
+
+        for subclass in cls.__subclasses__():
+            all_subclasses.append(subclass)
+            all_subclasses.extend(subclass._recursive_subclasses())
+
+        return all_subclasses
 
     @classmethod
     def from_hdf5(cls, group, **kwargs):
diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py
index 2e2b7f908..7e0289064 100644
--- a/openmc/mgxs/library.py
+++ b/openmc/mgxs/library.py
@@ -59,7 +59,7 @@ class Library(object):
         The spatial domain(s) for which MGXS in the Library are computed
     correction : {'P0', None}
         Apply the P0 correction to scattering matrices if set to 'P0'
-    scatter_format : {'legendre', or 'histogram'}
+    scatter_format : {'legendre', 'histogram'}
         Representation of the angular scattering distribution (default is
         'legendre')
     legendre_order : int
diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py
index db31e2f32..58b8a32c2 100644
--- a/openmc/mgxs/mgxs.py
+++ b/openmc/mgxs/mgxs.py
@@ -1529,9 +1529,19 @@ class MGXS(object):
         else:
             df = df.drop('score', axis=1)
 
+        # Determine if change-in-angle bins are included in the MGXS to
+        # properly tile the group boundaries
+        if 'mu low' in df:
+            # Find the length of the mu filters indirectly from the number
+            # of times the mu bins repeats.
+            num_mu = int(df.shape[0] /
+                         df[df['mu low'] == df['mu low'][0]].shape[0])
+        else:
+            num_mu = 1
+
         # Override energy groups bounds with indices
         all_groups = np.arange(self.num_groups, 0, -1, dtype=np.int)
-        all_groups = np.repeat(all_groups, len(query_nuclides))
+        all_groups = np.repeat(all_groups, len(query_nuclides) * num_mu)
         if 'energy low [eV]' in df and 'energyout low [eV]' in df:
             df.rename(columns={'energy low [eV]': 'group in'},
                       inplace=True)
@@ -3791,10 +3801,8 @@ class ScatterMatrixXS(MatrixMGXS):
             # than 1, so try each axis in axes one at a time, catching the
             # ValueError as needed.
             for axis in axes:
-                try:
+                if xs.shape[axis] == 1:
                     xs = np.squeeze(xs, axis=axis)
-                except ValueError:
-                    pass
 
         return xs
 
@@ -3874,9 +3882,12 @@ class ScatterMatrixXS(MatrixMGXS):
                 df = df[df['moment'] == 'P{}'.format(moment)]
 
         elif self.scatter_format == 'histogram':
-            # Add a change-in-angle (mu) column to dataframe
-            ###TODO NOT SURE I NEED TO DO THIS
-            pass
+            # Replace the mu low and mu high columns with a single mu bin
+            del df['mu high']
+            df.rename(columns={'mu low': 'mu bins'}, inplace=True)
+            bins = [i + 1 for i in range(self.histogram_bins)]
+            bins = np.tile(bins, int(df.shape[0] / len(bins)))
+            df['mu bins'] = bins
 
         return df
 
diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py
index 25452bdb4..347c015e7 100644
--- a/openmc/mgxs_library.py
+++ b/openmc/mgxs_library.py
@@ -1773,8 +1773,14 @@ class XSdata(object):
                                     np.sum(self._scatter_matrix[i][p, a, g_in, :, :],
                                            axis=1)
                         nz = np.nonzero(matrix)
-                        g_out_bounds[p, a, g_in, 0] = nz[0][0]
-                        g_out_bounds[p, a, g_in, 1] = nz[0][-1]
+                        # It is possible that there only zeros in matrix
+                        # and therefore nz will be empty, in that case set
+                        # g_out_bounds to 0s
+                        if len(nz[0]) == 0:
+                            g_out_bounds[p, a, g_in, :] = 0
+                        else:
+                            g_out_bounds[p, a, g_in, 0] = nz[0][0]
+                            g_out_bounds[p, a, g_in, 1] = nz[0][-1]
 
             # Now create the flattened scatter matrix array
             flat_scatt = []
diff --git a/tests/test_tallies/test_tallies.py b/tests/test_tallies/test_tallies.py
index 9e525d2cf..033c7b812 100644
--- a/tests/test_tallies/test_tallies.py
+++ b/tests/test_tallies/test_tallies.py
@@ -3,13 +3,14 @@
 import os
 import sys
 sys.path.insert(0, os.pardir)
-import numpy as np
+
 from testing_harness import PyAPITestHarness
 from openmc.filter import *
 from openmc import Mesh, Tally, Tallies
 from openmc.source import Source
 from openmc.stats import Box
 
+
 class TalliesTestHarness(PyAPITestHarness):
     def _build_inputs(self):
         # Build default materials/geometry