diff --git a/include/openmc/angle_energy.h b/include/openmc/angle_energy.h
index ac931b1b533..35e5b68b64b 100644
--- a/include/openmc/angle_energy.h
+++ b/include/openmc/angle_energy.h
@@ -16,6 +16,8 @@ class AngleEnergy {
 public:
   virtual void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const = 0;
+  virtual double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const = 0;
   virtual ~AngleEnergy() = default;
 };
 
diff --git a/include/openmc/chain.h b/include/openmc/chain.h
index a3bc6f3a364..c2a45f75adc 100644
--- a/include/openmc/chain.h
+++ b/include/openmc/chain.h
@@ -70,6 +70,8 @@ class DecayPhotonAngleEnergy : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   const Distribution* photon_energy_;
diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h
index 854cf7d7719..9f082edcdf0 100644
--- a/include/openmc/distribution.h
+++ b/include/openmc/distribution.h
@@ -9,6 +9,7 @@
 #include "pugixml.hpp"
 
 #include "openmc/constants.h"
+#include "openmc/error.h"
 #include "openmc/memory.h" // for unique_ptr
 #include "openmc/span.h"
 #include "openmc/vector.h" // for vector
@@ -23,6 +24,10 @@ class Distribution {
 public:
   virtual ~Distribution() = default;
   virtual double sample(uint64_t* seed) const = 0;
+  virtual double evaluate(double x) const
+  {
+    fatal_error("evaluate not available for this Distribution type");
+  }
 
   //! Return integral of distribution
   //! \return Integral of distribution
@@ -111,6 +116,7 @@ class Uniform : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double evaluate(double x) const override;
 
   double a() const { return a_; }
   double b() const { return b_; }
@@ -135,6 +141,7 @@ class PowerLaw : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double evaluate(double x) const override;
 
   double a() const { return std::pow(offset_, ninv_); }
   double b() const { return std::pow(offset_ + span_, ninv_); }
@@ -204,6 +211,7 @@ class Normal : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double evaluate(double x) const override;
 
   double mean_value() const { return mean_value_; }
   double std_dev() const { return std_dev_; }
@@ -227,6 +235,7 @@ class Tabular : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double evaluate(double x) const override;
 
   // properties
   vector<double>& x() { return x_; }
diff --git a/include/openmc/distribution_angle.h b/include/openmc/distribution_angle.h
index efd4e58425b..a7e021085f0 100644
--- a/include/openmc/distribution_angle.h
+++ b/include/openmc/distribution_angle.h
@@ -25,6 +25,7 @@ class AngleDistribution {
   //! \param[inout] seed pseudorandom number seed pointer
   //! \return Cosine of the angle in the range [-1,1]
   double sample(double E, uint64_t* seed) const;
+  double evaluate(double E, double mu) const;
 
   //! Determine whether angle distribution is empty
   //! \return Whether distribution is empty
diff --git a/include/openmc/distribution_multi.h b/include/openmc/distribution_multi.h
index 75126593f7d..965c832a330 100644
--- a/include/openmc/distribution_multi.h
+++ b/include/openmc/distribution_multi.h
@@ -28,6 +28,10 @@ class UnitSphereDistribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Direction sampled
   virtual Direction sample(uint64_t* seed) const = 0;
+  virtual double evaluate(Direction u) const
+  {
+    fatal_error("evaluate not available for this UnitSphereDistribution type");
+  }
 
   Direction u_ref_ {0.0, 0.0, 1.0}; //!< reference direction
 };
@@ -45,6 +49,7 @@ class PolarAzimuthal : public UnitSphereDistribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Direction sampled
   Direction sample(uint64_t* seed) const override;
+  double evaluate(Direction u) const override;
 
   // Observing pointers
   Distribution* mu() const { return mu_.get(); }
@@ -71,6 +76,7 @@ class Isotropic : public UnitSphereDistribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled direction
   Direction sample(uint64_t* seed) const override;
+  double evaluate(Direction u) const override;
 };
 
 //==============================================================================
