% gauss.m -- Compute int(f(x),x=-1..1) using Gaussian quadrature
%
%     s=gauss(k,f)
%
% Written Fall 2005 by Eric Olson for MATH/CS 466/666
%
% The table of roots and weights were computed using 32 digits of 
% precision by the Maple script legendre.mpl and then rounded to 
% double precision for use here as Matlab input.

function s=gauss(k,f)
    roots=[ 0 -0.5773502691896258 0.5773502691896258 -0.7745966692414834
        0 0.7745966692414834 -0.8611363115940526 -0.3399810435848563
        0.3399810435848563 0.8611363115940526];
    weights=[ 2 1 1 0.5555555555555556 0.8888888888888889 0.5555555555555556
        0.3478548451374539 0.6521451548625461 0.6521451548625461
        0.3478548451374539];
    n=1+k*(k-1)/2;
    if n<1 || n+k-1>length(roots)
        disp(sprintf('Gauss rule for %d points is not here!',k));
    end
    s=0.0;
    for i=n:n+k-1
        s=s+weights(i)*feval(f,roots(i));
    end
end
