[Bug] Incorrect interpolation in GateGenericSource::UpdateActivityWithTAC

[Marc-AntoineLeclerc]

GateGenericSource::UpdateActivityWithTAC in opengate/core/opengate_core/opengate_lib/GateGenericSource.cpp does not interpolate the TAC_activities correctly. There are two issues:

1. The index i is not calculated properly. It is one index too high.
2. The interpolation equation is flipped. w1 and w2 should be swapped.

To illustrate the issue, below is a short .cpp script that shows the wrong interpolation and compare it to a correct one.

#include
#include
#include
#include

// ------------------------------------------------------------
// GateGenericSource::UpdateActivityWithTAC (the original code)
// ------------------------------------------------------------
double bad_interpolate(const std::vector& times,
const std::vector& activities,
double time)
{
double fActivity = 0.0;

// Below/above the TAC
if (time < times.front() || time > times.back()) {
    return 0.0;
}

// Search for the time bin
const auto lower = std::lower_bound(times.begin(), times.end(), time);
const auto i = std::distance(times.begin(), lower);

// Last element ?
if (i >= times.size() - 1) {
    return activities.back();
}

// BAD linear interpolation
const double bin_time = times[i + 1] - times[i];
const double w1 = (time - times[i]) / bin_time;
const double w2 = (times[i + 1] - time) / bin_time;

fActivity = activities[i] * w1 + activities[i + 1] * w2;
return fActivity;

}

// ------------------------------------------------------------
// GOOD INTERPOLATION (correct bin selection + correct weights)
// ------------------------------------------------------------
double good_interpolate(const std::vector& times,
const std::vector& activities,
double time)
{
// Below/above the TAC (same behavior)
if (time < times.front() || time > times.back()) {
return 0.0;
}

// Find first element >= time
const auto lower = std::lower_bound(times.begin(), times.end(), time);
auto i = std::distance(times.begin(), lower);

// Exact match or first sample
if (i == 0) {
    return activities[0];
}

// CORRECTION: Move to lower bin edge.
i -= 1;

// Last usable element?
if (i >= times.size() - 1) {
    return activities.back();
}

// CORRECTION: Proper linear interpolation
const double dt = times[i + 1] - times[i];
const double t = (time - times[i]) / dt;

return activities[i] + t * (activities[i + 1] - activities[i]);

}

// ------------------------------------------------------------
// MAIN TEST PROGRAM
// ------------------------------------------------------------
int main() {
// Sample TAC
std::vector times = { 0, 1, 2, 3, 4, 5 };
std::vector activities = { 0, 2, 4, 3, 1, 0 };

std::cout << std::fixed << std::setprecision(4);

std::cout << "time\tgood\tbad\tdifference\n";
std::cout << "--------------------------------------\n";

// Evaluate the old and new interpolation code to reveal the differences
for (double t = 0.0; t <= 5.0; t += 0.1) {
    double g = good_interpolate(times, activities, t);
    double b = bad_interpolate(times, activities, t);
    double diff = b - g;

    std::cout << t << "\t" << g << "\t" << b << "\t" << diff << "\n";
}

return 0;

}

[dsarrut]

oh good catch, thanks ! Can you please propose a PR ?