diff --git a/include/openmc/reaction_product.h b/include/openmc/reaction_product.h
index 4fbbc1b626a..32e82524f51 100644
--- a/include/openmc/reaction_product.h
+++ b/include/openmc/reaction_product.h
@@ -49,6 +49,9 @@ class ReactionProduct {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(double E_in, double& E_out, double& mu, uint64_t* seed) const;
 
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const;
+
   ParticleType particle_;      //!< Particle type
   EmissionMode emission_mode_; //!< Emission mode
   double decay_rate_; //!< Decay rate (for delayed neutron precursors) in [1/s]
diff --git a/include/openmc/secondary_correlated.h b/include/openmc/secondary_correlated.h
index 6905c38e369..a4435171860 100644
--- a/include/openmc/secondary_correlated.h
+++ b/include/openmc/secondary_correlated.h
@@ -40,6 +40,8 @@ class CorrelatedAngleEnergy : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
   // energy property
   vector<double>& energy() { return energy_; }
diff --git a/include/openmc/secondary_kalbach.h b/include/openmc/secondary_kalbach.h
index 83806d35248..1c3e80c36ac 100644
--- a/include/openmc/secondary_kalbach.h
+++ b/include/openmc/secondary_kalbach.h
@@ -31,6 +31,8 @@ class KalbachMann : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   //! Outgoing energy/angle at a single incoming energy
diff --git a/include/openmc/secondary_nbody.h b/include/openmc/secondary_nbody.h
index efb4fd75ba1..82789bfa585 100644
--- a/include/openmc/secondary_nbody.h
+++ b/include/openmc/secondary_nbody.h
@@ -27,6 +27,8 @@ class NBodyPhaseSpace : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   int n_bodies_;      //!< Number of particles distributed
diff --git a/include/openmc/secondary_thermal.h b/include/openmc/secondary_thermal.h
index 5b18902afbb..d5b7f0190b6 100644
--- a/include/openmc/secondary_thermal.h
+++ b/include/openmc/secondary_thermal.h
@@ -32,6 +32,8 @@ class CoherentElasticAE : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section
@@ -55,6 +57,8 @@ class IncoherentElasticAE : public AngleEnergy {
   //! \param[inout] seed Pseudorandom number seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   double debye_waller_;
@@ -80,6 +84,8 @@ class IncoherentElasticAEDiscrete : public AngleEnergy {
   //! \param[inout] seed Pseudorandom number seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   const vector<double>& energy_;  //!< Energies at which cosines are tabulated
@@ -106,6 +112,8 @@ class IncoherentInelasticAEDiscrete : public AngleEnergy {
   //! \param[inout] seed Pseudorandom number seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   const vector<double>& energy_; //!< Incident energies
@@ -134,6 +142,8 @@ class IncoherentInelasticAE : public AngleEnergy {
   //! \param[inout] seed Pseudorandom number seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   //! Secondary energy/angle distribution
@@ -169,6 +179,8 @@ class MixedElasticAE : public AngleEnergy {
   //! \param[inout] seed Pseudorandom number seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   CoherentElasticAE coherent_dist_;         //!< Coherent distribution
diff --git a/include/openmc/secondary_uncorrelated.h b/include/openmc/secondary_uncorrelated.h
index 3afa3d9ceb7..d971231b277 100644
--- a/include/openmc/secondary_uncorrelated.h
+++ b/include/openmc/secondary_uncorrelated.h
@@ -31,6 +31,8 @@ class UncorrelatedAngleEnergy : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
   // Accessors
   AngleDistribution& angle() { return angle_; }
diff --git a/include/openmc/thermal.h b/include/openmc/thermal.h
index de0767d0af0..708d7264bc2 100644
--- a/include/openmc/thermal.h
+++ b/include/openmc/thermal.h
@@ -52,6 +52,8 @@ class ThermalData {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(const NuclideMicroXS& micro_xs, double E_in, double* E_out,
     double* mu, uint64_t* seed);
+  double sample_energy_and_pdf(const NuclideMicroXS& micro_xs, double E_in,
+    double mu, double& E_out, uint64_t* seed) const;
 
 private:
   struct Reaction {
diff --git a/src/chain.cpp b/src/chain.cpp
index e279d1f5916..17efdc69d3b 100644
--- a/src/chain.cpp
+++ b/src/chain.cpp
@@ -74,6 +74,13 @@ void DecayPhotonAngleEnergy::sample(
   mu = Uniform(-1., 1.).sample(seed);
 }
 
+double DecayPhotonAngleEnergy::sample_energy_and_pdf(
+  double E_in, double mu, double& E_out, uint64_t* seed) const
+{
+  E_out = photon_energy_->sample(seed);
+  return 0.5;
+}
+
 //==============================================================================
 // Global variables
 //==============================================================================
diff --git a/src/distribution.cpp b/src/distribution.cpp
index ca00cbde4eb..4e09a00b879 100644
--- a/src/distribution.cpp
+++ b/src/distribution.cpp
@@ -12,6 +12,7 @@
 #include "openmc/math_functions.h"
 #include "openmc/random_dist.h"
 #include "openmc/random_lcg.h"
+#include "openmc/search.h"
 #include "openmc/xml_interface.h"
 
 namespace openmc {
@@ -156,6 +157,14 @@ double Uniform::sample(uint64_t* seed) const
   return a_ + prn(seed) * (b_ - a_);
 }
 
+double Uniform::evaluate(double x) const
+{
+  if (x <= b_ && x >= a_)
+    return 1 / (b_ - a_);
+  else
+    return 0;
+}
+
 //==============================================================================
 // PowerLaw implementation
 //==============================================================================
@@ -182,6 +191,21 @@ double PowerLaw::sample(uint64_t* seed) const
   return std::pow(offset_ + prn(seed) * span_, ninv_);
 }
 
+double PowerLaw::evaluate(double x) const
+{
+  // Use accessors
+  double a_val = this->a();
+  double b_val = this->b();
+  double n_val = this->n();
+
+  if (x < a_val || x > b_val) {
+    return 0.0; // outside support
+  }
+
+  return (n_val + 1.0) * std::pow(x, n_val) /
+         (std::pow(b_val, n_val + 1.0) - std::pow(a_val, n_val + 1.0));
+}
+
 //==============================================================================
 // Maxwell implementation
 //==============================================================================
@@ -236,6 +260,13 @@ double Normal::sample(uint64_t* seed) const
   return normal_variate(mean_value_, std_dev_, seed);
 }
 
+double Normal::evaluate(double x) const
+{
+  double exponent = -0.5 * std::pow((x - mean_value_) / std_dev_, 2);
+  double coefficient = 1 / (std_dev_ * std::sqrt(2 * PI));
+  return coefficient * std::exp(exponent);
+}
+
 //==============================================================================
 // Tabular implementation
 //==============================================================================
@@ -356,6 +387,43 @@ double Tabular::sample(uint64_t* seed) const
   }
 }
 
+double Tabular::evaluate(double x) const
+{
+  // get PDF value at x
+
+  int i;
+  std::size_t n = x_.size();
+  if (x < x_[0]) {
+    return 0;
+  } else if (x > x_[n - 1]) {
+    return 0;
+  } else {
+    i = lower_bound_index(x_.begin(), x_.end(), x);
+  }
+
+  // Determine bounding PDF values
+  double x_i = x_[i];
+  double p_i = p_[i];
+
+  switch (interp_) {
+  case Interpolation::histogram:
+    // Histogram interpolation
+    return p_i;
+
+  case Interpolation::lin_lin: {
+    // Linear-linear interpolation
+    double x_i1 = x_[i + 1];
+    double p_i1 = p_[i + 1];
+
+    double m = (p_i1 - p_i) / (x_i1 - x_i);
+    return (m == 0.0) ? p_i : p_i + (x - x_i) * m;
+  }
+
+  default:
+    fatal_error("Unsupported interpolation type in PDF evaluation.");
+  }
+}
+
 //==============================================================================
 // Equiprobable implementation
 //==============================================================================
diff --git a/src/distribution_angle.cpp b/src/distribution_angle.cpp
index 5e92b794a72..7540db0a673 100644
--- a/src/distribution_angle.cpp
+++ b/src/distribution_angle.cpp
@@ -95,4 +95,32 @@ double AngleDistribution::sample(double E, uint64_t* seed) const
   return mu;
 }
 
+double AngleDistribution::evaluate(double E, double mu) const
+{
+  // Determine number of incoming energies
+  auto n = energy_.size();
+
+  // Find energy bin and calculate interpolation factor -- if the energy is
+  // outside the range of the tabulated energies, choose the first or last bins
+  int i;
+  double r;
+  if (E < energy_[0]) {
+    i = 0;
+    r = 0.0;
+  } else if (E > energy_[n - 1]) {
+    i = n - 2;
+    r = 1.0;
+  } else {
+    i = lower_bound_index(energy_.begin(), energy_.end(), E);
+    r = (E - energy_[i]) / (energy_[i + 1] - energy_[i]);
+  }
+
+  double pdf = 0.0;
+  if (r > 0.0)
+    pdf += r * distribution_[i + 1]->evaluate(mu);
+  if (r < 1.0)
+    pdf += (1.0 - r) * distribution_[i]->evaluate(mu);
+  return pdf;
+}
+
 } // namespace openmc
diff --git a/src/distribution_multi.cpp b/src/distribution_multi.cpp
index cdb33adc2ad..edd2dfe9163 100644
--- a/src/distribution_multi.cpp
+++ b/src/distribution_multi.cpp
@@ -44,6 +44,7 @@ UnitSphereDistribution::UnitSphereDistribution(pugi::xml_node node)
       fatal_error("Angular distribution reference direction must have "
                   "three parameters specified.");
     u_ref_ = Direction(u_ref.data());
+    u_ref_ /= u_ref_.norm();
   }
 }
 
