% Gauss_Legendre quadrature
clear; clc;

a = 0; b = pi;

f = @(x) exp(x).*cos(x);

for n = 1:6
    [Gauss_point_reference_1D, Gauss_coefficient_reference_1D] = generate_Gauss_reference_1D(n);
    [Gauss_point_local_1D, Gauss_coefficient_local_1D] = generate_Gauss_local_1D(a,b,Gauss_point_reference_1D, Gauss_coefficient_reference_1D);

    I =0.0;
    for k = 1:length(Gauss_point_local_1D)
        I = I + Gauss_coefficient_local_1D(k)*f(Gauss_point_local_1D(k));
    end
    fprintf("number of Gauss points=%d, I=%0.8f\n", length(Gauss_point_local_1D), I);
end


function [Gauss_point_local_1D, Gauss_coefficient_local_1D] = generate_Gauss_local_1D(a,b,Gauss_point_reference_1D, Gauss_coefficient_reference_1D)
Gauss_point_local_1D = (b-a)/2 * Gauss_point_reference_1D + (b+a)/2;
Gauss_coefficient_local_1D = (b-a)/2 * Gauss_coefficient_reference_1D;
end

function [Gauss_point_reference_1D, Gauss_coefficient_reference_1D] = generate_Gauss_reference_1D(n)
syms x;
Poly = diff((x^2-1)^(n+1),x,n+1) / factorial(n+1) / 2^(n+1);
Gauss_point_reference_1D = solve(Poly);
dPoly = diff(Poly, x, 1);
Gauss_coefficient_reference_1D = zeros(size(Gauss_point_reference_1D));
for k = 1:length(Gauss_coefficient_reference_1D)
    Gauss_coefficient_reference_1D(k) = 2 / ((1-Gauss_point_reference_1D(k).^2) .* subs(dPoly,x,Gauss_point_reference_1D(k)).^2);
end
end