diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp
index d0d0feb26b..863ffd8b08 100644
--- a/src/distribution_energy.cpp
+++ b/src/distribution_energy.cpp
@@ -182,11 +182,13 @@ double ContinuousTabular::sample(double E) const
 
   // Interpolation for energy E1 and EK
   int n_energy_out = distribution_[i].e_out.size();
-  double E_i_1 = distribution_[i].e_out[0];
+  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();
-  double E_i1_1 = distribution_[i+1].e_out[0];
+  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);
@@ -194,25 +196,38 @@ double ContinuousTabular::sample(double E) const
 
   // Determine outgoing energy bin
   n_energy_out = distribution_[l].e_out.size();
+  n_discrete = distribution_[l].n_discrete;
   double r1 = prn();
   double c_k = distribution_[l].c[0];
-  double c_k1;
-  int k;
-  for (k = 0; k < n_energy_out - 2; ++k) {
-    c_k1 = distribution_[l].c[k+1];
-    if (r1 < c_k1) break;
-    c_k = c_k1;
+  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;
+    }
   }
 
-  // Check to make sure 1 <= k <= NP - 1
-  k = std::max(0, std::min(k, n_energy_out - 2));
+  // 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 E_out;
   if (distribution_[l].interpolation == Interpolation::histogram) {
     // Histogram interpolation
-    if (p_l_k > 0.0 && k >= distribution_[l].n_discrete) {
+    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;
@@ -233,7 +248,7 @@ double ContinuousTabular::sample(double E) const
   }
 
   // Now interpolate between incident energy bins i and i + 1
-  if (!histogram_interp && n_energy_out > 1) {
+  if (!histogram_interp && n_energy_out > 1 && k >= n_discrete) {
     if (l == i) {
       return E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1);
     } else {
diff --git a/src/photon_header.F90 b/src/photon_header.F90
index bec559b995..a284c7cc02 100644
--- a/src/photon_header.F90
+++ b/src/photon_header.F90
@@ -163,14 +163,18 @@ contains
     call this % coherent_int_form_factor % from_hdf5(dset_id)
     call close_dataset(dset_id)
 
-    dset_id = open_dataset(rgroup, 'anomalous_real')
-    call this % coherent_anomalous_real % from_hdf5(dset_id)
-    call close_dataset(dset_id)
+    if (object_exists(group_id, 'anomalous_real')) then
+      dset_id = open_dataset(rgroup, 'anomalous_real')
+      call this % coherent_anomalous_real % from_hdf5(dset_id)
+      call close_dataset(dset_id)
+    end if
 
-    dset_id = open_dataset(rgroup, 'anomalous_imag')
-    call this % coherent_anomalous_imag % from_hdf5(dset_id)
-    call close_dataset(dset_id)
-    call close_group(rgroup)
+    if (object_exists(group_id, 'anomalous_imag')) then
+      dset_id = open_dataset(rgroup, 'anomalous_imag')
+      call this % coherent_anomalous_imag % from_hdf5(dset_id)
+      call close_dataset(dset_id)
+      call close_group(rgroup)
+    end if
 
     ! Read incoherent scattering
     rgroup = open_group(group_id, 'incoherent')
diff --git a/src/secondary_correlated.cpp b/src/secondary_correlated.cpp
index ec741a1e60..e238634c3d 100644
--- a/src/secondary_correlated.cpp
+++ b/src/secondary_correlated.cpp
@@ -177,11 +177,13 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu) const
 
   // Interpolation for energy E1 and EK
   int n_energy_out = distribution_[i].e_out.size();
-  double E_i_1 = distribution_[i].e_out[0];
+  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();
-  double E_i1_1 = distribution_[i+1].e_out[0];
+  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);
@@ -189,24 +191,37 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu) const
 
   // Determine outgoing energy bin
   n_energy_out = distribution_[l].e_out.size();
+  n_discrete = distribution_[l].n_discrete;
   double r1 = prn();
   double c_k = distribution_[l].c[0];
-  double c_k1;
-  int k;
-  for (k = 0; k < n_energy_out - 2; ++k) {
-    c_k1 = distribution_[l].c[k+1];
-    if (r1 < c_k1) break;
-    c_k = c_k1;
+  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;
+    }
   }
 
-  // Check to make sure 1 <= k <= NP - 1
-  k = std::max(0, std::min(k, n_energy_out - 2));
+  // 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 >= distribution_[l].n_discrete) {
+    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;
@@ -228,10 +243,12 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu) const
   }
 
   // Now interpolate between incident energy bins i and i + 1
-  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);
+  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);
+    }
   }
 
   // Find correlated angular distribution for closest outgoing energy bin
diff --git a/src/secondary_kalbach.cpp b/src/secondary_kalbach.cpp
index 0779230688..cec5344def 100644
--- a/src/secondary_kalbach.cpp
+++ b/src/secondary_kalbach.cpp
@@ -144,11 +144,13 @@ void KalbachMann::sample(double E_in, double& E_out, double& mu) const
 
   // Interpolation for energy E1 and EK
   int n_energy_out = distribution_[i].e_out.size();
-  double E_i_1 = distribution_[i].e_out[0];
+  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();
-  double E_i1_1 = distribution_[i+1].e_out[0];
+  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);
@@ -156,25 +158,38 @@ void KalbachMann::sample(double E_in, double& E_out, double& mu) const
 
   // Determine outgoing energy bin
   n_energy_out = distribution_[l].e_out.size();
+  n_discrete = distribution_[l].n_discrete;
   double r1 = prn();
   double c_k = distribution_[l].c[0];
-  double c_k1;
-  int k;
-  for (k = 0; k < n_energy_out - 2; ++k) {
-    c_k1 = distribution_[l].c[k+1];
-    if (r1 < c_k1) break;
-    c_k = c_k1;
+  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;
+    }
   }
 
-  // Check to make sure 1 <= k <= NP - 1
-  k = std::max(0, std::min(k, n_energy_out - 2));
+  // 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 >= distribution_[l].n_discrete) {
+    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;
@@ -205,10 +220,12 @@ void KalbachMann::sample(double E_in, double& E_out, double& mu) const
   }
 
   // Now interpolate between incident energy bins i and i + 1
-  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);
+  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);
+    }
   }
 
   // Sampled correlated angle from Kalbach-Mann parameters