function ch07_gauss_legendre_2d
%==========================================================================
% ch07_gauss_legendre_2d.m
%
% Two-dimensional integration using a Gauss-Legendre product rule.
%
% The program evaluates
%
%          1   1
%          /   /
%         |   |        1
%     I = |   |   -------------  dy dx
%         |   |   1 + x^2 + y^2
%         /   /
%          0   0
%
% An n-point Gauss-Legendre rule is used in each direction. The
% tensor-product rule therefore requires n^2 function evaluations.
%
% Based on the Fortran implementation by Alexander Godunov.
% MATLAB version prepared for the companion website, 2026.
%
% A plot of the integrand (or integration domain) is shown before the
% numerical calculation to make the character of the problem visible.
%==========================================================================

ax = 0.0;
bx = 1.0;
ay = 0.0;
by = 1.0;

reference = 6.3951035187031102e-01;
nvalues = [2, 4, 8, 16, 32];

% === Plot the integrand =================================================
xplot = linspace(ax,bx,120);
yplot = linspace(ay,by,120);
[Xplot,Yplot] = meshgrid(xplot,yplot);
Zplot = 1.0./(1.0 + Xplot.^2 + Yplot.^2);

figure('Name','Two-dimensional integrand')
surf(Xplot,Yplot,Zplot,'EdgeColor','none')
xlabel('x')
ylabel('y')
zlabel('f(x,y)')
title('Integrand: 1/(1+x^2+y^2)')
view(45,30)
grid on
% === End plot ===========================================================

fprintf('Two-dimensional Gauss-Legendre integration\n');
fprintf('Integral of 1/(1+x^2+y^2) over [0,1] x [0,1]\n');
fprintf('Reference value      = %15.7e\n\n', reference);

fprintf('        n     Function calls          Integral       Absolute error\n');
fprintf('-------------------------------------------------------------------\n');

for n = nvalues
    [integral, nfun] = gauss_legendre_2d(@f, ax, bx, ay, by, n);
    error_value = abs(integral-reference);

    fprintf('%9d%19d%20.7e%20.7e\n', ...
            n, nfun, integral, error_value);
end

end


function z = f(x,y)
%--------------------------------------------------------------------------
% Function to be integrated over the rectangular domain.
%--------------------------------------------------------------------------
z = 1.0/(1.0 + x*x + y*y);
end


function [integral, nfun] = gauss_legendre_2d(fun, ax, bx, ay, by, n)
%--------------------------------------------------------------------------
% Two-dimensional Gauss-Legendre tensor-product rule.
%--------------------------------------------------------------------------

if n < 1
    error('gauss_legendre_2d: n must be positive.');
end

[xnode,xweight] = gauss_legendre_nodes_weights(n,ax,bx);
[ynode,yweight] = gauss_legendre_nodes_weights(n,ay,by);

integral = 0.0;
nfun = 0;

for i = 1:n
    for j = 1:n
        integral = integral + ...
            xweight(i)*yweight(j)*fun(xnode(i),ynode(j));
        nfun = nfun + 1;
    end
end

end


function [node,weight] = gauss_legendre_nodes_weights(n,a,b)
%--------------------------------------------------------------------------
% Generate n Gauss-Legendre nodes and weights on [a,b].
%--------------------------------------------------------------------------

tol = 1.0e-14;
max_iter = 100;

if n < 1
    error('gauss_legendre_nodes_weights: n must be positive.');
end

node = zeros(1,n);
weight = zeros(1,n);

midpoint = (a+b)/2.0;
halfwidth = (b-a)/2.0;
m = floor((n+1)/2);

for i = 1:m

    z = cos(pi*(i-0.25)/(n+0.5));
    converged = false;

    for iter = 1:max_iter
        [pn,pnm1] = legendre_pair(n,z);

        if n == 1
            dp = 1.0;
        else
            dp = n*(z*pn-pnm1)/(z*z-1.0);
        end

        zold = z;
        z = zold - pn/dp;

        if abs(z-zold) <= tol
            converged = true;
            break;
        end
    end

    if ~converged
        error('Gauss-Legendre root iteration did not converge.');
    end

    [pn,pnm1] = legendre_pair(n,z);

    if n == 1
        dp = 1.0;
    else
        dp = n*(z*pn-pnm1)/(z*z-1.0);
    end

    w = 2.0/((1.0-z*z)*dp*dp);

    left = i;
    right = n+1-i;

    node(left) = midpoint-halfwidth*z;
    node(right) = midpoint+halfwidth*z;

    weight(left) = halfwidth*w;
    weight(right) = halfwidth*w;
end

end


function [pn,pnm1] = legendre_pair(n,z)
if n == 1
    pn = z;
    pnm1 = 1.0;
    return;
end

p0 = 1.0;
p1 = z;

for j = 2:n
    p2 = ((2.0*j-1.0)*z*p1 - (j-1.0)*p0)/j;
    p0 = p1;
    p1 = p2;
end

pn = p1;
pnm1 = p0;
end