%% Gaussian Quadrature for Double Integrals

% Compute tensor-product nodes and weights.
n = 7;
[x, c] = gaussquad(n);
y = x';
w = c * c';

%% Integrate exp(x + y) on [0, 1] by [0, 2]

x1 = (x + 1) / 2;
y1 = y + 1;
w1 = w / 2;
J1 = sum(exp(x1 + y1) .* w1, 'all');

% Compare with Simpson
J1s = simpsondouble(@(x, y) exp(x + y), ...
    @(x) 0, @(x) 2, 0, 1, 16, 16);
fprintf('Rectangle: Gaussian = %.12g, Simpson = %.12g\n', J1, J1s)

%% Integrate exp(x + y) for 0 <= x <= 1 and x <= y <= 1 + x / 2

x2 = (x + 1) / 2;
y2 = x2 + (y + 1) / 2 .* (1 - x2 / 2);
w2 = w .* (1 - x2 / 2) / 4;
J2 = sum(exp(x2 + y2) .* w2, 'all');

% Compare with Simpson
J2s = simpsondouble(@(x, y) exp(x + y), ...
    @(x) x, @(x) 1 + x / 2, 0, 1, 16, 16);
fprintf('Curved domain: Gaussian = %.12g, Simpson = %.12g\n', J2, J2s)

% Plot domain and Gaussian integration nodes
plot(x2, y2, '.k')
line([0, 1, 1, 0, 0], [0, 1, 1.5, 1, 0])
axis equal
