%------------------------------------------------------------
% Inner product
%   Computes the inner product of two functions over the
%   specified grid.
%------------------------------------------------------------
% Parameters:
%   f1(r), f2(r): function pointers accepting r = [x y z]
%   gp: "grid parameters"
%       integration grid in the following format
%            [x_lower x_upper x_steps;...
%             y_lower y_upper y_steps;...
%             z_lower z_upper z_steps];
%------------------------------------------------------------
function ip = inner_product(f1, f2, gp)

x_lower = gp(1,1); x_upper = gp(1,2); x_steps = gp(1,3);
y_lower = gp(2,1); y_upper = gp(2,2); y_steps = gp(2,3);
z_lower = gp(3,1); z_upper = gp(3,2); z_steps = gp(3,3);

dx = (x_upper - x_lower)/(x_steps-1);
dy = (y_upper - y_lower)/(y_steps-1);
dz = (z_upper - z_lower)/(z_steps-1);
dV = dx*dy*dz;

ip = 0;
for x = linspace(x_lower,x_upper,x_steps)
    for y = linspace(y_lower,y_upper,y_steps)
        for z = linspace(z_lower,z_upper,z_steps)
            ip = ip + dV * conj(f1([x y z])) * f2([x y z]);
        end
    end
end