subroutine qmacdonald(x,a,b,Kr,Ki,canr,cani)
      real*16 x,a,b,Kr,Ki,canr,cani,xg,wg,xgl,wgl,const,yr,yi,
     *h,c,zr,zi,yzrm,yzim,canr,cani
      common/gauss/xg(55),wg(55),xgl(60),wgl(60)
    1 format(/)
      if(a.lt.0.q0 .or. a.gt.2.q0) then
        write(*,2)
    2   format('a not in range')
        stop
      end if
      if(qabs(b).gt.10.q0) then
        write(*,3)
    3   format('b not in range')
        stop
      end if
      canr=0.q0
      cani=0.q0
      if(x.ge.1.q0) then
        const=.5q0*(2.q0/x)**a*qexp(1.0q0-.5q0*x)
        yr=0.q0
        yi=0.q0
        do 20 k=1,55
          h=1.q0+(.5q0*x-1.q0)*qexp(-xg(k))
          yr=yr+wg(k)*qexp(a*xg(k)-.25q0*x**2*qexp(-xg(k))
     *      /h)*h**(-a-1.q0)*(h**(2.q0*a)+(.5q0*x*
     *      qexp(-xg(k)))**(2.q0*a))*qcos(b*(xg(k)+
     *      qlog(2.q0*h/x)))
          yi=yi+wg(k)*qexp(a*xg(k)-.25q0*x**2*qexp(-xg(k))
     *      /h)*h**(-a-1.q0)*(h**(2.q0*a)-(.5q0*x*
     *      qexp(-xg(k)))**(2.q0*a))*qsin(b*(xg(k)+
     *      qlog(2.q0*h/x)))
   20   continue
        Kr=const*yr
        Ki=const*yi
      else
        const=.5q0*(2.q0/x)**a*qexp(.5q0)
        yr=0.q0
        yi=0.q0
        do 30 k=1,55
          h=1.q0-.5q0*qexp(-xg(k))
          yr=yr+wg(k)*qexp(a*xg(k)-.25q0*x**2*qexp(-xg(k))
     *      /h)*h**(-a-1.q0)*(h**(2.q0*a)+(.5q0*x*
     *      qexp(-xg(k)))**(2.q0*a))*qcos(b*(xg(k)+
     *      qlog(2.q0*h/x)))
          yi=yi+wg(k)*qexp(a*xg(k)-.25q0*x**2*qexp(-xg(k))
     *      /h)*h**(-a-1.q0)*(h**(2.q0*a)-(.5q0*x*
     *      qexp(-xg(k)))**(2.q0*a))*qsin(b*(xg(k)+
     *      qlog(2.q0*h/x)))
   30   continue
        yr=const*yr
        yi=const*yi
        c=qlog(1.q0/x)
        zr=0.q0
        zi=0.q0
        do 40 k=1,60
          zr=zr+wgl(k)*qexp(-x*qcosh(xgl(k)*c))*qcosh(a*    
     *      xgl(k)*c)*qcos(b*xgl(k)*c)
          zi=zi+wgl(k)*qexp(-x*qcosh(xgl(k)*c))*qsinh(a*    
     *      xgl(k)*c)*qsin(b*xgl(k)*c)
   40   continue
        zr=c*zr
        zi=c*zi
        yzrm=qabs(zr)
        if(qabs(yr).gt.yzrm) yzrm=qabs(yr)
        yzim=qabs(zi)
        if(qabs(yi).gt.yzim) yzim=qabs(yi)
        canr=qlog10(yzrm/qabs(yr+zr))
        if(qabs(yi+zi).gt.0.q0) cani=qlog10(yzim/qabs(yi+zi))
        Kr=yr+zr
        Ki=yi+zi
      end if
      end