anchor.diff:      ! Find energy index on unionized grid
       call find_energy_index(p)
   
  +    ! Check if there's an S(a,b) table for this material
  +    if (mat % has_sab_table) then
  +       sab_threshold = sab_tables(mat % sab_table) % threshold_inelastic
  +    else
  +       sab_threshold = ZERO
  +    end if
  +
       ! Add contribution from each nuclide in material
       do i = 1, mat % n_nuclides
          ! Determine microscopic cross sections for this nuclide
          index_nuclide = mat % nuclide(i)
  -       call calculate_nuclide_xs(p, index_nuclide)
  +
  +       ! Determine whether to use S(a,b) based on energy of particle and whether
  +       ! the nuclide matches the S(a,b) table
  +       if (p % E < sab_threshold .and. i == mat % sab_nuclide) then
  +          index_sab = mat % sab_table
  +       else
  +          index_sab = 0
  +       end if
  +
  +       ! Calculate microscopic cross section for this nuclide
  +       call calculate_nuclide_xs(p, index_nuclide, index_sab)
   
          ! Copy atom density of nuclide in material
          atom_density = mat % atom_density(i)
  
positive.code: 
      type(Particle), pointer :: p
  
      integer                 :: i             ! loop index over nuclides
      integer                 :: index_nuclide ! index into nuclides array
      real(8)                 :: atom_density  ! atom density of a nuclide
      type(Material), pointer :: mat           ! current material
  
      ! If the material is the same as the last material and the energy of the
      ! particle hasn't changed, we don't need to lookup cross sections again.
  
negative.code:             (ONE-f) * nuc % fission(IE) + f * nuc % fission(IE+1)
  
         ! Calculate microscopic nuclide nu-fission cross section
         nu = nu_total(nuc, p % E)
         micro_xs(i) % nu_fission = nu * micro_xs(i) % fission
      end if
  
      ! If there is S(a,b) data for this nuclide, we need to do a few
      ! things. Since the total cross section was based on non-S(a,b) data, we
      ! need to correct it by subtracting the non-S(a,b) elastic cross section and
      ! then add back in the calculated S(a,b) elastic+inelastic cross section.
  
      if (index_sab > 0) then
         ! Get pointer to S(a,b) table
         sab => sab_tables(index_sab)
  
         ! Get index and interpolation factor for inelastic grid
         if (p%E < sab % inelastic_e_in(1)) then
            IE_sab = 1
            f_sab = ZERO
         else
            IE_sab = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, p%E)
            f_sab = (p%E - sab%inelastic_e_in(IE_sab)) / & 
                 (sab%inelastic_e_in(IE_sab+1) - sab%inelastic_e_in(IE_sab))
         end if
  
         ! Calculate S(a,b) inelastic scattering cross section
         inelastic = (ONE-f_sab) * sab % inelastic_sigma(IE_sab) + f_sab * &
              sab % inelastic_sigma(IE_sab + 1)
  
         ! Check for elastic data
         if (p % E < sab % threshold_elastic) then
            ! Get index on elastic P grid
            IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E)
  
            ! Determine treatment of elastic scattering -- if NXS(5) = 4, then the
            ! cross section is given as P/E whereas if NXS(5) /= 4, then it is P.
            if (sab % elastic_type == 4) then
               elastic = sab % elastic_P(IE_sab) / p % E
            else
               ! Get interpolation factor for elastic grid
               f_sab = (p%E - sab%elastic_e_in(IE_sab))/(sab%elastic_e_in(IE_sab+1) - &
                    sab%elastic_e_in(IE_sab))
  
               ! Calculate S(a,b) elastic scattering cross section
               elastic = (ONE-f_sab) * sab % elastic_P(IE_sab) + f_sab * &
                    sab % elastic_P(IE_sab + 1)
            end if
         else
            elastic = ZERO
         end if
  
         ! Correct total and elastic cross sections
         micro_xs(i) % total = micro_xs(i) % total - micro_xs(i) % elastic &
              + inelastic + elastic
         micro_xs(i) % elastic = inelastic + elastic
      end if
  
    end subroutine calculate_nuclide_xs
  