function [i1,i2]=gacheb(func,n)
% Implements Gauss-Chebyshev integration.
%
% Example call: [i1,i2] = gacheb(func,n)
% 
% Integrates user defined function func (times a square root function) from -1 to +1 using n divisions.
% 
% i1 = \int f(x) / sqrt( 1 - x^2 ) dx 
% 
% i2 = \int sqrt( 1 - x^2 ) f(x) dx 
%
% See EPage 74 in the book.
%
% Written by:
% --
% John L. Weatherwax                2009-04-21
%
% email: wax@alum.mit.edu
%
% Please send comments and especially bug reports to the
% above email address.
%
%-----

ks = 1:n;

% For i1: 
% 
x_k   = cos( ( 2 * ks - 1 ) * pi / ( 2 * n ) );
f_x_k = feval(func,x_k);
i1 = ( pi / n ) * sum( f_x_k );

% For i2: 
%
x_k   = cos( ks * pi / ( n + 1 ) );
f_x_k = feval(func,x_k);
s2 = sin( ks * pi / ( n + 1 ) ).^2;
i2 = ( pi / ( n + 1 ) ) * sum( s2 .* f_x_k );