% TEST_GLOBATTO_RGT0 Moment test for generalized Gauss-Lobatto formula
%
f1='%5.0f %12.2e %12.2e %8.0f %8.0f\n';
disp('         maximum errors in moments for r up to 10')
disp('    N     maxerrleg   maxerrcheb   rmaxleg  rmaxcheb')
endl=-1; endr=1;
for N=3:7
%
% Legendre measure
%
  maxerrleg=0;
  for r=1:10
    mom=zeros(2*N+2*r,1); errmom=zeros(2*N+2*r,1);
    ab=r_jacobi(N+2*r);
    xw=globatto(N,ab,r,endl,endr);
    for k=1:2:2*N+2*r
      mom(k)=2/k;
    end
    for k=1:2*N+2*r
      s=0;
      for rho=1:r
        s=s+xw(rho,2)*fpow(xw(rho,1),k-1,rho-1);
        s=s+xw(r+rho,2)*fpow(xw(r+rho,1),k-1,rho-1);
      end
      s=s+sum(xw(2*r+1:2*r+N,2).*xw(2*r+1:2*r+N,1).^(k-1));
      errmom(k)=abs(mom(k)-s);
    end
    merr=max(errmom);
    if merr>maxerrleg
       maxerrleg=merr;   
       rmaxleg=r;
    end
  end
%
% Chebyshev measure
%
  maxerrcheb=0;
  for r=1:10
    mom=zeros(2*N+2*r,1); errmom=zeros(2*N+2*r,1);
    ab=r_jacobi(N+2*r,-.5);
    xw=globatto(N,ab,r,endl,endr);
    mom(1)=pi;
    for k=3:2:2*N+2*r
      mom(k)=(pi/(k-1))*factorial(k-2)/(2^(k-3)...
        *(factorial((k-3)/2))^2);
    end
    for k=1:2*N+2*r
      s=0;
      for rho=1:r
        s=s+xw(rho,2)*fpow(xw(rho,1),k-1,rho-1);
        s=s+xw(r+rho,2)*fpow(xw(r+rho,1),k-1,rho-1);
      end
      s=s+sum(xw(2*r+1:2*r+N,2).*xw(2*r+1:2*r+N,1).^(k-1));
      errmom(k)=abs(mom(k)-s);
    end
    merr=max(errmom);
    if merr>maxerrcheb
       maxerrcheb=merr;   
       rmaxcheb=r;
    end
  end
  fprintf(f1,N,maxerrleg,maxerrcheb,rmaxleg,rmaxcheb)
end