@@ -65,6 +66,7 @@ PolarAzimuthal::PolarAzimuthal(pugi::xml_node node)
       fatal_error("Angular distribution reference v direction must have "
                   "three parameters specified.");
     v_ref_ = Direction(v_ref.data());
+    v_ref_ /= v_ref_.norm();
   }
   w_ref_ = u_ref_.cross(v_ref_);
   if (check_for_node(node, "mu")) {
@@ -99,6 +101,13 @@ Direction PolarAzimuthal::sample(uint64_t* seed) const
   return mu * u_ref_ + f * std::cos(phi) * v_ref_ + f * std::sin(phi) * w_ref_;
 }
 
+double PolarAzimuthal::evaluate(Direction u) const
+{
+  double mu = u.dot(u_ref_);
+  double phi = std::acos(u.dot(v_ref_) / std::sqrt(1 - mu * mu));
+  return mu_->evaluate(mu) * phi_->evaluate(phi);
+}
+
 //==============================================================================
 // Isotropic implementation
 //==============================================================================
@@ -116,6 +125,11 @@ Direction Isotropic::sample(uint64_t* seed) const
   return isotropic_direction(seed);
 }
 
+double Isotropic::evaluate(Direction u) const
+{
+  return 1.0 / (4.0 * PI);
+}
+
 //==============================================================================
 // Monodirectional implementation
 //==============================================================================
