function z=gquad(fname,a,b,n,alg,varargin)
% gquad - Numerical integration of f(x) over [a;b] using
%         Gaussian quadrature.
%
% out=gquad(fname,a,b,n,alg,x1,x2,...,xn)
%
% Input:  fname - name of the m-file implementing f(x)
%         a - lower bound
%         b - upper bound
%         n - number of nodes
%         alg - method
%               'c':  Gauss-Chebyshev
%               'l':  Gauss-Legendre
%               'l':  Simpson's rule
%               't':  Trapezoid rule
%               'gl': Gauss-Lobatto
%         x1,x2,...,xn - arguments passed to fname
%
% Output: Integral of f(x)
%
% (c) 2012 by Marco Maffezzoli. Ver. 1.0.2.
%

n=max(n,3);
switch alg
    case {'c','l'}
        [m,w]=qnodes(n,alg);
        w=(b-a)/2*w;
        x=(m+1)*((b-a)/2)+a;
    case {'s','t'}
        h=(b-a)/(floor(n)-1);
        x=(a:h:b).';
    case {'gl'}
        n=max(3*round((n-1)/3),3)+1;
        h=(b-a)/(n-1);
        d=(3/realsqrt(5)-1)*h/2;
        x=(a:h:b).';
        x(2:3:n-2)=x(2:3:n-2)-d;
        x(3:3:n-1)=x(3:3:n-1)+d;
    otherwise
        error('Choose an algorithm, please!');
end
f=feval(fname,x,varargin{:});
f=f(:);
switch alg
    case 'c' % Gauss-Chebyshev
        z=w'*(f.*realsqrt(1-m.^2));
    case 'l' % Gauss-Legendre
        z=w'*f;
    case 's' % Simpson's rule
        f1=2*sum(f(3:2:n-2),1);
        f2=4*sum(f(2:2:n-1),1);
        z=h/3*(f(1)+f1+f2+f(n));
    case 't' % Trapezoid rule
        z=h*(sum(f,1)-(f(1)+f(n))/2);
    case 'gl' % Gauss-Lobatto rule
        f1=2*sum(f(4:3:n-3),1);
        f2=5*sum(f([3:3:n-1,2:3:n-2]),1);
        z=h/4*(f(1)+f1+f2+f(n));
end

end