% KLT  Kontorovich-Lebedev transform
%
function [Fr,Fi,n1,n2,n3]=KLT(fKL,K0f,nu,eps0)
nur=real(nu); nui=imag(nu); 
n1=20; n2=30; n3=50;
nuilarge=0; nut=nu; nurt=nur; nuit=nui;
Gnut=exp(mfun('lnGAMMA',nut));
Gmnut=exp(mfun('lnGAMMA',-nut));
if nui>10
  nuilarge=1; nut=nur+10*i;
  nurt=real(nut); nuit=imag(nut);
  Gnut=exp(mfun('lnGAMMA',nut));
  Gmnut=exp(mfun('lnGAMMA',-nut));
end
ablag=r_laguerre(400); ableg=r_jacobi01(200);
%
% first integral
%
I0=1; I1=0; 
while abs(I1-I0)>eps0*max(abs(I1),1)
  n1=n1+5;
  if n1>185
    disp('n1 exceeds 185')
    break
  end
  I0=I1;
  xw=gauss(n1,ablag); x=xw(:,1); w=xw(:,2);
  [Kr,Ki]=macdonald(2+x,nurt,nuit); K=Kr+Ki*i;
  f=feval(fKL,2+x);
  I1=sum(w.*exp(x).*K.*f);
end
y1=I1;
%
% second integral
%
I0=1; I1=0; 
while abs(I1-I0)>eps0*max(abs(I1),1)
  n2=n2+5;
  if n2>200
    disp('n2 exceeds 200')
    break
  end
  I0=I1;
  xw=gauss(n2,ableg); x=xw(:,1); w=xw(:,2);
  K0=.5*((1./x).^nut*Gnut+(1./x).^(-nut)*Gmnut);
  [Kr,Ki]=macdonald(2*x,nurt,nuit); K=Kr+Ki*i;
  f=feval(fKL,2*x);
  I1=2*sum(w.*(K-K0).*f);
end
y2=I1;
%
% third integral
%
I0=1; I1=0; 
while abs(I1-I0)>eps0*max(abs(I1),1)
  if n3<=150
    n3=n3+20;
  else
    n3=n3+40;
  end
  if n3>350
%  if n3>175
    n3=350;
%    n3=175;
    disp('n3 exceeds 350')
%    disp('n3 exceeds 175')
    break
  end
  I0=I1;
  xw=gauss(n3,ablag); x=xw(:,1); w=xw(:,2);
  K0=.5*(exp(nut*x)*Gnut+exp(-nut*x)*Gmnut);
  f=feval(fKL,2*exp(-x));
  I1=2*sum(w.*K0.*f);
%
% alternate evaluation of the third integral
%
%  g=feval(K0f,x,nut,Gnut,Gmnut);
%  I1=2*sum(w.*g);
end
y3=I1;
if nuilarge==1
  Gnu=exp(mfun('lnGAMMA',nu));
  Gmnu=exp(mfun('lnGAMMA',-nu));
%
% first integral
%
  I0=1; I1=y1; n1=n1-5;
  while abs(I1-I0)>eps0*max(abs(I1),1)
    n1=n1+5;
    if n1>350
      disp('n1 exceeds 350')
      break
    end
    I0=I1;
    xw=gauss(n1,ablag); x=xw(:,1); w=xw(:,2);
    [Kr,Ki]=macdonald(2+x,nur,nui); K=Kr+Ki*i;
    f=feval(fKL,2+x);
    I1=sum(w.*exp(x).*K.*f);
  end
  y1=I1;
%
% second integral
%
  I0=1; I1=y2; n2=n2-5;
  while abs(I1-I0)>eps0*max(abs(I1),1)
    n2=n2+5;
    if n2>200
      disp('n2 exceeds 200')
      break
    end
    I0=I1;
    xw=gauss(n2,ableg); x=xw(:,1); w=xw(:,2);
    K0=.5*((1./x).^nu*Gnu+(1./x).^(-nu)*Gmnu);
    [Kr,Ki]=macdonald(2*x,nur,nui); K=Kr+Ki*i;
    f=feval(fKL,2*x);
    I1=2*sum(w.*(K-K0).*f);
  end
  y2=I1;
%
% third integral
%
  I0=1; I1=y3; n3=n3-5;
  while abs(I1-I0)>eps0*max(abs(I1),1)
    n3=n3+5;
    if n3>400
%    if n3>180
      disp('n3 exceeds 400')
%      disp('n3 exceeds 180')
      break
    end
    I0=I1;
    xw=gauss(n3,ablag); x=xw(:,1); w=xw(:,2);
    K0=.5*(exp(nu*x)*Gnu+exp(-nu*x)*Gmnu);
    f=feval(fKL,2*exp(-x));
    I1=2*sum(w.*K0.*f);
%
% alternate evaluation of the third integral
%
%    g=feval(K0f,x,nu,Gnu,Gmnu);
%    I1=2*sum(w.*g);
  end
  y3=I1;
end
F=y1+y2+y3; Fr=real(F); Fi=imag(F);