diff --git a/src/reaction_product.cpp b/src/reaction_product.cpp
index 3ba2c0cfb05..36b92e0b913 100644
--- a/src/reaction_product.cpp
+++ b/src/reaction_product.cpp
@@ -129,4 +129,33 @@ void ReactionProduct::sample(
   }
 }
 
+double ReactionProduct::sample_energy_and_pdf(
+  double E_in, double mu, double& E_out, uint64_t* seed) const
+{
+
+  int distribution_index;
+  auto n = applicability_.size();
+  if (n > 1) {
+    double prob = 0.0;
+    double c = prn(seed);
+    for (int i = 0; i < n; ++i) {
+      // Determine probability that i-th energy distribution is sampled
+      prob += applicability_[i](E_in);
+
+      // If i-th distribution is sampled, sample energy from the distribution
+      if (c <= prob) {
+        // distribution_[i]->sample(E_in, E_out, mu, seed);
+        distribution_index = i;
+        break;
+      }
+    }
+  } else {
+    distribution_index = 0;
+  }
+  // now extract pdf
+
+  AngleEnergy* angleEnergyPtr = distribution_[distribution_index].get();
+  return angleEnergyPtr->sample_energy_and_pdf(E_in, mu, E_out, seed);
+}
+
 } // namespace openmc
diff --git a/src/secondary_correlated.cpp b/src/secondary_correlated.cpp
index 0e4891dd3e7..1e648bce003 100644
--- a/src/secondary_correlated.cpp
+++ b/src/secondary_correlated.cpp
@@ -263,4 +263,116 @@ void CorrelatedAngleEnergy::sample(
   }
 }
 
+double CorrelatedAngleEnergy::sample_energy_and_pdf(
+  double E_in, double mu, double& E_out, uint64_t* seed) const
+{
+  // Find energy bin and calculate interpolation factor -- if the energy is
+  // outside the range of the tabulated energies, choose the first or last bins
+  auto n_energy_in = energy_.size();
+  int i;
+  double r;
+  if (E_in < energy_[0]) {
+    i = 0;
+    r = 0.0;
+  } else if (E_in > energy_[n_energy_in - 1]) {
+    i = n_energy_in - 2;
+    r = 1.0;
+  } else {
+    i = lower_bound_index(energy_.begin(), energy_.end(), E_in);
+    r = (E_in - energy_[i]) / (energy_[i + 1] - energy_[i]);
+  }
+
+  // Sample between the ith and [i+1]th bin
+  int l = r > prn(seed) ? i + 1 : i;
+
+  // Interpolation for energy E1 and EK
+  int n_energy_out = distribution_[i].e_out.size();
+  int n_discrete = distribution_[i].n_discrete;
+  double E_i_1 = distribution_[i].e_out[n_discrete];
+  double E_i_K = distribution_[i].e_out[n_energy_out - 1];
+
+  n_energy_out = distribution_[i + 1].e_out.size();
+  n_discrete = distribution_[i + 1].n_discrete;
+  double E_i1_1 = distribution_[i + 1].e_out[n_discrete];
+  double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1];
+
+  double E_1 = E_i_1 + r * (E_i1_1 - E_i_1);
+  double E_K = E_i_K + r * (E_i1_K - E_i_K);
+
+  // Determine outgoing energy bin
+  n_energy_out = distribution_[l].e_out.size();
+  n_discrete = distribution_[l].n_discrete;
+  double r1 = prn(seed);
+  double c_k = distribution_[l].c[0];
+  int k = 0;
+  int end = n_energy_out - 2;
+
+  // Discrete portion
+  for (int j = 0; j < n_discrete; ++j) {
+    k = j;
+    c_k = distribution_[l].c[k];
+    if (r1 < c_k) {
+      end = j;
+      break;
+    }
+  }
+
+  // Continuous portion
+  double c_k1;
+  for (int j = n_discrete; j < end; ++j) {
+    k = j;
+    c_k1 = distribution_[l].c[k + 1];
+    if (r1 < c_k1)
+      break;
+    k = j + 1;
+    c_k = c_k1;
+  }
+
+  double E_l_k = distribution_[l].e_out[k];
+  double p_l_k = distribution_[l].p[k];
+  if (distribution_[l].interpolation == Interpolation::histogram) {
+    // Histogram interpolation
+    if (p_l_k > 0.0 && k >= n_discrete) {
+      E_out = E_l_k + (r1 - c_k) / p_l_k;
+    } else {
+      E_out = E_l_k;
+    }
+
+  } else if (distribution_[l].interpolation == Interpolation::lin_lin) {
+    // Linear-linear interpolation
+    double E_l_k1 = distribution_[l].e_out[k + 1];
+    double p_l_k1 = distribution_[l].p[k + 1];
+
+    double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k);
+    if (frac == 0.0) {
+      E_out = E_l_k + (r1 - c_k) / p_l_k;
+    } else {
+      E_out =
+        E_l_k +
+        (std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) -
+          p_l_k) /
+          frac;
+    }
+  }
+
+  // Now interpolate between incident energy bins i and i + 1
+  if (k >= n_discrete) {
+    if (l == i) {
+      E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);
+    } else {
+      E_out = E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);
+    }
+  }
+
+  double pdf; // assuming the data in lab frame!
+  if (r1 - c_k < c_k1 - r1 ||
+      distribution_[l].interpolation == Interpolation::histogram) {
+    pdf = distribution_[l].angle[k]->evaluate(mu);
+  } else {
+    pdf = distribution_[l].angle[k + 1]->evaluate(mu);
+  }
+
+  return pdf;
+}
+
 } // namespace openmc
