% TESTKLT Test of KL-transform
%
f1='%8.2f %17.9e %17.9e %10.2e %10.2e %4.0f %4.0f %4.0f\n';
eps0=.5e-10;
fprintf('eps0=%12.4e\n',eps0)
disp('     b           Fr                Fi            err        eri      ni   n2   n3')
global a lambda 
for a=0;
  for lambda=3/4;
    for b=[.2 .5 1 4 7 10 15]
      nu=b*i;
      [Fr,Fi,n1,n2,n3]=KLT(@fKL,@K0f,nu,eps0);
%
% exact answers
%
      ex=.5*pi/cosh(.5*pi*b);
%      ex=.5*pi*b/sinh(.5*b*pi);
%      ex=2^(lambda-1)*exp(2*real(mfun('lnGAMMA',(lambda+1+i*b)/2)));
%      ex=(pi/sin(a))*sinh(a*b)/sinh(pi*b);
%      ex=(pi/sinh(a))*sin(a*b)/sinh(pi*b);
%      ex=pi*b/sinh(pi*b);
%      ex=(2^lambda)*(gamma(lambda+1)/gamma(2*lambda+2)) ...
%         *exp(2*real(mfun('lnGAMMA',lambda+1+i*b)));
%      ex=sqrt(pi/2)*pi/cosh(pi*b);
%      K=macdonald(a,0,b);
%      ex=(pi/cosh(pi*b))*K/sqrt(a);
%      K=macdonald(1/(8*a^2),0,b/2);
%      ex=(sqrt(pi)/(4*a))*exp(1/(8*a^2))*K/cosh(pi*b/2);
%      K=macdonald(a,0,b);
%      ex=K^2;
%      P=mehler(cosh(a),b,0,eps0);
%      ex=sqrt(pi/2)*(pi/cosh(pi*b))*P;
%      K=macdonald(a,0,b);
%      ex=(pi/cosh(pi*b))*K/sqrt(a);
%      ex=(pi/(2*sqrt(2)))*cos(b*acosh(sqrt(2))) ...
%         /cosh(pi*b/2);
%      ex=(pi/(2*sqrt(2)))*sin(b*acosh(sqrt(2))) ...
%         /sinh(pi*b/2);
%      ex=(pi/(2*cosh(a/2)))*cos(b*a/2)/cosh(pi*b/2);
%      ex=(pi/(2*cosh(a/2)))*sin(b*a/2)/sinh(pi*b/2);
%      ex=(pi/2)*sin(a*b/2)/(b*cosh(pi*b/2));
%      ex=.5*pi*(1-cos(a*b/2))/(b*sinh(pi*b/2));
%      ex=.5*pi*(b*(1-cos(.5*a*b)*cosh(a/2))+ ...
%          sin(.5*a*b)*sinh(a/2))/(b*(1+b^2)* ...
%          cosh(.5*pi*b));
%      ex=(pi/sinh(pi*b))*(cosh(lambda*b)*sin(a*b) ...
%         *cos(lambda)*sinh(a)/(sinh(a)^2*cos(lambda)^2+ ...
%         cosh(a)^2*sin(lambda)^2)+sinh(lambda*b)*cos(a*b) ...
%         *sin(lambda)*cosh(a)/(sinh(a)^2*cos(lambda)^2+ ...
%         cosh(a)^2*sin(lambda)^2));
%      ex=(pi/(2*cosh(.5*pi*b)))*mehler(cosh(a),.5*b,0,eps0);
%      ex=.5*pi^2/(cosh(pi*b)+cos(lambda*pi));
%      ex=pi*(sqrt(2)*cosh(.5*pi*b)-1)/(2*cosh(pi*b) ...
%         *cosh(.5*pi*b));
      err=abs(Fr-ex);
      eri=abs(Fi-imag(ex));
      fprintf(f1,b,Fr,Fi,err,eri,n1,n2,n3)
    end
  end
end