function ch07_nested_simpson_2d
%==========================================================================
% ch07_nested_simpson_2d.m
%
% Two-dimensional integration using nested Simpson integration.
%
% The program evaluates
%
%          1   sin(x)
%          /      /
%         |      |       x^2
%     I = |      |   --------- dy dx .
%         |      |    y^2 + 2
%         /      /
%          0      0
%
% The two-dimensional problem is reduced to repeated one-dimensional
% integrations. The same Simpson refinement routine is used for the
% inner and outer integrations.
%
% The number of intervals is doubled until
%
%       |I_(2n)-I_n|/15 < tolerance.
%
% This is automatic global refinement, not local adaptive subdivision.
%
% Based on code written by Alexander Godunov, October 2009.
% 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.
%==========================================================================

a = 0.0;
b = 1.0;

eps_outer = 1.0e-7;
eps_inner = 1.0e-8;

reference = 1.0344649764293148e-1;

% === Plot the integrand and integration domain ==========================
% The upper boundary of the inner integral is y = sin(x).  Values outside
% the physical integration domain are masked so the surface shows both
% the integrand and the variable boundary.
xplot = linspace(a,b,140);
yplot = linspace(0.0,sin(b),140);
[Xplot,Yplot] = meshgrid(xplot,yplot);

Zplot = Xplot.^2./(Yplot.^2 + 2.0);
Zplot(Yplot > sin(Xplot)) = NaN;

figure('Name','Nested two-dimensional integrand')
surf(Xplot,Yplot,Zplot,'EdgeColor','none')
xlabel('x')
ylabel('y')
zlabel('f(x,y)')
title('Integrand over 0 <= y <= sin(x)')
view(45,30)
grid on
% === End plot ===========================================================

nfun_2d = 0;

[integral,error_estimate,n_used] = ...
    simpson_refine(@inner_integral,a,b,eps_outer);

absolute_error = abs(integral-reference);

fprintf('Two-dimensional nested Simpson integration\n');
fprintf('Integral of x^2/(y^2+2), 0 <= x <= 1,\n');
fprintf('with 0 <= y <= sin(x)\n');
fprintf('Reference value      = %15.7e\n', reference);
fprintf('Calculated value     = %15.7e\n', integral);
fprintf('Estimated outer error= %12.4e\n', error_estimate);
fprintf('Absolute error       = %12.4e\n', absolute_error);
fprintf('Outer intervals      = %10d\n', n_used);
fprintf('2D function calls    = %10d\n', nfun_2d);


    function value = integrand(x,y)
        %------------------------------------------------------------------
        % Original two-dimensional integrand.
        %------------------------------------------------------------------
        nfun_2d = nfun_2d + 1;
        value = x*x/(y*y+2.0);
    end


    function value = y_lower(x)
        %#ok<INUSD>
        value = 0.0;
    end


    function value = y_upper(x)
        value = sin(x);
    end


    function value = inner_integral(x)
        %------------------------------------------------------------------
        % Evaluate the inner integral F(x).
        %------------------------------------------------------------------
        c = y_lower(x);
        d = y_upper(x);

        [value,~,~] = simpson_refine(@(y) integrand(x,y), ...
                                     c,d,eps_inner);
    end

end


function [integral,error_estimate,n_used] = simpson_refine(fun,a,b,eps)
%--------------------------------------------------------------------------
% One-dimensional Simpson integration with automatic global refinement.
%
% Composite Simpson estimates are generated for
%
%       n = 2, 4, 8, 16, ...
%
% and refinement stops when
%
%       |I_(2n)-I_n|/15 <= eps.
%--------------------------------------------------------------------------

nmax = 1048576;

if eps <= 0.0
    error('simpson_refine: tolerance must be positive.');
end

if abs(b-a) <= realmin
    integral = 0.0;
    error_estimate = 0.0;
    n_used = 0;
    return;
end

h = (b-a)/2.0;
sn = h*(fun(a) + 4.0*fun(a+h) + fun(b))/3.0;

n = 4;

while n <= nmax

    h = (b-a)/n;

    s2n = fun(a) + fun(b);

    for i = 1:n-1
        x = a+i*h;

        if mod(i,2) == 0
            s2n = s2n + 2.0*fun(x);
        else
            s2n = s2n + 4.0*fun(x);
        end
    end

    s2n = h*s2n/3.0;

    error_estimate = abs(s2n-sn)/15.0;

    if error_estimate <= eps
        integral = s2n;
        n_used = n;
        return;
    end

    sn = s2n;
    n = 2*n;
end

integral = sn;
n_used = n/2;

warning('simpson_refine: maximum number of intervals reached.');

end