diff --git a/src/secondary_kalbach.cpp b/src/secondary_kalbach.cpp
index 4a7878343c6..b8ebe29b8ca 100644
--- a/src/secondary_kalbach.cpp
+++ b/src/secondary_kalbach.cpp
@@ -237,5 +237,126 @@ void KalbachMann::sample(
     mu = std::log(r1 * std::exp(km_a) + (1.0 - r1) * std::exp(-km_a)) / km_a;
   }
 }
+double KalbachMann::sample_energy_and_pdf(
+  double E_in, double mu, double& E_out, uint64_t* seed) const
+{
+  // Find energy bin and calculate interpolation factor -- if the energy is
+  // outside the range of the tabulated energies, choose the first or last bins
+
+  // double E_out;
+  auto n_energy_in = energy_.size();
+  int i;
+  double r;
+  if (E_in < energy_[0]) {
+    i = 0;
+    r = 0.0;
+  } else if (E_in > energy_[n_energy_in - 1]) {
+    i = n_energy_in - 2;
+    r = 1.0;
+  } else {
+    i = lower_bound_index(energy_.begin(), energy_.end(), E_in);
+    r = (E_in - energy_[i]) / (energy_[i + 1] - energy_[i]);
+  }
+
+  // Sample between the ith and [i+1]th bin
+  int l = r > prn(seed) ? i + 1 : i;
+
+  // Interpolation for energy E1 and EK
+  int n_energy_out = distribution_[i].e_out.size();
+  int n_discrete = distribution_[i].n_discrete;
+  double E_i_1 = distribution_[i].e_out[n_discrete];
+  double E_i_K = distribution_[i].e_out[n_energy_out - 1];
+
+  n_energy_out = distribution_[i + 1].e_out.size();
+  n_discrete = distribution_[i + 1].n_discrete;
+  double E_i1_1 = distribution_[i + 1].e_out[n_discrete];
+  double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1];
+
+  double E_1 = E_i_1 + r * (E_i1_1 - E_i_1);
+  double E_K = E_i_K + r * (E_i1_K - E_i_K);
+
+  // Determine outgoing energy bin
+  n_energy_out = distribution_[l].e_out.size();
+  n_discrete = distribution_[l].n_discrete;
+  double r1 = prn(seed);
+  double c_k = distribution_[l].c[0];
+  int k = 0;
+  int end = n_energy_out - 2;
+
+  // Discrete portion
+  for (int j = 0; j < n_discrete; ++j) {
+    k = j;
+    c_k = distribution_[l].c[k];
+    if (r1 < c_k) {
+      end = j;
+      break;
+    }
+  }
+
+  // Continuous portion
+  double c_k1;
+  for (int j = n_discrete; j < end; ++j) {
+    k = j;
+    c_k1 = distribution_[l].c[k + 1];
+    if (r1 < c_k1)
+      break;
+    k = j + 1;
+    c_k = c_k1;
+  }
+
+  double E_l_k = distribution_[l].e_out[k];
+  double p_l_k = distribution_[l].p[k];
+  double km_r, km_a;
+  if (distribution_[l].interpolation == Interpolation::histogram) {
+    // Histogram interpolation
+    if (p_l_k > 0.0 && k >= n_discrete) {
+      E_out = E_l_k + (r1 - c_k) / p_l_k;
+    } else {
+      E_out = E_l_k;
+    }
+    // Determine Kalbach-Mann parameters
+    km_r = distribution_[l].r[k];
+    km_a = distribution_[l].a[k];
+
+  } else {
+    // Linear-linear interpolation
+    double E_l_k1 = distribution_[l].e_out[k + 1];
+    double p_l_k1 = distribution_[l].p[k + 1];
+
+    double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k);
+    if (frac == 0.0) {
+      E_out = E_l_k + (r1 - c_k) / p_l_k;
+    } else {
+      E_out =
+        E_l_k +
+        (std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) -
+          p_l_k) /
+          frac;
+    }
+    //  Linear-linear interpolation
+    // Determine Kalbach-Mann parameters
+    km_r = distribution_[l].r[k] +
+           (E_out - E_l_k) / (E_l_k1 - E_l_k) *
+             (distribution_[l].r[k + 1] - distribution_[l].r[k]);
+    km_a = distribution_[l].a[k] +
+           (E_out - E_l_k) / (E_l_k1 - E_l_k) *
+             (distribution_[l].a[k + 1] - distribution_[l].a[k]);
+  }
+
+  // Now interpolate between incident energy bins i and i + 1
+  if (k >= n_discrete) {
+    if (l == i) {
+      E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);
+    } else {
+      E_out = E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);
+    }
+  }
+
+  // https://docs.openmc.org/en/latest/methods/neutron_physics.html#equation-KM-pdf-angle
+  double pdf = km_a / (2 * std::sinh(km_a)) *
+               (std::cosh(km_a * mu) + km_r * std::sinh(km_a * mu));
+
+  return pdf;
+}
 
 } // namespace openmc
