% check symbolic integration of tricubic
clear
syms a00 a01 a02 a03
syms a10 a11 a12 a13
syms a20 a21 a22 a23
syms a30 a31 a32 a33
syms x y;
syms sx sy;
syms f00 f01 f10 f11
syms fx00 fx01 fx10 fx11
syms fy00 fy01 fy10 fy11
syms fxy00 fxy01 fxy10 fxy11

%% tricubic polynomial
F = [f00 f01 fy00 fy01;
         f10 f11 fy10 fy11;
         fx00 fx01 fxy00 fxy01;
         fx10 fx11 fxy10 fxy11]; 
M = [1 0 0 0; 0 0 1 0; -3 3 -2 -1; 2 -2 1 1];
A = M*F*M.';
%A = [a00 a01 a02 a03; a10 a11 a12 a13; a20 a21 a22 a23; a30 a31 a32 a33];
xhat = x/sx;
yhat = y/sy;

f = sym(0);
for i = 1 : 4
    for j = 1 : 4
        f = f + A(i,j) * xhat^(i-1) * yhat^(j-1);
    end
end

%% integrate on interval 0,sx, 0,sy
fint = int(int(f, x, 0, sx), y, 0, sy);
fint = simplify(fint);

syms dx dy;
fint_x = int(int(f, x, dx, sx+dx), y, 0, sy);
fint_x = subs(diff(fint_x, dx), dx, 0);
fint_y = int(int(f, x, 0, sx), y, dy, dy+sy);
fint_y = subs(diff(fint_y, dy), dy, 0);


simplify(fint/sx/sy)
simplify(fint_x/sy)
simplify(fint_y/sx)

%%
[c,ai]=coeffs(fint,A);
[cx, aix] = coeffs(fint_x,A);
[cy, aiy] = coeffs(fint_y,A);

pretty(fint)
disp([ai; c/sx/sy].');

pretty(fint_x);
disp([aix; cx/sy].');

pretty(fint_y);
disp([aiy; cy/sx].');

%% check integration coefficients
c_ref=sx*sy*[1, 1/2, 1/3, 1/4, 1/2, 1/4, 1/6, 1/8, 1/3, 1/6, 1/9, 1/12, 1/4, 1/8, 1/12, 1/16].';
cx_ref=sy*[0; 1; 1; 1; 0; 1/2; 1/2; 1/2; 0; 1/3; 1/3; 1/3; 0; 1/4; 1/4; 1/4];
cy_ref=sx*[0; 0; 0; 0; 1; 1/2; 1/3; 1/4; 1; 1/2; 1/3; 1/4; 1; 1/2; 1/3; 1/4];

simplify(c_ref.'*A(:) - fint)
simplify(cx_ref.'*A(:) - fint_x)
simplify(cy_ref.'*A(:) - fint_y)