diff --git a/include/openmc/endf.h b/include/openmc/endf.h
index 7ec078e93f..27958db5c5 100644
--- a/include/openmc/endf.h
+++ b/include/openmc/endf.h
@@ -57,7 +57,7 @@ public:
   //! Evaluate the polynomials
   //! \param[in] x independent variable
   //! \return Polynomial evaluated at x
-  double operator()(double x) const;
+  double operator()(double x) const override;
 private:
   std::vector<double> coef_; //!< Polynomial coefficients
 };
@@ -77,7 +77,7 @@ public:
   //! Evaluate the tabulated function
   //! \param[in] x independent variable
   //! \return Function evaluated at x
-  double operator()(double x) const;
+  double operator()(double x) const override;
 
   // Accessors
   const std::vector<double>& x() const { return x_; }
@@ -96,13 +96,32 @@ private:
 //==============================================================================
 
 class CoherentElasticXS : public Function1D {
+public:
   explicit CoherentElasticXS(hid_t dset);
-  double operator()(double E) const;
+
+  double operator()(double E) const override;
+
+  const std::vector<double>& bragg_edges() const { return bragg_edges_; }
+  const std::vector<double>& factors() const { return factors_; }
 private:
   std::vector<double> bragg_edges_; //!< Bragg edges in [eV]
   std::vector<double> factors_;     //!< Partial sums of structure factors [eV-b]
 };
 
+//==============================================================================
+//! Incoherent elastic scattering cross section
+//==============================================================================
+
+class IncoherentElasticXS : public Function1D {
+public:
+  explicit IncoherentElasticXS(hid_t dset);
+
+  double operator()(double E) const override;
+private:
+  double bound_xs_; //!< Characteristic bound xs in [b]
+  double debye_waller_; //!< Debye-Waller integral divided by atomic mass in [eV^-1]
+};
+
 //! Read 1D function from HDF5 dataset
 //! \param[in] group HDF5 group containing dataset
 //! \param[in] name Name of dataset
diff --git a/include/openmc/secondary_correlated.h b/include/openmc/secondary_correlated.h
index 8ba18d3d15..603d8c414b 100644
--- a/include/openmc/secondary_correlated.h
+++ b/include/openmc/secondary_correlated.h
@@ -38,7 +38,7 @@ public:
   //! \param[in] E_in Incoming energy in [eV]
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[out] mu Outgoing cosine with respect to current direction
-  void sample(double E_in, double& E_out, double& mu) const;
+  void sample(double E_in, double& E_out, double& mu) const override;
 
   // energy property
   std::vector<double>& energy() { return energy_; }