diff --git a/src/secondary_nbody.cpp b/src/secondary_nbody.cpp
index da0bb81c471..bc170f1d8b8 100644
--- a/src/secondary_nbody.cpp
+++ b/src/secondary_nbody.cpp
@@ -59,7 +59,7 @@ void NBodyPhaseSpace::sample(
         std::log(r5) * std::pow(std::cos(PI / 2.0 * r6), 2);
     break;
   default:
-    throw std::runtime_error {"N-body phase space with >5 bodies."};
+    fatal_error("N-body phase space with >5 bodies.");
   }
 
   // Now determine v and E_out
@@ -67,4 +67,46 @@ void NBodyPhaseSpace::sample(
   E_out = E_max * v;
 }
 
+double NBodyPhaseSpace::sample_energy_and_pdf(
+  double E_in, double mu, double& E_out, uint64_t* seed) const
+{
+  // Determine E_max parameter
+  double Ap = mass_ratio_;
+  double E_max = (Ap - 1.0) / Ap * (A_ / (A_ + 1.0) * E_in + Q_);
+
+  // x is essentially a Maxwellian distribution
+  double x = maxwell_spectrum(1.0, seed);
+
+  double y;
+  double r1, r2, r3, r4, r5, r6;
+  switch (n_bodies_) {
+  case 3:
+    y = maxwell_spectrum(1.0, seed);
+    break;
+  case 4:
+    r1 = prn(seed);
+    r2 = prn(seed);
+    r3 = prn(seed);
+    y = -std::log(r1 * r2 * r3);
+    break;
+  case 5:
+    r1 = prn(seed);
+    r2 = prn(seed);
+    r3 = prn(seed);
+    r4 = prn(seed);
+    r5 = prn(seed);
+    r6 = prn(seed);
+    y = -std::log(r1 * r2 * r3 * r4) -
+        std::log(r5) * std::pow(std::cos(PI / 2.0 * r6), 2);
+    break;
+  default:
+    fatal_error("N-body phase space with >5 bodies.");
+  }
+
+  // Now determine v and E_out
+  double v = x / (x + y);
+  E_out = E_max * v;
+  return 0.5;
+}
+
 } // namespace openmc
diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp
index 8b9e8737c66..0c26e4d2d82 100644
--- a/src/secondary_thermal.cpp
+++ b/src/secondary_thermal.cpp
@@ -3,6 +3,7 @@
 #include "openmc/hdf5_interface.h"
 #include "openmc/random_lcg.h"
 #include "openmc/search.h"
+#include "openmc/vector.h"
 
 #include "xtensor/xview.hpp"
 
@@ -25,6 +26,53 @@ void get_energy_index(
   }
 }
 
+double get_pdf_discrete(
+  const vector<double>& mu, const vector<double>& w, double mu_0)
+{
+  // Make sure mu is in range [-1,1]
+  if (std::abs(mu_0) > 1.0)
+    mu_0 = std::copysign(1.0, mu_0);
+  double a0;
+  double a1;
+  double b0;
+  double b1;
+  int32_t ai = -1;
+  int32_t bi = -1;
+  if (mu_0 > mu[0]) {
+    ai = lower_bound_index(mu.begin(), mu.end(), mu_0);
+    a0 = mu[ai];
+    a1 = (ai > 1) ? mu[ai - 1] : -1.0;
+  } else {
+    a0 = -1.0;
+    a1 = -1.0;
+  }
+  if (mu_0 < mu[mu.size() - 1]) {
+    bi = upper_bound_index(mu.begin(), mu.end(), mu_0);
+    b0 = mu[bi];
+    b1 = (bi < mu.size() - 1) ? mu[bi + 1] : 1.0;
+  } else {
+    b0 = 1.0;
+    b1 = 1.0;
+  }
+
+  //  Calculate Delta_a and Delta_b
+  double delta_a = 0.5 * std::min(b0 - a0, a0 - a1);
+  double delta_b = 0.5 * std::min(b1 - b0, b0 - a0);
+
+  if (mu_0 < a0 + delta_a)
+    return w[ai] / (2.0 * delta_a);
+  else if (mu_0 + delta_b < b0)
+    return w[bi] / (2.0 * delta_b);
+  else
+    return 0.0;
+}
+
+double get_pdf_discrete(const vector<double>& mu, double mu_0)
+{
+  vector<double> w(mu.size(), 1.0 / mu.size());
+  return get_pdf_discrete(mu, w, mu_0);
+}
+
 //==============================================================================
 // CoherentElasticAE implementation
 //==============================================================================
