diff --git a/openmc/filter.py b/openmc/filter.py
index d11af96f52..91ceb2931a 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 c03aa78771..96bfa6fab2 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 65802f2ef0..cafda1ecf4 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 35dff1223a..5504b6e7cb 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 f9daa7a228..6922e77ac7 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]))