function [fint, fint_x, fint_y] = integrate_bicubic(phi, dphi, d2phi, x_lim, y_lim)
    % integrate a function phi(r) over the interval
    %       x_lim(1) < x < x_lim(2)
    %       y_lim(1) < y < y_lim(2)
    % where r^2 = x^2 + y^2
    % using bicubic interpolation
    %
    % https://en.wikipedia.org/wiki/Bicubic_interpolation
    % 
    % see also maths_bicubic.m
    
    %% get the tricubic polynomial coefficients on the rescaled interval
    % f(x,y) -> f(x^, y^)
    % 0 < x^ < 1
    % 0 < y^ < 1
    
    % evaluate the matrix F
    r00 = sqrt(x_lim(1)^2 + y_lim(1)^2);
    r01 = sqrt(x_lim(1)^2 + y_lim(2)^2);
    r10 = sqrt(x_lim(2)^2 + y_lim(1)^2);
    r11 = sqrt(x_lim(2)^2 + y_lim(2)^2);
    f00 = phi(r00);
    f01 = phi(r01);
    f10 = phi(r10);
    f11 = phi(r11);

    fx00= dphi(r00) * x_lim(1)/r00;
    fx01= dphi(r01) * x_lim(1)/r01;
    fx10= dphi(r10) * x_lim(2)/r10;
    fx11= dphi(r11) * x_lim(2)/r11;
    fy00= dphi(r00) * y_lim(1)/r00;
    fy01= dphi(r01) * y_lim(2)/r01;
    fy10= dphi(r10) * y_lim(1)/r10;
    fy11= dphi(r11) * y_lim(2)/r11;
    fxy00= d2phi(r00) * x_lim(1)*y_lim(1)/r00^2;
    fxy01= d2phi(r01) * x_lim(1)*y_lim(2)/r01^2;
    fxy10= d2phi(r10) * x_lim(2)*y_lim(1)/r10^2;
    fxy11= d2phi(r11) * x_lim(2)*y_lim(2)/r11^2;
    
    % scale factor for integral
    sx= x_lim(2) - x_lim(1);
    sy= y_lim(2) - y_lim(1);
    
    fint = (sx*sy*(36*f00 + 36*f01 + 36*f10 + 36*f11 + 6*fx00 + 6*fx01 - 6*fx10 - 6*fx11 + fxy00 - fxy01 - fxy10 + fxy11 + 6*fy00 - 6*fy01 + 6*fy10 - 6*fy11))/144;
    fint_x = sy*(f10/2 - f01/2 - f00/2 + f11/2 - fy00/12 + fy01/12 + fy10/12 - fy11/12);
    fint_y = sx*(f01/2 - f00/2 - f10/2 + f11/2 - fx00/12 + fx01/12 + fx10/12 - fx11/12);
end