I work as a freelance developer, and my current client has given me a program to complete.
I need to code a bicubic spline interpolation in C++.
I have to compare my values to those of an internal Python program that uses the scipy.RectBivariateSpline function.

Despite thoroughly reading up on the topic on several websites, I am unable to manage it.

Here are the details of my problem:
I receive a file that contains power values in the form of a matrix.
Each value is associated with an x-coordinate and a y-coordinate.
The row index indicates the x-coordinate, and the column index indicates the y-coordinate.
The matrix is not large (maximum 600×600).
The step between each x/y-coordinate value is regular.
Example:
x-coordinates: -5.0, -4.0, -3.0, -2.0, etc.
y-coordinates: -6.0, -4.0, -2.0, 0.0, 2.0, etc.

The calculation to be performed is as follows:
I have a point P(x, y), and I need to interpolate its power value based on the power values of the matrix cells in which it is located.

Localization:
By searching the matrix for the first x-coordinate index greater than the x-coordinate of P, I can deduce the pair of indices:
absSupIdx (upper x-coordinate index)
absInfIdx (lower x-coordinate index = absSupIdx – 1).
Similarly for the y-coordinate:
ordSupIdx (upper y-coordinate index)
ordInfIdx (lower y-coordinate index = ordSupIdx – 1).

I thus obtain 2 pairs of indices: (absInfIdx, absSupIdx) and (ordInfIdx, ordSupIdx).
To get a 4×4 matrix necessary for bicubic spline interpolation, I extend to the neighboring cells of the previous 4 cells to obtain a 4×4 matrix.
For each cell in this matrix, I can associate an x-coordinate, a y-coordinate, and a power value.

This is where it gets complicated!
I can’t figure out how to determine the coefficient matrix with the function values, first and second derivatives, etc.
If I can correctly determine this matrix, I can then invert it, which I know how to do, and calculate the interpolated value.
If the inverse coefficient matrix can be defined directly, it saves me from having to code the inversion and the process will be faster.

If we start with the following example, how can we interpolate the power value at the given point P(7.099, -1.8)?
Grey cells are row/column indices and blue cells are x-coordinates and y-coordinates.
Initial matrix (simplified to values close to P)

Based on my localization algorithm, point P is located within the following indices: (330, 331) for the x-coordinate and (210, 211) for the y-coordinate.
This gives the cells outlined in red and the expanded matrix outlined in blue in the previous image.

Could someone explain to me (in detail if possible with the values from the example) the calculations to perform to define the coefficients ai,j in the following formula (from the site https://en.wikipedia.org/wiki/Bicubic_interpolation):

I also wonder if we should normalize the x-axis and y-axis values so that the cell in the extended matrix at indices (0, 0) corresponds to x-axis/y-axis = (0.0, 0.0)?
All other values are also normalized as well as the values of P.

For information, the value interpolated for this example by scipy.RectBivariateSpline is 4831,71270119

Any help is welcome and I thank you in advance if you could direct me to the right methodology because there I am really stuck on the subject!

Regards.