% RMKLINVT Inverse of real modified Kontorovich-Lebedev transform
%
function [Phi,err0,n]=RMKLinvT(fMKLinv,x,B,phiB0,phiB1,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(fMKLinv,xg);
    Kr=macdonald(x,1/2,xg);
    Phi=Phi+wg*Kr*f;
  end
  Phi=B*Phi;
end
err0=(.535*x^(-1/4)*phiB0+sqrt(2*pi/x)*exp(-x)* ...
      phiB1)/abs(Phi);