diff --git a/src/thermal.cpp b/src/thermal.cpp
index 48479bbd3a..4168d53200 100644
--- a/src/thermal.cpp
+++ b/src/thermal.cpp
@@ -490,7 +490,7 @@ ThermalData::sample(const NuclideMicroXS* micro_xs, double E,
       double E_i_J = inelastic_data_[i].e_out(n - 1);
 
       // Determine endpoints on grid i + 1
-      n = inelastic_data_[i].e_out.size();
+      n = inelastic_data_[i + 1].e_out.size();
       double E_i1_1 = inelastic_data_[i + 1].e_out(0);
       double E_i1_J = inelastic_data_[i + 1].e_out(n - 1);
 
@@ -504,7 +504,7 @@ ThermalData::sample(const NuclideMicroXS* micro_xs, double E,
       double c_j = inelastic_data_[l].e_out_cdf[0];
       double c_j1;
       std::size_t j;
-      for (j = 0; j < n - 2; ++j) {
+      for (j = 0; j < n - 1; ++j) {
         c_j1 = inelastic_data_[l].e_out_cdf[j + 1];
         if (r1 < c_j1) break;
         c_j = c_j1;
@@ -532,9 +532,9 @@ ThermalData::sample(const NuclideMicroXS* micro_xs, double E,
 
       // Now interpolate between incident energy bins i and i + 1
       if (l == i) {
-        *E_out = E_1 + (E - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1);
+        *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 - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1);
+        *E_out = E_1 + (*E_out - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1);
       }
 
       // Sample outgoing cosine bin
@@ -560,9 +560,9 @@ ThermalData::sample(const NuclideMicroXS* micro_xs, double E,
       // Determine (k+1)th mu value
       double mu_right;
       if (k == n_inelastic_mu_ - 1) {
-        mu_right = 1.0 - *mu;
+        mu_right = 1.0;
       } else {
-        mu_right = mu_l(j, k+1) + f*(mu_l(j+1, k+1) - mu_l(j, k+1)) - *mu;
+        mu_right = mu_l(j, k+1) + f*(mu_l(j+1, k+1) - mu_l(j, k+1));
       }
 
       // Smear angle