c
c
      program qOPbess
c
c
      external qquad
      real*16 qnu,x0,a1(100),b1(100),a2(100),b2(100),
     *a3(100),b3(100),a4(100),b4(100),e(100),q1mach,
     *qepsma,dendl(4),dendr(4),qeps,dxp(1),dyp(1),
     *dxfer(100),dwfer(100),dalpha(40),dbeta(40),dbe(40),
     *dx(100),dw(100),dxm(400),dwm(400),dp0(400),dp1(400),
     *dp2(400)
      logical finld,finrd
      common/q/qnu,x0,a1,b1,a2,b2,a3,b3,a4,b4,e
c
c This routine generates in quadruple precision the first 40 
c recurrence coefficients of the orthogonal polynomials for 
c the K-Bessel weight function. Printed are the quadruple-precision
c values of the alpha- and beta-coefficients. The routine
c requires the routines q1mach, qrecur, qgauss, qmcdis, qqgp,
c qsti, and qlancz, which are quadruple-precision versions of
c the respective (similarly named) ORTHPOL routines.
c
      nmax=40
      qnu=1.q0/3.q0
      if(qnu.lt..05q0 .or.qn.gt..95q0) then
        write(*,1)
    1   format(1x,'nu out of range')
        stop
      end if
      x0=1.q0
      call qrecur(100,6,0.q0,-qnu,a1,b1,ier1)
      call qrecur(100,6,0.q0,qnu,a2,b2,ier2)
      call qrecur(100,1,0.q0,0.q0,a3,b3,ier3)
      call qrecur(100,7,0.q0,0.q0,a4,b4,ier4)
      if(ier1.ne.0 .or. ier2.ne.0 .or. ier3.ne.0 .or.
     *  ier4.ne.0) then
        write(*,2) ier1,ier2,ier3,ier4
    2   format(1x,'ier1 in qrecur=',i5,'  ier2=',i5,'  ier3=',i5,
     *    '  ier4=',i5)
        stop
      end if
      finld=.true.
      finrd=.false.
      qepsma=q1mach(3)
c
c qepsma is the machine quadruple precision.
c
      iq=1
      idelta=1
      irout=1
      mcd=4
      mp=0
      ncpmd=100
c
c Set up the partition for the discretization of the inner product.
c
      dendl(1)=0.q0
      dendl(2)=1.q0
      dendl(3)=2.q0
      dendl(4)=3.q0
      dendr(1)=1.q0
      dendr(2)=2.q0
      dendr(3)=3.q0
      qeps=1.q5*qepsma
c
c Compute the desired recursion coefficients by the multiple-component
c discretization procedure. 
c
      call qmcdis(nmax,ncpmd,mcd,mp,dxp,dyp,qquad,qeps,iq,idelta,irout,
     *  finld,finrd,dendl,dendr,dxfer,dwfer,dalpha,dbeta,ncapd,kountd,
     *  ierrd,ied,dbe,dx,dw,dxm,dwm,dp0,dp1,dp2)
