Fix several bugs in continuous S(a,b) treatment

This commit is contained in:
Paul Romano 2018-10-02 08:46:37 -05:00
parent 04000c5839
commit 3b7a157719

View file

@ -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