@@ -55,6 +103,39 @@ void CoherentElasticAE::sample(
   mu = 1.0 - 2.0 * energies[k] / E_in;
 }
 
+double CoherentElasticAE::sample_energy_and_pdf(
+  double E_in, double mu, double& E_out, uint64_t* seed) const
+{
+  // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7-1)
+
+  double pdf;
+  E_out = E_in;
+  const auto& energies {xs_.bragg_edges()};
+  const auto& factors = xs_.factors();
+
+  if (E_in < energies.front() || E_in > energies.back()) {
+    return 0;
+  }
+
+  const int i = lower_bound_index(energies.begin(), energies.end(), E_in);
+  vector<double> energies_cut(energies.begin(), energies.begin() + i + 1);
+  vector<double> factors_cut(factors.begin(), factors.begin() + i + 1);
+
+  vector<double> mu_vector_rev;
+  std::transform(energies_cut.begin(), energies_cut.end(),
+    std::back_inserter(mu_vector_rev),
+    [E_in](double Ei) { return 1 - 2 * Ei / E_in; });
+  vector<double> mu_vector(mu_vector_rev.rbegin(), mu_vector_rev.rend());
+
+  auto f = xt::adapt(factors_cut, {
+                                    factors_cut.size(),
+                                  });
+  auto weights = xt::diff(f);
+  weights /= xt::sum(weights);
+  vector<double> w(weights.begin(), weights.end());
+  return get_pdf_discrete(mu_vector, w, mu);
+}
+
 //==============================================================================
 // IncoherentElasticAE implementation
 //==============================================================================
@@ -74,6 +155,19 @@ void IncoherentElasticAE::sample(
   // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7.4)
   E_out = E_in;
 }
+double IncoherentElasticAE::sample_energy_and_pdf(
+  double E_in, double mu, double& E_out, uint64_t* seed) const
+{
+  // Sample angle by inverting the distribution in ENDF-102, Eq. 7.4
+  double c = 2 * E_in * debye_waller_;
+  E_out = E_in;
+
+  double A = c / (1 - std::exp(-2.0 * c)); // normalization factor
+  double pdf = A * std::exp(-c * (1 - mu));
+  return pdf;
+
+  // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7.4)
+}
 
 //==============================================================================
 // IncoherentElasticAEDiscrete implementation
@@ -129,6 +223,27 @@ void IncoherentElasticAEDiscrete::sample(
   E_out = E_in;
 }
 
+double IncoherentElasticAEDiscrete::sample_energy_and_pdf(
+  double E_in, double mu, double& E_out, uint64_t* seed) const
+{
+  // Get index and interpolation factor for elastic grid
+  int i;
+  double f;
+  get_energy_index(energy_, E_in, i, f);
+  // Energy doesn't change in elastic scattering
+  E_out = E_in;
+  int n_mu = mu_out_.shape()[1];
+
+  std::vector<double> mu_vector;
+
+  for (int k = 0; k < n_mu; ++k) {
+    double mu_k = mu_out_(i, k) + f * (mu_out_(i + 1, k) - mu_out_(i, k));
+    mu_vector.push_back(mu_k);
+  }
+
+  return get_pdf_discrete(mu_vector, mu);
+}
+
 //==============================================================================
 // IncoherentInelasticAEDiscrete implementation
 //==============================================================================
@@ -202,6 +317,58 @@ void IncoherentInelasticAEDiscrete::sample(
   mu = (1 - f) * mu_ijk + f * mu_i1jk;
 }
 