diff --git a/include/openmc/secondary_kalbach.h b/include/openmc/secondary_kalbach.h
index a898db67c8..0217c7e776 100644
--- a/include/openmc/secondary_kalbach.h
+++ b/include/openmc/secondary_kalbach.h
@@ -29,7 +29,7 @@ public:
   //! \param[in] E_in Incoming energy in [eV]
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[out] mu Outgoing cosine with respect to current direction
-  void sample(double E_in, double& E_out, double& mu) const;
+  void sample(double E_in, double& E_out, double& mu) const override;
 private:
   //! Outgoing energy/angle at a single incoming energy
   struct KMTable {
diff --git a/include/openmc/secondary_nbody.h b/include/openmc/secondary_nbody.h
index 8509b07143..eb90bd8a04 100644
--- a/include/openmc/secondary_nbody.h
+++ b/include/openmc/secondary_nbody.h
@@ -24,7 +24,7 @@ public:
   //! \param[in] E_in Incoming energy in [eV]
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[out] mu Outgoing cosine with respect to current direction
-  void sample(double E_in, double& E_out, double& mu) const;
+  void sample(double E_in, double& E_out, double& mu) const override;
 private:
   int n_bodies_; //!< Number of particles distributed
   double mass_ratio_; //!< Total mass of particles [neutron mass]
diff --git a/include/openmc/secondary_thermal.h b/include/openmc/secondary_thermal.h
new file mode 100644
index 0000000000..507176c2e3
--- /dev/null
+++ b/include/openmc/secondary_thermal.h
@@ -0,0 +1,113 @@
+//! \file secondary_thermal.h
+//! Angle-energy distributions for thermal scattering
+
+#ifndef OPENMC_SECONDARY_THERMAL_H
+#define OPENMC_SECONDARY_THERMAL_H
+
+#include "openmc/angle_energy.h"
+#include "openmc/endf.h"
+#include "openmc/secondary_correlated.h"
+
+#include <hdf5.h>
+#include "xtensor/xtensor.hpp"
+
+#include <vector>
+
+namespace openmc {
+
+//==============================================================================
+//! Coherent elastic scattering angle-energy distribution
+//==============================================================================
+
+class CoherentElasticAE : public AngleEnergy {
+public:
+  //! Construct from a coherent elastic scattering cross section
+  //
+  //! \param[in] xs Coherent elastic scattering cross section
+  explicit CoherentElasticAE(const CoherentElasticXS& xs);
+
+
+  //! Sample distribution for an angle and energy
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[out] mu Outgoing cosine with respect to current direction
+  void sample(double E_in, double& E_out, double& mu) const override;
+private:
+  const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section
+};
+
+//==============================================================================
+//! Incoherent elastic scattering angle-energy distribution
+//==============================================================================
+
+class IncoherentElasticAE : public AngleEnergy {
+public:
+  explicit IncoherentElasticAE(hid_t group);
+
+  //! Sample distribution for an angle and energy
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[out] mu Outgoing cosine with respect to current direction
+  void sample(double E_in, double& E_out, double& mu) const override;
+private:
+  double debye_waller_;
+};
+
+//==============================================================================
+//! Incoherent elastic scattering angle-energy distribution (discrete)
+//==============================================================================
+
+class IncoherentElasticAEDiscrete : public AngleEnergy {
+public:
+  explicit IncoherentElasticAEDiscrete(hid_t group, const std::vector<double>& energy);
+
+  void sample(double E_in, double& E_out, double& mu) const override;
+private:
+  const std::vector<double>& energy_; //!< Incoherent inelastic scattering cross section
+  xt::xtensor<double, 2> mu_out_; //!< Cosines for each incident energy
+};
+
+//==============================================================================
+//! Incoherent inelastic scattering angle-energy distribution (discrete)
+//==============================================================================
+
+class IncoherentInelasticAEDiscrete : public AngleEnergy {
+public:
+  explicit IncoherentInelasticAEDiscrete(hid_t group, const std::vector<double>& energy);
+
+  void sample(double E_in, double& E_out, double& mu) const override;
+private:
+  const std::vector<double>& energy_; //!< Incident energies
+  xt::xtensor<double, 2> energy_out_; //!< Outgoing energies for each incident energy
+  xt::xtensor<double, 3> mu_out_; //!< Outgoing cosines for each incident/outgoing energy
+  bool skewed_; //!< Whether outgoing energy distribution is skewed
+};
+
+//==============================================================================
+//! Incoherent inelastic scattering angle-energy distribution (discrete)
+//==============================================================================
+
+class IncoherentInelasticAEContinuous : public AngleEnergy {
+public:
+  explicit IncoherentInelasticAEContinuous(hid_t group);
+
+  void sample(double E_in, double& E_out, double& mu) const override;
+private:
+  //! Secondary energy/angle distribution
+  struct DistEnergySab {
+    std::size_t n_e_out; //!< Number of outgoing energies
+    xt::xtensor<double, 1> e_out;     //!< Outgoing energies
+    xt::xtensor<double, 1> e_out_pdf; //!< Probability density function
+    xt::xtensor<double, 1> e_out_cdf; //!< Cumulative distribution function
+    xt::xtensor<double, 2> mu; //!< Equiprobable angles at each outgoing energy
+  };
+
+  std::vector<double> energy_; //!< Incident energies
+  std::vector<DistEnergySab> distribution_; //!< Secondary angle-energy at
+                                            //!< each incident energy
+};
+
+
+} // namespace openmc
+
+#endif // OPENMC_SECONDARY_THERMAL_H
diff --git a/include/openmc/secondary_uncorrelated.h b/include/openmc/secondary_uncorrelated.h
index e895a17bdc..e430c75c45 100644
--- a/include/openmc/secondary_uncorrelated.h
+++ b/include/openmc/secondary_uncorrelated.h
@@ -29,7 +29,7 @@ public:
   //! \param[in] E_in Incoming energy in [eV]
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[out] mu Outgoing cosine with respect to current direction
-  void sample(double E_in, double& E_out, double& mu) const;
+  void sample(double E_in, double& E_out, double& mu) const override;
 
   // Accessors
   AngleDistribution& angle() { return angle_; }
diff --git a/openmc/data/angle_energy.py b/openmc/data/angle_energy.py
index 19dddfe19e..5fb176f681 100644
--- a/openmc/data/angle_energy.py
+++ b/openmc/data/angle_energy.py
@@ -43,6 +43,8 @@ class AngleEnergy(EqualityMixin, metaclass=ABCMeta):
             return openmc.data.IncoherentElasticAEDiscrete.from_hdf5(group)
         elif dist_type == 'incoherent_inelastic_discrete':
             return openmc.data.IncoherentInelasticAEDiscrete.from_hdf5(group)
+        elif dist_type == 'incoherent_inelastic_continuous':
+            return openmc.data.IncoherentInelasticAEContinuous.from_hdf5(group)
 
     @staticmethod
     def from_ace(ace, location_dist, location_start, rx=None):
diff --git a/openmc/data/correlated.py b/openmc/data/correlated.py
index b760f2e2fd..567d23abd4 100644
--- a/openmc/data/correlated.py
+++ b/openmc/data/correlated.py
@@ -44,6 +44,8 @@ class CorrelatedAngleEnergy(AngleEnergy):
 
     """
 
+    _name = 'correlated'
+
     def __init__(self, breakpoints, interpolation, energy, energy_out, mu):
         super().__init__()
         self.breakpoints = breakpoints
@@ -111,7 +113,7 @@ class CorrelatedAngleEnergy(AngleEnergy):
             HDF5 group to write to
 
         """
-        group.attrs['type'] = np.string_('correlated')
+        group.attrs['type'] = np.string_(self._name)
 
         dset = group.create_dataset('energy', data=self.energy)
         dset.attrs['interpolation'] = np.vstack((self.breakpoints,
diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py
index 0366af1a66..eaf8b49bce 100644
--- a/openmc/data/thermal.py
+++ b/openmc/data/thermal.py
@@ -248,14 +248,14 @@ class IncoherentElastic(Function1D):
     Parameters
     ----------
     xs : float
-        Characteristic cross section in [b]
+        Characteristic bound cross section in [b]
     debye_waller : float
         Debye-Waller integral in [eV\ :math:`^{-1}`]
 
     Attributes
     ----------
     xs : float
-        Characteristic cross section in [b]
+        Characteristic bound cross section in [b]
     debye_waller : float
         Debye-Waller integral in [eV\ :math:`^{-1}`]
 
@@ -282,7 +282,7 @@ class IncoherentElastic(Function1D):
         dataset = group.create_dataset(name)
         dataset.attrs['type'] = np.string_('incoherent')
         dataset.attrs['debye_waller'] = self.debye_waller
-        dataset.attrs['characteristic_xs'] = self.xs
+        dataset.attrs['bound_xs'] = self.xs
 
     @classmethod
     def from_hdf5(cls, dataset):
@@ -826,10 +826,10 @@ class ThermalScattering(EqualityMixin):
             elif lhtr == 2:
                 # incoherent elastic
                 params, W = endf.get_tab1_record(file_obj)
-                characteristic_xs = params[0]
+                bound_xs = params[0]
                 for T, debye_waller in zip(W.x, W.y):
                     strT = _temperature_str(T)
-                    xs[strT] = IncoherentElastic(characteristic_xs, debye_waller)
+                    xs[strT] = IncoherentElastic(bound_xs, debye_waller)
                     distribution[strT] = IncoherentElasticAE(debye_waller)
 
             elastic = ThermalScatteringReaction(xs, distribution)
diff --git a/openmc/data/thermal_angle_energy.py b/openmc/data/thermal_angle_energy.py
index e5b77b597a..d810044768 100644
--- a/openmc/data/thermal_angle_energy.py
+++ b/openmc/data/thermal_angle_energy.py
@@ -136,5 +136,6 @@ class IncoherentInelasticAEDiscrete(AngleEnergy):
         skewed = bool(group['skewed'])
         return cls(energy_out, mu_out, skewed)
 
+
 class IncoherentInelasticAEContinuous(CorrelatedAngleEnergy):
-    pass
+    _name = 'incoherent_inelastic_continuous'
diff --git a/src/endf.cpp b/src/endf.cpp
index 464bc974c2..b10bbd5502 100644
--- a/src/endf.cpp
+++ b/src/endf.cpp
@@ -234,4 +234,21 @@ double CoherentElasticXS::operator()(double E) const
   }
 }
 
+//==============================================================================
+// IncoherentElasticXS implementation
+//==============================================================================
+
+IncoherentElasticXS::IncoherentElasticXS(hid_t dset)
+{
+  read_attribute(dset, "bound_xs", bound_xs_);
+  read_attribute(dset, "debye_waller", debye_waller_);
+}
+
+double IncoherentElasticXS::operator()(double E) const
+{
+  // Determine cross section using ENDF-102, Eq. (7.5)
+  double W = debye_waller_;
+  return bound_xs_ / 2.0 * ((1 - std::exp(-4.0*E*W))/(2.0*E*W));
+}
+
 } // namespace openmc
diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp
new file mode 100644
index 0000000000..cc2fce961c
--- /dev/null
+++ b/src/secondary_thermal.cpp
@@ -0,0 +1,327 @@
+#include "openmc/secondary_thermal.h"
+
+#include "openmc/hdf5_interface.h"
+#include "openmc/random_lcg.h"
+#include "openmc/search.h"
+
+#include "xtensor/xview.hpp"
+
+#include <cmath> // for log, exp
+
+namespace openmc {
+
+// Helper function to get index on incident energy grid
+void
+get_energy_index(const std::vector<double>& energies, double E, int& i, double& f)
+{
+  // Get index and interpolation factor for elastic grid
+  i = 0;
+  f = 0.0;
+  if (E >= energies.front()) {
+    i = lower_bound_index(energies.begin(), energies.end(), E);
+    f = (E - energies[i]) / (energies[i+1] - energies[i]);
+  }
+}
+
+//==============================================================================
+// CoherentElasticAE implementation
+//==============================================================================
+
+CoherentElasticAE::CoherentElasticAE(const CoherentElasticXS& xs)
+  : xs_{xs}
+{ }
+
+void
+CoherentElasticAE::sample(double E_in, double& E_out, double& mu) const
+{
+  // Get index and interpolation factor for elastic grid
+  int i;
+  double f;
+  const auto& energies {xs_.bragg_edges()};
+  get_energy_index(energies, E_in, i, f);
+
+  // Sample a Bragg edge between 1 and i
+  const auto& factors = xs_.factors();
+  double prob = prn() * factors[i+1];
+  int k = 0;
+  if (prob >= factors.front()) {
+    k = lower_bound_index(factors.begin(), factors.begin() + (i+1), prob);
+  }
+
+  // Characteristic scattering cosine for this Bragg edge (ENDF-102, Eq. 7-2)
+  mu = 1.0 - 2.0*energies[k] / E_in;
+
+  // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7-1)
+  E_out = E_in;
+}
+
+//==============================================================================
+// IncoherentElasticAE implementation
+//==============================================================================
+
+IncoherentElasticAE::IncoherentElasticAE(hid_t group)
+{
+  read_attribute(group, "debye_waller", debye_waller_);
+}
+
+void
+IncoherentElasticAE::sample(double E_in, double& E_out, double& mu) const
+{
+  // Sample angle by inverting the distribution in ENDF-102, Eq. 7.4
+  double c = 2 * E_in * debye_waller_;
+  mu = std::log(1.0 + prn()*(std::exp(2.0*c) - 1))/c - 1.0;
+
+  // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7.4)
+  E_out = E_in;
+}
+
+//==============================================================================
+// IncoherentElasticAEDiscrete implementation
+//==============================================================================
+
+IncoherentElasticAEDiscrete::IncoherentElasticAEDiscrete(hid_t group,
+  const std::vector<double>& energy)
+  : energy_{energy}
+{
+  read_dataset(group, "mu_out", mu_out_);
+}
+
+void
+IncoherentElasticAEDiscrete::sample(double E_in, double& E_out, double& mu) const
+{
+  // Get index and interpolation factor for elastic grid
+  int i;
+  double f;
+  get_energy_index(energy_, E_in, i, f);
+
+  // Interpolate between two discrete cosines corresponding to neighboring
+  // incoming energies.
+
+  // Sample outgoing cosine bin
+  int k = prn() * mu_out_.shape()[1];
+
+  // Determine outgoing cosine corresponding to E_in[i] and E_in[i+1]
+  double mu_ik  = mu_out_(i, k);
+  double mu_i1k = mu_out_(i+1, k);
+
+  // Cosine of angle between incoming and outgoing neutron
+  mu = (1 - f)*mu_ik + f*mu_i1k;
+
+  // Energy doesn't change in elastic scattering
+  E_out = E_in;
+}
+
+//==============================================================================
+// IncoherentInelasticAEDiscrete implementation
+//==============================================================================
+
+IncoherentInelasticAEDiscrete::IncoherentInelasticAEDiscrete(hid_t group,
+  const std::vector<double>& energy)
+  : energy_{energy}
+{
+  read_dataset(group, "energy_out", energy_out_);
+  read_dataset(group, "mu_out", mu_out_);
+  read_dataset(group, "skewed", skewed_);
+}
+
+void
+IncoherentInelasticAEDiscrete::sample(double E_in, double& E_out, double& mu) const
+{
+  // Get index and interpolation factor for inelastic grid
+  int i;
+  double f;
+  get_energy_index(energy_, E_in, i, f);
+
+  // Now that we have an incoming energy bin, we need to determine the outgoing
+  // energy bin. This will depend on whether the outgoing energy distribution is
+  // skewed. If it is skewed, then the first two and last two bins have lower
+  // probabilities than the other bins (0.1 for the first and last bins and 0.4
+  // for the second and second to last bins, relative to a normal bin
+  // probability of 1). Otherwise, each bin is equally probable.
+
+  int j;
+  int n = energy_out_.shape()[1];
+  if (!skewed_) {
+    // All bins equally likely
+    j = prn() * n;
+  } else {
+    // Distribution skewed away from edge points
+    double r = prn() * (n - 3);
+    if (r > 1.0) {
+      // equally likely N-4 middle bins
+      j = r + 1;
+    } else if (r > 0.6) {
+      // second to last bin has relative probability of 0.4
+      j = n - 2;
+    } else if (r > 0.5) {
+      // last bin has relative probability of 0.1
+      j = n - 1;
+    } else if (r > 0.1) {
+      // second bin has relative probability of 0.4
+      j = 1;
+    } else {
+      // first bin has relative probability of 0.1
+      j = 0;
+    }
+  }
+
+  // Determine outgoing energy corresponding to E_in[i] and E_in[i+1]
+  double E_ij  = energy_out_(i, j);
+  double E_i1j = energy_out_(i+1, j);
+
+  // Outgoing energy
+  E_out = (1 - f)*E_ij + f*E_i1j;
+
+  // Sample outgoing cosine bin
+  int m = mu_out_.shape()[2];
+  int k = prn() * m;
+
+  // Determine outgoing cosine corresponding to E_in[i] and E_in[i+1]
+  double mu_ijk  = mu_out_(i, j, k);
+  double mu_i1jk = mu_out_(i+1, j, k);
+
+  // Cosine of angle between incoming and outgoing neutron
+  mu = (1 - f)*mu_ijk + f*mu_i1jk;
+}
+
+//==============================================================================
+// IncoherentInelasticAEContinuous implementation
+//==============================================================================
+
+IncoherentInelasticAEContinuous::IncoherentInelasticAEContinuous(hid_t group)
+{
+  // Read correlated angle-energy distribution
+  CorrelatedAngleEnergy dist {group};
+
+  // Copy incident energies
+  energy_ = dist.energy();
+
+  // Convert to S(a,b) native format
+  for (const auto& edist : dist.distribution()) {
+    // Create temporary distribution
+    DistEnergySab d;
+
+    // Copy outgoing energy distribution
+    d.n_e_out = edist.e_out.size();
+    d.e_out = edist.e_out;
+    d.e_out_pdf = edist.p;
+    d.e_out_cdf = edist.c;
+
+    for (int j = 0; j < d.n_e_out; ++j) {
+      auto adist = dynamic_cast<Tabular*>(edist.angle[j].get());
+      if (adist) {
+        // On first pass, allocate space for angles
+        if (j == 0) {
+          auto n_mu = adist->x().size();
+          d.mu = xt::empty<double>({d.n_e_out, n_mu});
+        }
+
+        // Copy outgoing angles
+        auto mu_j = xt::view(d.mu, j);
+        std::copy(adist->x().begin(), adist->x().end(), mu_j.begin());
+      }
+    }
+
+    distribution_.emplace_back(std::move(d));
+  }
+
+}
+
+void
+IncoherentInelasticAEContinuous::sample(double E_in, double& E_out, double& mu) const
+{
+  // Get index and interpolation factor for inelastic grid
+  int i;
+  double f;
+  get_energy_index(energy_, E_in, i, f);
+
+  // Sample between ith and [i+1]th bin
+  int l = f > prn() ? i + 1 : i;
+
+  // Determine endpoints on grid i
+  auto n = distribution_[i].e_out.size();
+  double E_i_1 = distribution_[i].e_out(0);
+  double E_i_J = distribution_[i].e_out(n - 1);
+
+  // Determine endpoints on grid i + 1
+  n = distribution_[i + 1].e_out.size();
+  double E_i1_1 = distribution_[i + 1].e_out(0);
+  double E_i1_J = distribution_[i + 1].e_out(n - 1);
+
+  double E_1 = E_i_1 + f * (E_i1_1 - E_i_1);
+  double E_J = E_i_J + f * (E_i1_J - E_i_J);
+
+  // Determine outgoing energy bin
+  // (First reset n_energy_out to the right value)
+  n = distribution_[l].n_e_out;
+  double r1 = prn();
+  double c_j = distribution_[l].e_out_cdf[0];
+  double c_j1;
+  std::size_t j;
+  for (j = 0; j < n - 1; ++j) {
+    c_j1 = distribution_[l].e_out_cdf[j + 1];
+    if (r1 < c_j1) break;
+    c_j = c_j1;
+  }
+
+  // check to make sure j is <= n_energy_out - 2
+  j = std::min(j, n - 2);
+
+  // Get the data to interpolate between
+  double E_l_j = distribution_[l].e_out[j];
+  double p_l_j = distribution_[l].e_out_pdf[j];
+
+  // Next part assumes linear-linear interpolation in standard
+  double E_l_j1 = distribution_[l].e_out[j + 1];
+  double p_l_j1 = distribution_[l].e_out_pdf[j + 1];
+
+  // Find secondary energy (variable E)
+  double frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j);
+  if (frac == 0.0) {
+    E_out = E_l_j + (r1 - c_j) / p_l_j;
+  } else {
+    E_out = E_l_j + (std::sqrt(std::max(0.0, p_l_j*p_l_j +
+          2.0*frac*(r1 - c_j))) - p_l_j) / frac;
+  }
+
+  // Now interpolate between incident energy bins i and i + 1
+  if (l == i) {
+    E_out = E_1 + (E_out - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1);
+  } else {
+    E_out = E_1 + (E_out - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1);
+  }
+
+  // Sample outgoing cosine bin
+  int n_mu = distribution_[l].mu.shape()[1];
+  std::size_t k = prn() * n_mu;
+
+  // Rather than use the sampled discrete mu directly, it is smeared over
+  // a bin of width min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the
+  // discrete mu value itself.
+  const auto& mu_l = distribution_[l].mu;
+  f = (r1 - c_j)/(c_j1 - c_j);
+
+  // Determine (k-1)th mu value
+  double mu_left;
+  if (k == 0) {
+    mu_left = -1.0;
+  } else {
+    mu_left = mu_l(j, k-1) + f*(mu_l(j+1, k-1) - mu_l(j, k-1));
+  }
+
+  // Determine kth mu value
+  mu = mu_l(j, k) + f*(mu_l(j+1, k) - mu_l(j, k));
+
+  // Determine (k+1)th mu value
+  double mu_right;
+  if (k == n_mu - 1) {
+    mu_right = 1.0;
+  } else {
+    mu_right = mu_l(j, k+1) + f*(mu_l(j+1, k+1) - mu_l(j, k+1));
+  }
+
+  // Smear angle
+  mu += std::min(mu - mu_left, mu_right - mu)*(prn() - 0.5);
+}
+
+} // namespace openmc