% TESTKLINV Test of inverse KLtransform using Gauss integration
%
f1='%8.2f %5.0f %10.2e %16.8e %16.8e %10.2e %10.2e\n';
f2='%8.2f %5.0f %10.2e %20.12e %10.2e %10.2e\n';
f3='%46.12e\n';
global a b
a=1;
eps0=.5e-8;
%eps0=.5e-12;
fprintf('eps0=%12.4e\n',eps0)
disp('     x       n     err0          Phi             exact         abserr     relerr')
%disp('     x       n     err0        approx/exact       abserr     relerr')
B=24; 
%phiB=(2/pi)*exp(-pi*B/2);
%phiB=(2/pi)*exp(-pi*B/2)*(B+2/pi);
%phiB=(2/pi)^(3/2)*(a^2-1)^(-1/4)*sqrt(B)*exp(-pi*B/2);
%phiB=(1.825/pi)*a^(-1/4)*exp(-pi*B)*(B+1/pi);
phiB=(2/pi)*1.825*a^(-1/4)*exp(-pi*B/4)*(B+4/pi);
%phiB=.5*1.825*b^(-1/4)*exp(-pi*B)*(exp(a*B)/(pi-a) ...
%     +1/pi);
%phiB=0;
%phiB=(1.825/pi)*a^(-1/4)*exp(-B*pi/2)*(B+2/pi);
%phiB=(1.825^2/pi)*(a*b)^(-1/4)*exp(-B*pi/2) ...
%     *(B+2/pi);
%phiB=1.825*a^(-1/4)*exp(-B*pi/2)*((B+2/pi)^2+4/pi^2);
%for x=[.5:.5:10]
for x=[.5 3 5.5 8 10]
%for x=.5
  [Phi,err0,n]=KLinvT(@fKLinv,x,B,phiB,eps0);
%
% exact answers
%
%  ex=.5*pi*x*sinh(a)*exp(-x*cosh(a));
%  ex=.5*pi*exp(-x*cosh(a));
%  ex=sqrt(.5*pi*x)*exp(-a*x);
%  ex=.5*pi*sqrt(a*x)*exp(-a-x)/(a+x);
   ex=(pi/2)^(3/2)*x*exp(-a-x^2/(8*a))/sqrt(a);
%  ex=.5*pi*mfun('BesselK',0,sqrt(x^2+b^2+2*b*x ...
%     *cos(a)));  
%  ex=.5*pi^2*mfun('BesselI',1,x)*mfun('BesselK',1,a);
%  ex=.25*(pi/2)^(3/2)*a*exp(-x-a^2/(8*x))/sqrt(x);
%  ex=.25*pi^2*exp(-.5*x*(a/b+b/a+a*b/x^2));
%  v=sqrt(x^2+4*a*b);
%  ex=.5*pi^2*x*exp(-(a+b)*v/(2*sqrt(a*b)))/v;
%   v=sqrt(x^2+a^2);
%   ex=2*pi^2*(a*x/(2*v))^2*besselk(2,v);
  abserr=abs(Phi-ex); relerr=abserr/abs(ex);
  fprintf(f1,x,n,err0,Phi,ex,abserr,relerr)
%  fprintf(f2,x,n,err0,Phi,abserr,relerr)
  fprintf(f3,ex)
end