c
c Print the results.
c
      write(*,3)
    3 format(/5x,'k',13x,'dalpha(k)',29x,'dbeta(k)'/)
      do 20 k=1,nmax
        km1=k-1
        write(*,4) km1,dalpha(k),dbeta(k)
    4   format(1x,i5,2d36.28)
   20 continue
      stop
      end 

      subroutine qquad(n,dx,dw,i,ierr)
      real*16 dx(n),dw(n),qnu,x0,a1(100),b1(100),a2(100),
     *b2(100),a3(100),b3(100),a4(100),b4(100),e(100),z(87),
     *w(87),pi,c,eps,q1mach,x1,c1,qgamma,x,x2,sum1,t1,qm,
     *sum,t,qbessK
      common/q/qnu,x0,a1,b1,a2,b2,a3,b3,a4,b4,e
      common/qa/z,w
      pi=4.q0*qatan(1.q0)
      c=2.q0*qcos(.5q0*qnu*pi)/pi
      eps=q1mach(3)
      x1=10.q0
      if(i.eq.1) then
        c1=pi*((.25q0*x0)**(1.q0-qnu))/(qsin(qnu*pi)*qgamma(1.q0-
     *    qnu,ie))
        call qgauss(n,a1,b1,eps,dx,dw,ierg,e)
        do 20 k=1,n
          dx(k)=.5q0*x0*(1.q0+dx(k))
          x=dx(k)
          m=0
          x2=x**2
          sum1=1.q0
          t1=1.q0
   10     m=m+1
          qm=qreal(m)
          sum=sum1
          t=t1
          t1=x2*t/(4.q0*qm*(qm-qnu))
          sum1=sum+t1
          if(sum1.ne.sum) goto 10
          dw(k)=c*c1*sum1*dw(k)
   20   continue
      else if(i.eq.2) then
        c1=pi*((.25q0*x0)**(1.q0+qnu))/(qsin(qnu*pi)*qgamma(1.q0+
     *    qnu,ie))
        call qgauss(n,a2,b2,eps,dx,dw,ierg,e)
        do 40 k=1,n
          dx(k)=.5q0*x0*(1.q0+dx(k))
          x=dx(k)
          m=0
          x2=x**2
          sum1=1.q0
          t1=1.q0
   30     m=m+1
          qm=qreal(m)
          sum=sum1
          t=t1
          t1=x2*t/(4.q0*qm*(qm+qnu))
          sum1=sum+t1
          if(sum1.ne.sum) goto 30
          dw(k)=-c*c1*sum1*dw(k)
   40   continue
      else if(i.eq.3) then
        call qgauss(n,a3,b3,eps,dx,dw,ierg,e)
        do 50 k=1,n
          dx(k)=.5q0*((x1-x0)*dx(k)+x1+x0)
          dw(k)=.5q0*c*(x1-x0)*qexp(-dx(k))*qbessK(dx(k),qnu)*dw(k)
   50   continue
      else 
        c1=qexp(-x1)
        call qgauss(n,a4,b4,eps,dx,dw,ierg,e)
        do 60 k=1,n
          dx(k)=x1+dx(k)
          dw(k)=c*c1*qbessK(dx(k),qnu)*dw(k)
   60   continue
      end if
      return
      end

      real*16 function qbessK(x,qnu)
c
c This computes exp(x)K_nu(x) in quadruple precision for 
c positive x and 0.05 <= nu <= 0.95.
c
      real*16 x,qnu,z(87),w(87),eps,q1mach,al,a(87),
     *b(87),e(87),pi,c1,qser,qgamma,c2,sum
      common/qa/z,w
      if(x.le.0.q0) then
        write(*,1)
    1   format(1x,'invalid x-argument in qbessK')
        return
      end if
      eps=q1mach(3)
      al=qnu-.5q0
      call qrecur(87,7,al,0.q0,a,b,ierr)
      call qgauss(87,a,b,eps,z,w,ierg,e)
      pi=4.q0*qatan(1.q0)
      c1=pi/(2.q0*qsin(qnu*pi))
      if(x.le.2.q0) then
        qbessK=c1*qexp(x)*(((.5q0*x)**(-qnu))*qser(x,-qnu)/
     *    qgamma(1.q0-qnu,ie)-((.5q0*x)**qnu)*qser(x,qnu)/
     *    qgamma(1.q0+qnu,ie))
      else
        n=87
        c2=qsqrt(pi)/((2.q0**qnu)*qgamma(qnu+.5q0,ie))
        sum=0.q0
        do 10 k=1,n
          sum=sum+w(k)*(2.q0+z(k)/x)**al
   10   continue
        qbessK=c2*sum/qsqrt(x)
      end if
      return
      end

      real*16 function qser(x,qnu)
      real*16 x,qnu,x2,sum1,t1,qk,sum,t
      k=0
      x2=x**2
      sum1=1.q0
      t1=1.q0
   10 k=k+1
      qk=qreal(k)
      sum=sum1
      t=t1
      t1=x2*t/(4.q0*qk*(qnu+qk))
      sum1=sum+t1
      if(sum1.ne.sum) goto 10
      qser=sum1
      return
      end
  
      real*16 function qwf(x,i)
      real*16 x
      return
      end

