function [x,w,q]=qgausl(a,b,n,f,par);
%QGAUSL  Numerical evaluation of integrals (quadrature).
%        [x,w]=qgausl(a,b,n) computes the abscissas and weights 
%        of the Gauss-Legendre n-point quadrature formula with
%        integration limits (a,b).  
%        [x,w,q]=qgausl(a,b,n,f,par) computes the integral
%                  b
%          q =   \int   f(x)  ~=  \sum w[i] f(x[i],par)
%                  a                i   
%        in addition to the abscissas and weights.  'f' is a 
%        string containing the name of the function.
%

%        Ellen McGrattan, 12-5-91
%        Revised, 12-9-95, ERM

% recursively compute coefficients of Legendre polynomials
p        = zeros(n+1);
p(1,n+1) = 1; 
p(2,n)   = 1;
for i=2:n;  
  p(i+1,:) = ( (2*i-1)*[p(i,2:n+1),0]-(i-1)*p(i-1,1:n+1) )/i;
end;

% the roots of the polynomials are the abscissas
x        = sort(roots(p(n+1,:)));
A        = p(1:n,2:n+1)*( (ones(n,1)*x').^([n-1:-1:0]'*ones(1,n)) );

% the weights are derived by noting that the quadrature formula must hold
% for the first n-1 Legendre polynomials -- and with these weights the next
% n hold exactly
w        = A \ [b-a;zeros(n-1,1)];

% map abscissas from [-1,1] to [a,b]
x        = 0.5 * (a *(1-x) + b * (1+x));

% compute the integral if all input and output arguments are given
if nargin>3 & nargout>2;    
  q      = w'*feval(f,x,par);
end;

