% QUADBESS A special quadrature routine used in Table2_19.
%
%    This uses approximation #26 for R(t) and approximation #111 for S(t)
%    from A.E. Russon and J.M. Blair, "Rational function minimax approximations
%    for the modified Bessel functions K_0(x) and K_1(x)", AECL-3461, Atomic
%    Energy of Canada Limited, Chalk River, Ontario, October 1969.
% 
function xw=quadbess(N,i)
global a abjac abjaclog ablag
Rnum=[1.1826895958076190115e3 2.8266004990762514731e3...
  2.108241953516787094e2 4.6678516537403681...
  3.52613614881237e-2];
Rden=[1.0201622820945108630e4 -1.68009578767071511e2 1];
Snum=[1.64399227014181315109e2 2.187632726229839382872e3...
  1.0504822074208022979599e4 2.3001392658611199096172e4...
  2.3945021887112558158991e4 1.1252648779364127795637e4...
  2.031037078778050605207e3 9.2077106373360784288e1];
Sden=[1.31171605042540629856e2 1.761874827596169863623e3...
  8.592646696965562512522e3 1.9312262387540974113366e4...
  2.1029558009481191713216e4 1.0710652023789478437783e4...
  2.256560775942557871261e3 1.46075579596262287743e2 1];
if i==1
  xw=gauss(N,abjac);
  x=xw(:,1); num=zeros(size(x)); den=zeros(size(x));
  for k=1:5
    num=num+Rnum(k)*(x.^(2*k-2));
  end
  for k=1:3
    den=den+Rden(k)*(x.^(2*k-2));
  end
  R=num./den;
  xw(:,2)=R.*xw(:,2);
elseif i==2
  xw=gauss(N,abjaclog);
  x=xw(:,1);
  xw(:,2)=bessel(0,j*x).*xw(:,2);
else
  xw=gauss(N,ablag);
  xw(:,1)=1+xw(:,1);
  x=xw(:,1); num=zeros(size(x)); den=zeros(size(x));
  for k=1:8
    num=num+Snum(k)*(x.^(-k+1));
    den=den+Sden(k)*(x.^(-k+1));
  end
  den=den+Sden(9)*(x.^(-8));
  S=num./den;
  xw(:,2)=S.*(x.^(a-1/2)).*xw(:,2)/exp(1);
end
