diff --git a/NuclearData/Reactions/thermalScattReactionCE/thermalScatterElastic_class.f90 b/NuclearData/Reactions/thermalScattReactionCE/thermalScatterElastic_class.f90 index 5dc831a2b..7ee1d231b 100644 --- a/NuclearData/Reactions/thermalScattReactionCE/thermalScatterElastic_class.f90 +++ b/NuclearData/Reactions/thermalScattReactionCE/thermalScatterElastic_class.f90 @@ -234,26 +234,24 @@ subroutine sampleOut(self, mu, phi, E_out, E_in, rand, lambda) ! Get energy indexes if (E_in > self % eIn(size(self % eIn))) then - l1 = size(self % eIn) - 1 + l2 = size(self % eIn) else - l1 = binarySearch(self % eIn, E_in) + l2 = binarySearch(self % eIn, E_in) end if - l2 = l1 + 1 - ! Sample a Bragg edge at a smaller energy than E_in - prob = self % pValues(1:l2) - prob(1) = ZERO - r = rand % get() * prob(l2) + prob = [ZERO, self % pValues(1:l2)] + r = rand % get() * prob(l2 + 1) k = binarySearch(prob, r) E2 = self % eIn(k) ! Compute angle - mu = 1 - TWO * E2/E_in + mu = ONE - TWO * E2/E_in if (abs(mu) > ONE) then call fatalError(Here,'Failed to get angle'//numToChar(mu)) end if + end if ! Sample phi diff --git a/NuclearData/ceNeutronData/aceDatabase/thermalScatteringData_class.f90 b/NuclearData/ceNeutronData/aceDatabase/thermalScatteringData_class.f90 index 7132c39e4..79f52276f 100644 --- a/NuclearData/ceNeutronData/aceDatabase/thermalScatteringData_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/thermalScatteringData_class.f90 @@ -203,13 +203,13 @@ pure function getInelXS(self, E) result(val) f = (E - E_low) / (E_top - E_low) end associate - val = self % inelastic % xs(idx + 1) * f + (ONE-f) * self % inelastic % xs(idx) + val = self % inelastic % xs(idx + 1) * f + (ONE - f) * self % inelastic % xs(idx) end function getInelXS !! !! Calculate elastic cross section, returns cross section values - !! It's called when S(a,b) tablea are switched on, even if there's no thermal + !! It's called when S(a,b) tables are switched on, even if there's no thermal !! elastic scattering !! !! Args: @@ -255,7 +255,7 @@ pure function getElXS(self, E) result(val) E_low => self % elastic % eGrid(idx)) f = (E - E_low) / (E_top - E_low) end associate - val = self % elastic % xs(idx + 1)*f + (ONE-f) * self % elastic % xs(idx) + val = self % elastic % xs(idx + 1) * f + (ONE - f) * self % elastic % xs(idx) end if end if