+double IncoherentInelasticAEDiscrete::sample_energy_and_pdf(
+  double E_in, double mu, double& E_out, uint64_t* seed) const
+{
+  // Get index and interpolation factor for inelastic grid
+  int i;
+  double f;
+  get_energy_index(energy_, E_in, i, f);
+  int j;
+  int n = energy_out_.shape()[1];
+  if (!skewed_) {
+    // All bins equally likely
+    j = prn(seed) * n;
+  } else {
+    // Distribution skewed away from edge points
+    double r = prn(seed) * (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;
+  int m = mu_out_.shape()[2];
+  std::vector<double> mu_vector;
+
+  for (int k = 0; k < m; ++k) {
+    double mu_ijk = mu_out_(i, j, k);
+    double mu_i1jk = mu_out_(i + 1, j, k);
+    double mu_k = (1 - f) * mu_ijk + f * mu_i1jk;
+    mu_vector.push_back(mu_k);
+  }
+
+  return get_pdf_discrete(mu_vector, mu);
+}
+
 //==============================================================================
 // IncoherentInelasticAE implementation
 //==============================================================================
@@ -331,6 +498,76 @@ void IncoherentInelasticAE::sample(
   mu += std::min(mu - mu_left, mu_right - mu) * (prn(seed) - 0.5);
 }
 
+double IncoherentInelasticAE::sample_energy_and_pdf(
+  double E_in, double mu, double& E_out, uint64_t* seed) const
+{
+  // Get index and interpolation factor for inelastic grid
+  int i;
+  double f;
+  get_energy_index(energy_, E_in, i, f);
+
+  // Pick closer energy based on interpolation factor
+  int l = f > 0.5 ? i + 1 : i;
+
+  // Determine outgoing energy bin
+  // (First reset n_energy_out to the right value)
+  auto n = distribution_[l].n_e_out;
+  double r1 = prn(seed);
+  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;
+  }
+
+  // Adjustment of outgoing energy
+  double E_l = energy_[l];
+  if (E_out < 0.5 * E_l) {
+    E_out *= 2.0 * E_in / E_l - 1.0;
+  } else {
+    E_out += E_in - E_l;
+  }
+
+  // Sample outgoing cosine bin
+  int n_mu = distribution_[l].mu.shape()[1];
+  const auto& mu_l = distribution_[l].mu;
+  f = (r1 - c_j) / (c_j1 - c_j);
+
+  std::vector<double> mu_vector;
+
+  for (int k = 0; k < n_mu; ++k) {
+    double mu_k = mu_l(j, k) + f * (mu_l(j + 1, k) - mu_l(j, k));
+    mu_vector.push_back(mu_k);
+  }
+
+  return get_pdf_discrete(mu_vector, mu);
+}
+
 //==============================================================================
 // MixedElasticAE implementation
 //==============================================================================
@@ -367,4 +604,18 @@ void MixedElasticAE::sample(
   }
 }
 
+double MixedElasticAE::sample_energy_and_pdf(
+  double E_in, double mu, double& E_out, uint64_t* seed) const
+{
+  // Evaluate coherent and incoherent elastic cross sections
+  double xs_coh = coherent_xs_(E_in);
+  double xs_incoh = incoherent_xs_(E_in);
+
+  if (prn(seed) * (xs_coh + xs_incoh) < xs_coh) {
+    return coherent_dist_.sample_energy_and_pdf(E_in, mu, E_out, seed);
+  } else {
+    return incoherent_dist_->sample_energy_and_pdf(E_in, mu, E_out, seed);
+  }
+}
+
 } // namespace openmc
diff --git a/src/secondary_uncorrelated.cpp b/src/secondary_uncorrelated.cpp
index 5cbb76fb9b9..27adf5a172e 100644
--- a/src/secondary_uncorrelated.cpp
+++ b/src/secondary_uncorrelated.cpp
@@ -65,4 +65,19 @@ void UncorrelatedAngleEnergy::sample(
   E_out = energy_->sample(E_in, seed);
 }
 
+double UncorrelatedAngleEnergy::sample_energy_and_pdf(
+  double E_in, double mu, double& E_out, uint64_t* seed) const
+{
+  // Sample outgoing energy
+  E_out = energy_->sample(E_in, seed);
+  double pdf;
+  if (!angle_.empty()) {
+    pdf = angle_.evaluate(E_in, mu);
+  } else {
+    // no angle distribution given => assume isotropic for all energies
+    pdf = 0.5;
+  }
+  return pdf;
+}
+
 } // namespace openmc
diff --git a/src/thermal.cpp b/src/thermal.cpp
index cbe0983ed65..8791db5c745 100644
--- a/src/thermal.cpp
+++ b/src/thermal.cpp
@@ -313,6 +313,20 @@ void ThermalData::sample(const NuclideMicroXS& micro_xs, double E,
     *mu = std::copysign(1.0, *mu);
 }
 
+double ThermalData::sample_energy_and_pdf(const NuclideMicroXS& micro_xs,
+  double E_in, double mu, double& E_out, uint64_t* seed) const
+{
+  AngleEnergy* angleEnergyPtr;
+
+  if (prn(seed) < micro_xs.thermal_elastic / micro_xs.thermal) {
+    angleEnergyPtr = elastic_.distribution.get();
+  } else {
+    angleEnergyPtr = inelastic_.distribution.get();
+  }
+
+  return angleEnergyPtr->sample_energy_and_pdf(E_in, mu, E_out, seed);
+}
+
 void free_memory_thermal()
 {
   data::thermal_scatt.clear();