refactor: Fix definition and usage of alpha (and delta)

Signed-off-by: Sietze van Buuren <[email protected]>

Author: swvanbuuren

docs/theory/cubic.md

@@ -70,8 +70,8 @@ in which the expressions $\delta^{(j,i)}_k$ read
 |   | $i=0$ | $i=1$ | $i=2$ | $i=3$ |
 |--:|:-----:|:-----:|:-----:|:-----:|
 | $j=0$ | ${}^1x_k^2({}^1x_k-3{}^0x_k)$ | $+6{}^0x_k{}^1x_k$ | $-3({}^0x_k+{}^1x_k)$ | $+2$ |
-| $j=1$ | ${}^0x_k^2(3{}^1x_k-{}^0x_k)$ | $-6{}^0x_k{}^1x_k$ | $+3({}^0x_k+{}^1x_k)$ | $-2$ |
-| $j=2$ | $-{}^0x_k{}^1x_k^2$ | ${}^1x_k(2{}^0x_k+{}^1x_k)$ | $-({}^0x_k+2{}^1x_k)$ | $+1$ |
+| $j=1$ | $-{}^0x_k{}^1x_k^2$ | ${}^1x_k(2{}^0x_k+{}^1x_k)$ | $-({}^0x_k+2{}^1x_k)$ | $+1$ |
+| $j=2$ | ${}^0x_k^2(3{}^1x_k-{}^0x_k)$ | $-6{}^0x_k{}^1x_k$ | $+3({}^0x_k+{}^1x_k)$ | $-2$ |
 | $j=3$ | $-{}^1x_k{}^0x_k^2$ | ${}^0x_k({}^0x_k+2{}^1x_k)$ | $-(2{}^0x_k+{}^1x_k)$ | $+1$ |
 
 Note that expression for $\alpha_i^{3,j}$ can be written even shorter by introducing a binary number $ij$, i.e. $00, 01, 10, 11 = 0, 1, 2, 3$:

include/cubic_spline.hpp

@@ -154,8 +154,8 @@ class CubicCell2D
         const T x02 = x0*x0;
         const T x12 = x1*x1;
         return {{{x12*(x1 - 3.0*x0), +6.0*x0*x1,           -3.0*(x0 + x1), +2.0},
-                 {x02*(3.0*x1 - x0), -6.0*x0*x1,           +3.0*(x0 + x1), -2.0},
                  {-x0*x12,           x1*(2.0*x0 + x1),     -(x0 + 2.0*x1), +1.0},
+                 {x02*(3.0*x1 - x0), -6.0*x0*x1,           +3.0*(x0 + x1), -2.0},
                  {-x1*x02,           x0*(x0 + 2.0*x1),     -(2.0*x0 + x1), +1.0}}};
     }
 
@@ -180,8 +180,8 @@ class CubicCell2D
                     {
                         std::size_t j1 = (j >> 0) & 1;
                         std::size_t j2 = (j >> 1) & 1;
-                        std::size_t alpha_index1 = (j1 << 1) | i1;
-                        std::size_t alpha_index2 = (j2 << 1) | i2;
+                        std::size_t alpha_index1 = (i1 << 1) | j1;
+                        std::size_t alpha_index2 = (i2 << 1) | j2;
                         double prod_h = (j1 ? h0 : 1.0)*(j2 ? h1 : 1.0);
                         coeffs[k][l] += prod_h*F(i1,i2,j1,j2)*alphas[0][alpha_index1][k]*alphas[1][alpha_index2][l];
                     }