% TEST_GLOBATTO_REQ2 Test of globatto.m against known Gauss- 
%    Lobatto formulae with double endpoints 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; endr=1;
for N=[10 20 40 80 160]
%
% Legendre measure
%
  ab=r_jacobi(N+2*r);
  xw=globatto(N,ab,r,endl,endr);
  abex=r_jacobi(N,2);
  xwex=gauss(N,abex);
  xwex(:,2)=xwex(:,2)./((1-xwex(:,1).^2).^2);
  errinleg=max(abs((xw(5:N+4,2)-xwex(:,2))./xwex(:,2)));
  den=(N+1)*(N+2)*(N+3)*(N+4);
  w1=8*(2*N^2+10*N+9)/(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+2*r,-.5); 
  xw=globatto(N,ab,r,endl,endr);
  abex=r_jacobi(N,1.5);
  xwex=gauss(N,abex);
  xwex(:,2)=xwex(:,2)./((1-xwex(:,1).^2).^2);
  errincheb=max(abs((xw(5:N+4,2)-xwex(:,2))./xwex(:,2)));
  den=(N+1)*(N+2)*(N+3);
  w1=3*pi*(3*N^2+12*N+10)/(10*den);
  w2=3*pi/(4*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
