Simplify get_energy_index and use it in place of duplicated lookup logic

[GuySten]

Simplify get_energy_index and use it in place of duplicated lookup logic

Description

get_energy_index returns an index into an incident energy grid along with an
interpolation factor. Its current form nests the whole body inside a condition
and decides i and f in two separate places:

i = 0;
f = 0.0;
if (E >= energies.front()) {
  i = lower_bound_index(energies.begin(), energies.end(), E);
  if (i + 1 < energies.size())
    f = (E - energies[i]) / (energies[i + 1] - energies[i]);
}

Cleaned up here:

• The out-of-range case becomes an early return instead of wrapping the body.
• i + 1 < energies.size() compares an int against a std::size_t. Comparing
  against a local const int n instead removes the signed/unsigned mismatch.
• The branch for an energy above the top of the grid sets i and f together,
  rather than leaving i at n - 1 and only suppressing f.
  Every caller of this function interpolates between the distributions at i and
  i + 1, so keeping i inside the topmost interval and setting f = 1.0 selects
  the distribution at the highest tabulated energy. That lets two other pieces of
  code lose logic they no longer need.

ContinuousTabular::sample

This function did the same lookup inline, with exactly the clamping behaviour
get_energy_index now has:

auto n_energy_in = energy_.size();
if (E < energy_[0]) {
  i = 0;  r = 0.0;
} else if (E > energy_[n_energy_in - 1]) {
  i = n_energy_in - 2;  r = 1.0;
} else {
  i = lower_bound_index(energy_.begin(), energy_.end(), E);
  r = (E - energy_[i]) / (energy_[i + 1] - energy_[i]);
}

Twelve lines replaced by a call to the shared helper. The two implementations
agree case by case, including E exactly at the first or last grid point, so
this removes a duplicate rather than changing anything. search.h becomes
unused in that translation unit and is dropped.

AngleDistribution::evaluate

double pdf = 0.0;
if (r > 0.0)
  pdf += r * distribution_[i + 1]->evaluate(mu);
if (r < 1.0)
  pdf += (1.0 - r) * distribution_[i]->evaluate(mu);
return pdf;

The r > 0.0 test was doing double duty: skipping a zero-weight term, and
keeping distribution_[i + 1] in range when i came back as n - 1. With the
index clamped, only the first job remains, and it no longer earns two branches
in an expression that is a plain linear blend:

return r * distribution_[i + 1]->evaluate(mu) +
       (1.0 - r) * distribution_[i]->evaluate(mu);

Behaviour

Results are unchanged. Because f moves from 0.0 to 1.0 at the same time as
i moves from n - 1 to n - 2, every caller lands on the same distribution:

Caller pattern	Before (i = n-1, f = 0)	After (i = n-2, f = 1)
l = f > prn(seed) ? i + 1 : i	n - 1	n - 1
l = f > 0.5 ? i + 1 : i	n - 1	n - 1
X_i + f * (X_{i+1} - X_i)	X_{n-1}	X_{n-1}
if (r > 0.0) ... [i+1]; if (r < 1.0) ... [i]	[n-1]	[n-1]

The number of prn(seed) calls is unchanged on every path, so the random number
stream is untouched. Energies at or below the top of the grid take exactly the
same path as before.

Two consequences worth stating explicitly:

With i left at n - 1, callers reading element i + 1 were reading one past
the end. In secondary_correlated.cpp and secondary_kalbach.cpp that is
distribution_[i + 1].e_out, forming a reference to an object outside the
vector. The value was always multiplied by f = 0 and so never affected
results, but the access itself was invalid. Keeping i in range removes it.

In AngleDistribution::evaluate, dropping the guards means both terms are now
always computed. For r exactly 0.0 or 1.0 the zero-weight term is
evaluated and multiplied by zero, where before it was skipped. This is the same
exposure the function already has for every interior energy, where both terms
are computed regardless, so it introduces no new failure mode — but it is a
behaviour difference if a Tabular::evaluate could ever return a non-finite
value.

The n < 2 guard in get_energy_index is defensive: without it the clamp would
produce i = -1 on a single-point grid. Such a grid has no interval to
interpolate over and callers would still be out of range, so this only prevents
the new branch from making that case worse.

Testing

No behaviour change, so existing regression tests cover this. Happy to add a
unit test for get_energy_index in tests/cpp_unit_tests/test_math.cpp if
reviewers would like one.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)