% TABLE2_25 Computations for Table 2.25.
%
f1='%6.3f %4.1f %8.1e %8.1e %8.1e %8.1e %8.1e\n';
f2='%20.1e %8.1e %8.1e %8.1e %8.1e\n';
Fex=[.2801911343790928e2 .2988804854276546e2 -.2324134092608720e2...
.3122746491469022e9;.3104991751202353e3 .3127405406398423e3...
.2624298811059665e3 .3119650155664194e9;.3137904492745048e4...
.3140187710458022e4 .3090172669532817e4 .3119647409171416e9];
N=80; ab0=r_hermite(N); 
fprintf('\n')
disp('        Gauss quadrature for the integral F(x,y)')
fprintf('\n')
disp('   y     x     n=5     n=10     n=20     n=40     n=80')
ky=0;
for y=[.1 .01 .001]
  ky=ky+1;
  rho0=complex(0,-pi*erfcx(y));
  ab=chri6(N,ab0,y,eps0,nu0,numax,rho0,2);
  kx=0;
  for x=[1 2 5 10]
    kx=kx+1;
    kn=0;
    err=zeros(5); err0=zeros(5);
    for n=[5 10 20 40 80]
      kn=kn+1;
      xw=gauss(n,ab); xw0=gauss(n,ab0);
      F=sum(xw(:,2).*exp(-x*xw(:,1)).*cos(xw(:,1)));
      F0=sum(xw0(:,2).*exp(-x*xw0(:,1)).*cos(xw0(:,1))...
         ./((xw0(:,1)).^2+y^2));
      err(kn)=abs((F-Fex(ky,kx))/Fex(ky,kx));
      err0(kn)=abs((F0-Fex(ky,kx))/Fex(ky,kx));
    end
    fprintf(f1,y,x,err(1),err(2),err(3),err(4),err(5))
    fprintf(f2,err0(1),err0(2),err0(3),err0(4),err0(5))
  end
end
