Fix bug in S(a,b) coherent elastic scattering angle sampling

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

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