% KLINVT Inverse Kontorovich-Lebedev transform
%
function [Phi,err0,n]=KLinvT(fKLinv,x,B,phiB,eps0)
ab=r_jacobi01(200);
n=0; Phi=1; Phi0=0;
while abs(Phi-Phi0)>eps0*abs(Phi)
  n=n+1; Phi0=Phi;
  if n>200
    disp('n exceeds 200')
    break
  end
  xw=gauss(n,ab);
  Phi=0;
  for k=1:n
    xg=B*xw(k,1); wg=xw(k,2);
    f=feval(fKLinv,xg);
    Phi=Phi+wg*macdonald(x,0,xg)*f;
  end
  Phi=B*Phi;
end
err0=1.825*x^(-1/4)*phiB/abs(Phi);