% TESTRMKLINV Test of inverse real modified KLtransform using Gauss integration
%
%f1='%8.2f %5.0f %10.2e %16.8e %16.8e %10.2e %10.2e\n';
f1='%8.2f %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 y
%a=.5*pi; 
a=pi/3;
eps0=.5e-8;
%eps0=.5e-12;
fprintf('eps0=%12.4e\n',eps0)
%disp('     x       n     err0          Phi             exact         abserr     relerr')
disp('     x        y      n      err0          Phi             exact         abserr     relerr')          
%disp('     x       n     err0        approx/exact       abserr     relerr')
B=10; 
abserm=0; relerm=0;
for x=[1:.2:2 2.5:.5:5 6:10]
  for y=[1:.2:2 2.5:.5:5 6:10]
    if y<=x
      phiB0=(4/pi)*(exp(-2*a*B)/(1-exp(-2*a*B)))*((.535/ ...
        (2*a))*y^(-3/4)*((B+1/(2*a))^2+1/(4*a^2))+(1/ ...
        (2*a+pi/2))*sqrt(2*pi/y)*exp(-y)*exp(-pi*B/2) ...
        *(B+1/(2*a+pi/2))); 
      phiB1=(4/pi)*(exp(-(2*a+pi/2)*B)/(1-exp(-a*B))) ...
        *((.535/(2*a+pi/2))*y^(-3/4)*(B+1/(2*a+pi/2)) ...
        +(1/(2*a+pi))*sqrt(2*pi/y)*exp(-y)*exp(-pi*B/2));
      [Phi,err0,n]=RMKLinvT(@fMKLinv,x,B,phiB0,phiB1,eps0);
%      ex=besselk(0,x+y)+besselk(1,x+y);
      temp=sqrt(x^2+y^2+x*y);
      ex=sqrt(3)*besselk(0,temp)+sqrt(3)*(x+y)*besselk(1,temp)/temp;
      abserr=abs(Phi-ex); relerr=abserr/abs(ex);
      if abserr>abserm, abserm=abserr; end
      if relerr>relerm, relerm=relerr; end
      fprintf(f1,x,y,n,err0,Phi,ex,abserr,relerr)
%      fprintf(f1,x,n,err0,Phi,ex,abserr,relerr)
%      fprintf(f2,x,n,err0,Phi,abserr,relerr)
%      fprintf(f3,ex)
    end
  end
end
[abserm relerm]