% MACDONALD The Macdonald function K_{a+ib}(x) for vector-valued x
%    with positive components and scalar-valued a with 0<=a<2 and
%    b with |b|<=10. The output Kr, Ki are the real and imaginary
%    part, respectively, of K_{a+ib}(x), and canr, cani the degree 
%    of cancellation occurring in the real and imaginary part of K
%    when the argument value is less than 1 and larger than 0.0001.
%
function [Kr,Ki,canr,cani]=macdonald(x,a,b)
s=find(x<=0); if isempty(s)==0, error('some x are nonpositive'), end
if (a<0|a>=2), error('a not in range'), end
Kr=zeros(size(x)); Ki=zeros(size(x)); 
canr=zeros(size(x));
cani=zeros(size(x));
if abs(b)>10
  K=mfun('BesselK',a+b*i,x);
  Kr=real(K); Ki=imag(K);
  return
end
high=find(x>=1); low=find(x<1&x>0.01);
lowlow=find(x<=0.01);
xh=x(high); xl=x(low); xll=x(lowlow);
load -ascii xwmacdonald;
xg=xwmacdonald(:,1); wg=xwmacdonald(:,2);
if isempty(high)==0
  const=.5*(2./xh).^a.*exp(1-.5*xh);
  yr=zeros(size(xh)); yi=zeros(size(xh));
  for nu=1:30
    h=1+(.5*xh-1)*exp(-xg(nu));
    yr=yr+wg(nu)*exp(a*xg(nu)-.25*xh.^2*exp(-xg(nu))./h) ...
    .*h.^(-a-1).*(h.^(2*a)+(.5*xh*exp(-xg(nu))).^(2*a)) ...
    .*cos(b*(xg(nu)+log(2*h./xh)));
    yi=yi+wg(nu)*exp(a*xg(nu)-.25*xh.^2*exp(-xg(nu))./h) ...
    .*h.^(-a-1).*(h.^(2*a)-(.5*xh*exp(-xg(nu))).^(2*a)) ...
    .*sin(b*(xg(nu)+log(2*h./xh)));
  end
  Kr(high)=const.*yr; Ki(high)=const.*yi;
  if a==0|b==0, Ki(high)=zeros(size(xh)); end
end
if isempty(low)==0
  load -ascii xwleg01;
  xgl=xwleg01(:,1); wgl=xwleg01(:,2);
  const=.5*(2./xl).^a*exp(.5);
  yr=zeros(size(xl)); yi=zeros(size(xl));
  for nu=1:30
    h1=1-.5*exp(-xg(nu));
    yr=yr+wg(nu)*exp(a*xg(nu)-.25*xl.^2*exp(-xg(nu))/h1) ...
    *h1^(-a-1).*(h1^(2*a)+(.5*xl*exp(-xg(nu))).^(2*a)) ...
    .*cos(b*(xg(nu)+log(2*h1./xl)));
    yi=yi+wg(nu)*exp(a*xg(nu)-.25*xl.^2*exp(-xg(nu))/h1) ...
    *h1^(-a-1).*(h1^(2*a)-(.5*xl*exp(-xg(nu))).^(2*a)) ...
    .*sin(b*(xg(nu)+log(2*h1./xl)));
  end
  yr=const.*yr; yi=const.*yi;
  c=log(1./xl);
  zr=zeros(size(xl)); zi=zeros(size(xl));
  for nu=1:30
    zr=zr+wgl(nu)*exp(-xl.*cosh(xgl(nu).*c)) ...
    .*cosh(a*xgl(nu).*c).*cos(b*xgl(nu).*c);
    zi=zi+wgl(nu)*exp(-xl.*cosh(xgl(nu).*c)) ...
    .*sinh(a*xgl(nu).*c).*sin(b*xgl(nu).*c);
  end
  zr=c.*zr; zi=c.*zi;
  Kr(low)=yr+zr; Ki(low)=yi+zi; 
  canr(low)=log10(max(abs(yr),abs(zr))./abs(yr+zr));
  if abs(yi+zi)>0, cani(low)=log10(max(abs(yi),abs(zi)) ...
  ./abs(yi+zi)); end
  if a==0|b==0, Ki(low)=zeros(size(xl)); end
end
if isempty(lowlow)==0
  ga=.57721566490153286; nu=a+b*i;
  if nu==0
    Kr(lowlow)=(log(2./xll)-ga).*(1+.25*xll.^2 ...
    +.015625*xll.^4)+.25*xll.^2+3*xll.^4/128;
  elseif nu==1
    Kr(lowlow)=1./xll-.5*xll.*log(2./xll).*(1+.125*xll.^2 ...
       +xll.^4/192)-.25*xll.*(1-2*ga+.125*(2.5-2*ga) ...
       *xll.^2+(10/3-2*ga)*xll.^4/192);
  else
    g1=exp(mfun('lnGAMMA',nu));
    gm1=exp(mfun('lnGAMMA',-nu));
    K=.5*(xll/2).^(-nu).*(g1*(1+.25*xll.^2/(1-nu) ...
      +.03125*xll.^4/((1-nu)*(2-nu))) ...
      +(xll/2).^(2*nu)*gm1.*(1+.25*xll.^2/(1+nu) ...
      +.03125*xll.^4/((1+nu)*(2+nu))));
    Kr(lowlow)=real(K); Ki(lowlow)=imag(K); 
  end
  if a==0|b==0, Ki(lowlow)=zeros(size(xll)); end
end