% TEST_GRADAU_REQ2 Test of gradau.m against known Gauss-Radau formulae 
%    with double endpoint for the Legendre and Chebyshev measures
%
f1='%5.0f %10.2e %10.2e %10.2e %10.2e\n';
disp('           maximum errors in weights for r=2')
disp('   N    errinleg   errexleg   errincheb  errexcheb')
r=2; endl=-1; 
for N=[10 20 40 80 160]
%
% Legendre measure
%
  ab=r_jacobi(N+r);
  xw=gradau(N,ab,r,endl);
  abex=r_jacobi(N,0,2);
  xwex=gauss(N,abex);
  xwex(:,2)=xwex(:,2)./((1+xwex(:,1)).^2);
  errinleg=max(abs((xw(3:N+2,2)-xwex(:,2))./xwex(:,2)));
  den=(N+1)^2*(N+2)^2;
  w1=8*(2*N^2+6*N+3)/(3*den);
  w2=8/den;
  err1=abs((xw(1,2)-w1)/w1);
  err2=abs((xw(2,2)-w2)/w2);
  errexleg=max([err1 err2]);
%
% Chebyshev measure
%
  ab=r_jacobi(N+r,-.5); 
  xw=gradau(N,ab,r,endl);
  abex=r_jacobi(N,-.5,1.5);
  xwex=gauss(N,abex);
  xwex(:,2)=xwex(:,2)./((1+xwex(:,1)).^2);
  errincheb=max(abs((xw(3:N+2,2)-xwex(:,2))./xwex(:,2)));
  den=(N+1)*(2*N+1)*(2*N+3);
  w1=3*pi*(6*N^2+12*N+5)/(5*den);
  w2=3*pi/den;
  err1=abs((xw(1,2)-w1)/w1);
  err2=abs((xw(2,2)-w2)/w2);
  errexcheb=max([err1 err2]);
  fprintf(f1,N,errinleg,errexleg,errincheb,errexcheb)
end
