% TESTSR_RIMKLTP  Test the routines sr_RMKLTp.m and sr_IMKLTp.m
%
f1='%8.0f %22.14e %22.14e %22.14e\n';
N=40;
a=0;
for b=[.01 .1 .5 1 3 5 10 15]
  dig=70
%  dig=75
%  dig=80
%  ab=sr_RMKLTp(dig,N,a,b);
%  ab=sr_SQRMKLTp(dig,N,a,b);
  lambda=.75;
  ab=sr_MLRMKLTp(dig,N,a,b,lambda);
%  ab=sr_IMKLTp(dig,N,a,b);
%  ab=sr_IMKLTp(dig,N,a,b);
  ab=vpa(ab,16);
  ab=double(ab);
%
% Example 3.5 (activate sr_RMKLTp on line 10)
%
%  abl=r_jacobi01(N);
%  c=2*sqrt(pi/(b*sinh(pi*b)));
%  g=imag(mfun('lnGAMMA',a+b*i));
%  for n=3:3:15
%    xw=gauss(n,ab);
%    x=xw(:,1); w=xw(:,2);
%    xwl=gauss(n,abl);
%    xl=xwl(:,1); wl=xwl(:,2);
%    I3=c*sum(w.*exp(-(2*x).^2)-wl.*exp(-(2*xl).^2));
%    RI0=sum(w);
%    RI0ex=1/(a+1)+((a+1)*cos(g)-b*sin(g))/ ...
%      ((a+1)^2+b^2);
%    fprintf(f1,n,I3,RI0,RI0ex)
%  end
%  fprintf('\n')
%
% Example 3.6 (activate sr_SQRMKLTp on line 11)
%
%  abl=r_jacobi01(N,0,1/2);
  abj=r_jacobi01(N,0,-lambda);
  c=2*sqrt(pi/(b*sinh(pi*b)));
  g=imag(mfun('lnGAMMA',a+b*i));
  for n=5:5:40
    xw=gauss(n,ab);
    x=xw(:,1); w=xw(:,2);
%    xwl=gauss(n,abl);
    xwj=gauss(n,abj);
%    xl=xwl(:,1); wl=xwl(:,2);
    xj=xwj(:,1); wj=xwj(:,2);
%    ab0=sr_RMKLTp(dig,1,a,b);    
%    ab0=vpa(ab0);
%    ab0=double(ab0);
%    I3=c*(ab0(1,2)-1-sqrt(2)*sum(w.*erf(sqrt(2*x))./ ...
%       sqrt(2*x)-wl.*erf(sqrt(2*xl))./sqrt(2*xl)));
    I3=c*sum(w.*x.^lambda.*besselk(lambda,2*x)-wj.* ...
       xj.^lambda.*besselk(lambda,2*xj));
%    SQI0=sum(w);
    MLI0=sum(w);
%    SQI0ex=1/(1.5-a)+((1.5-a)*cos(g)-b*sin(g))/ ...
      ((1.5-a)^2+b^2);
    MLI0ex=1/(1-a-lambda)+((1-a-lambda)*cos(g)-b*sin(g))/ ...
      ((1-a-lambda)^2+b^2);
%    fprintf(f1,n,I3,SQI0,SQI0ex)
    fprintf(f1,n,I3,MLI0,MLI0ex)
  end
end