function [x,w,q]=qgausc(a,b,n,f,par);
%QGAUSC  Numerical evaluation of integrals (quadrature).
%        [x,w]=qgausc(a,b,n) computes the abscissas and weights 
%        of the Gauss-Chebyshev n-point quadrature formula with
%        integration limits (a,b).  
%        [x,w,q]=qgausc(a,b,n,f,par) computes the integral
%                  b
%          q =   \int f(x) dx   ~=  \sum w[i] f(x[i],par)* sqrt(1-y[i]^2)
%                  a                  i   
%
%        where y[i] = (2*x[i]-a-b)/(b-a) in addition to the abscissas 
%        and weights.  'f' is a string containing the name of the function.
%


% recursively compute coefficients of Chebyshev polynomials
p        = zeros(n+1);
p(1,n+1) = 1; 
p(2,n)   = 1;
for i=2:n;  
  p(i+1,:) = 2*[p(i,2:n+1),0]-p(i-1,1:n+1) ;
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 \ [asin(1)-asin(-1);zeros(n-1,1)];
w        = 0.5*(b-a)*w;

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

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


