function [ Ih ] = gauss_legendre_simple_quadrature( fun, a, b, n )
% GAUSS_LEGENDRE_SIMPLE_QUADRATURE approximate the integral of a function in
% the interval [a,b] by means of the simple Gauss-Legendre quadrature formula
%  [ Ih ] = gauss_legendre_simple_quadrature( fun, a, b, n )
%  Inputs: fun = function handle, 
%          a,b = extrema of the interval [a,b]
%          n + 1 = number of quadrature nodes and weights
%  Output: Ih = approximate value of the integral
%

% reference nodes and weights
switch n
    
    case 0
        y_ref = 0;
        alpha_ref = 2;
    
    case 1
        y_ref = [ -1/sqrt(3), 1/sqrt(3) ]; 
        alpha_ref = [ 1, 1 ];        
        
    case 2
        y_ref = [ -sqrt(15)/5, 0, sqrt(15)/5 ]; 
        alpha_ref = [ 5/9, 8/9, 5/9 ];                
    
    otherwise
        error('n must be 0, 1, or 2');
        
end

% nodes and weigths rescaled in the interval [a,b]
y_rescaled = ( a + b ) / 2 + ( b - a ) / 2 * y_ref;
alpha_rescaled = ( b - a ) / 2 * alpha_ref;

% Integral
Ih = sum( alpha_rescaled .* fun( y_rescaled ) ); 

return