%============================================================
% Two electron integrals
%   Computes the two electron integral.
%------------------------------------------------------------
% Note the use of two fudge-factors:
%   
%------------------------------------------------------------
% Parameters:
%   f1(r), f2(r), f3(r), f4(r): function pointers accepting 
%                               r = [x y z]. Interpreted as:
%                               <f1 f2|1/r{12}|f3 f4>
%   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 integral = two_elec_integral(f1,f2,f3,f4,gp)

% Fudge factors:
eps = 1; % Prevents division by zeros that occur when using
               % identical grids for r1 and r2
               
div = 1;     % Makes the integration grid coarser by the factor
               % `div'. This is necessary because we don't want to
               % wait three minutes for a single two-electron integral.

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

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

dV2 = (dx*dy*dz)^2; % Note that the two electron integral is 
                    % a six-dimensional volume integral

% A small fudge parameter to prevent divide-by-zeros
eps = exp(-3);

integral = 0;
%------------------------------------------------------------
% BEGIN: "Outer" integral over r1
%------------------------------------------------------------
for x1 = linspace(x_lower, x_upper, x_steps)
    for y1 = linspace(y_lower, y_upper, y_steps)
        for z1 = linspace(z_lower, z_upper, z_steps)
            
        r1 = [x1 y1 z1];
        %------------------------------------------------------------
        % BEGIN: "Inner" integral over r2
        %------------------------------------------------------------
        for x2 = linspace(x_lower+eps, x_upper+eps, x_steps)
            for y2 = linspace(y_lower+eps, y_upper+eps, y_steps)
                for z2 = linspace(z_lower+eps, z_upper+eps, z_steps)
                    r2 = [x2 y2 z2];
                    r12inv   = 1/norm(r1-r2);
                    integral = integral + ...
                               dV2 * conj( f1(r1)*f2(r2) ) * ...
                               r12inv * f3(r1) * f4(r2);
                end
            end
        end
        
        %------------------------------------------------------------
        % END: "Inner" integral over r2
        %------------------------------------------------------------
        end
    end
end

