% EXAMPLE2_38 Computations for Example 2.38.
%
eps0=1e-3; exact=2+pi;
%
% Fejer
%
fprintf('\n')
fprintf('        Fejer quadrature')
fprintf('\n')
N=1000; n=685; int=0;
while abs(int-exact)>eps0
  n=n+1; if n>N, disp('n exceeds N'), break, end
  xw=fejer(n);
  int=sum(xw(:,2).*(1+1./sqrt(1-xw(:,1).^2)))
end
n
%
% Gauss-Legendre
%
fprintf('\n')
fprintf('   Gauss-Legendre quadrature')
fprintf('\n')
N=700; n=690; int=0;
ab=r_jacobi(N);
while abs(int-exact)>eps0
  n=n+10; if n>N, disp('n exceeds 700'), break, end
  xw=gauss(n,ab);
  int=sum(xw(:,2).*(1+1./sqrt(1-xw(:,1).^2))) 
end
n
%
% Gauss Chebyshev
%
fprintf('\n')
fprintf('   Gauss-Chebyshev quadrature')
fprintf('\n')
N=700; n=25; int=0;
ab=r_jacobi(N,-1/2);
while abs(int-exact)>eps0
  n=n+1; if n>N, disp('n exceeds 700'), break, end
  xw=gauss(n,ab);
  int=sum(xw(:,2).*(1+sqrt(1-xw(:,1).^2)))
end
n
fprintf('\n')
eps0=1e-5;
n=280; int=0;
while abs(int-exact)>eps0
  n=n+1; if n>N, error('n exceeds 700'), end
  xw=gauss(n,ab);
  int=sum(xw(:,2).*(1+sqrt(1-xw(:,1).^2)))
end
n
