% TABLE2_24  Computations for Table2.24.
%
f1='%18.4f %6.0f %5.0f %5.0f %5.0f %12.3e %11.3e\n';
fprintf('\n')
disp('           Modification by a quadratic divisor using gchri2 and chri5')
fprintf('\n')
N=40; numax=450; eps0=100*eps;
tic
%for a=-.8:.4:.8
%  for b=a:.4:.81
for a=0
  for b=0
    disp('              x       nu0    nu  nus0  nus0      erra        errb')
    rr=[1.05 1.1625 1.275 1.3875 1.5];
    for ir=1:5
      r=rr(ir);
      th=linspace(0,pi,21);
      nu0av=0; nus0av=0;
      nuav=0; nusav=0;
      errav=zeros(1,2);
      for iz=2:20
        z(iz)=.5*(r*exp(complex(0,th(iz)))+...
        exp(-complex(0,th(iz)))/r);
        x=real(z(iz)); y=imag(z(iz));
        nu0=nu0jac(2*N-1,z(iz),eps0);
        nus0=nu0jac(N,z(iz),eps0);
        ab0=r_jacobi(numax,a,b);
        [abgc,rho,nu]=gchri2(N,ab0,x,y,eps0,nu0,numax);
        [abcs,nus]=chri5(N,ab0,z(iz),eps0,nus0,numax,rho(1),1);
        erra=abs(abgc(:,1)-abcs(:,1));
        errb=abs((abgc(:,2)-abcs(:,2))./abcs(:,2));
        err=[erra errb];
        errm=max(err);
        nu0av=nu0av+nu0; nuav=nuav+nu;
        nus0av=nus0av+nus0; nusav=nusav+nus;
        errav=errav+errm;
      end
      nu0av=nu0av/19; nuav=nuav/19;
      nus0av=nus0av/19; nusav=nusav/19;
      errav=errav/19;
      fprintf(f1,r,nu0av,nuav,nus0av,nusav,errav)
    end
    fprintf('\n')
  end
end
toc
