Linear interpolation

Linear interpolation is implemented in GammaLib through the GNodeArray class. This class contains a collection of nodes $x_i$ that may be used to describe a functional relation $y_i=f(x_i)$. The following code illustrates how the GNodeArray class is used (see examples/cpp/interpolate/interpolate.cpp for the source code):

C++

double x_i[] = {1.0, 4.0, 6.0};
double y_i[] = {8.0, 7.0, 2.0};
GNodeArray nodes(3, x_i);
for (double x = 0; x < 10.0; x += 0.5) {
    nodes.set_value(x);
    double y = y_i[nodes.inx_left()]  * nodes.wgt_left() + y_i[nodes.inx_right()] * nodes.wgt_right();
    std::cout << "x=" << x << " : y=" << y << std::endl;
}

Python

x_i   = [1.0, 4.0, 6.0]
y_i   = [8.0, 7.0, 2.0]
nodes = gammalib.GNodeArray(x_i)
for x in range(20):
    nodes.set_value(x*0.5)
    y = y_i[nodes.inx_left()]  * nodes.wgt_left() + y_i[nodes.inx_right()] * nodes.wgt_right()
    print('x=%f : y=%f' % (x, y))

In line 1, the nodes $x_i$ at which the function values $y_i$ are given are declared, the actual function values $y_i$ are declared in line 2. In line 3, a node array is constructed from the node values. Note that the actual function values are not part of the node array, only the node values are in fact used by the GNodeArray class.

In lines 4-8, the function is interpolated at a number of values in the interval $[0,10[$. In line 5, the $x$ value is set at which the interpolation should be done. The interpolation is then done in line 6 using the formula

$$y=y_{i_{\mathrm{l}\mathrm{e}\mathrm{f}\mathrm{t}}}\times w_{i_{\mathrm{l}\mathrm{e}\mathrm{f}\mathrm{t}}}+y_{i_{\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}}\times w_{i_{\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}}$$

where $i_{\mathrm{l}\mathrm{e}\mathrm{f}\mathrm{t}}$ and $i_{\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}$ are the node indices that encompass the $x$ value, and $w_{i_{\mathrm{l}\mathrm{e}\mathrm{f}\mathrm{t}}}$ and $w_{i_{\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}}$ are the weights with which the function values $y_{i_{\mathrm{l}\mathrm{e}\mathrm{f}\mathrm{t}}}$ and $y_{i_{\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}}$ need to be multiplied to obtain the interpolated value $y$. Note that

$$w_{i_{\mathrm{l}\mathrm{e}\mathrm{f}\mathrm{t}}}+w_{i_{\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}}=1$$

The method also works for extrapolation. For $x<x_0$, $i_{\mathrm{l}\mathrm{e}\mathrm{f}\mathrm{t}}=0$ and $i_{\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}=1$, while for $x>x_{i_{\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{t}}}$, $i_{\mathrm{l}\mathrm{e}\mathrm{f}\mathrm{t}}=i_{\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{t}}-1$ and $i_{\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}=i_{\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{t}}$ (where $i_{\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{t}}$ is the index of the last node, which is $2$ in the example above). The weights are set so that $y$ is extrapolated linearly.

It is obvious that GNodeArray needs at least 2 node values to operate.