% EXAMPLE2_29 Numerical illustration for Example 2.29.
%
% See also  W. Gautschi, ``Algorithm 726: ORTHPOL -- a package of routines
% for generating orthogonal polynomials and Gauss-type quadrature rules'',
% ACM Trans. Math. Software 20 (1994), 21-62, Example 3.1.
%
f1='%9.3f %4.0f %20.15f %9.2e\n';
f2='%14.0f %20.15f %9.2e\n';
fprintf('\n')
disp('   Recursion coefficients for an elliptic weight')
disp('   function computed via Chebyshev moments')
fprintf('\n')
disp('     k2      n         beta            magb')
load -ascii coeffell;
bexact=coeffell;
N=80; abm=r_jacobi(2*N-1,-1/2);
for k2=[.1 .5 .9 .999]
  mom=mm_ell(N,k2);
  [ab,normsq]=chebyshev(N,mom,abm);
  if k2==.1
    magb=abs(ab(:,2)-bexact(1:N))/eps;
    n=[1 2 6 12];
    for k=1:4
      if k==1
        fprintf(f1,k2,n(k)-1,ab(n(k),2),magb(n(k)))
      else
        fprintf(f2,n(k)-1,ab(n(k),2),magb(n(k)))
      end
    end
    fprintf('\n')
  elseif k2==.5
    magb=abs(ab(:,2)-bexact(N+1:2*N))/eps;
    n=[1 2 9 21];
    for k=1:4
      if k==1
        fprintf(f1,k2,n(k)-1,ab(n(k),2),magb(n(k)))
      else
        fprintf(f2,n(k)-1,ab(n(k),2),magb(n(k)))
      end
    end
    fprintf('\n')
  elseif k2==.9
    magb=abs(ab(:,2)-bexact(2*N+1:3*N))/eps;
    n=[1 2 20 44];
    for k=1:4
      if k==1
        fprintf(f1,k2,n(k)-1,ab(n(k),2),magb(n(k)))
      else
        fprintf(f2,n(k)-1,ab(n(k),2),magb(n(k)))
      end
    end
    fprintf('\n')
  else
    magb=abs(ab(:,2)-bexact(3*N+1:4*N))/eps;
    n=[1 2 20 44 80];
    for k=1:5
      if k==1
        fprintf(f1,k2,n(k)-1,ab(n(k),2),magb(n(k)))
      else
        fprintf(f2,n(k)-1,ab(n(k),2),magb(n(k)))
      end
    end
  end
end