
C     OBTAINED FROM ILKKA KARASALO 8/10 1990

C     Sven Ivansson, July 2002: BESSJ(N,X) and BESSH(N,X) 


C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     DBES0RAT
C     Bessel function J0 of a REAL*8 argument
C     J.F. Hart, et.al: 'Computer Approximations', 
C     John Wiley, 1968, p. 146-147.
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      real*8 function dbes0rat(x)
      implicit real*8 (a-h,o-z)

c     0 <= X <= 8:
c     Coefficients of rational approximation P0(X**2)/Q0(X**2)
c     to J0(X),  (Table 5845 p. 300).
      parameter (
     *  p0 = +0.1859623176218978D+18,
     *  p1 = -0.4414582939181598D+17,
     *  p2 = +0.2334489171877870D+16,
     *  p3 = -0.4776555944267359D+14,
     *  p4 = +0.4621722250317180D+12,
     *  p5 = -0.2271490439553603D+10,
     *  p6 = +0.5513584564770752D+7,
     *  p7 = -0.5292617130384557D+4,

     *  q0 = +0.1859623176218977D+18,
     *  q1 = +0.2344750013658997D+16,
     *  q2 = +0.1501546244976975D+14,
     *  q3 = +0.6439867453513326D+11,
     *  q4 = +0.2042514835213436D+9,
     *  q5 = +0.4940307949181397D+6,
     *  q6 = +0.8847203675617550D+3,
     *  q7 = +1.0000000000000000D+0
     *   )


c     X > 8:
c     Coefficients of rational approximations of P0(X) and Q0(X)
c     in J0(X) = SQRT(2/PI/X)*(P0(X)*COS(X0) - Q0(X)*SIN(X0)) 
c     where X0 = X - PI/4:

c     Rational approximation of P0(X) = R(Z**2)/S(Z**2) where
c     Z = 8/X,    (Table 6546, p.321).
      parameter (
     *  r0 = +0.8554822541506661D+4,
     *  r1 = +0.8894437532960619D+4,
     *  r2 = +0.2204501043965180D+4,
     *  r3 = +0.1286775857487142D+3,
     *  r4 = +0.9004793474802880D+0,

     *  s0 = +0.8554822541506663D+4,
     *  s1 = +0.8903836141709595D+4,
     *  s2 = +0.2214048851914710D+4,
     *  s3 = +0.1308849004999239D+3,
     *  s4 = +1.0000000000000000D+0
     *   )

c     Rational approximation of Q0(X) = Z*U(Z**2)/V(Z**2), where
c     Z = 8/X,    (Table 6946, p.328).
      parameter (
     *  u0 = -0.3751053495495711D+2,
     *  u1 = -0.4609382681462517D+2,
     *  u2 = -0.1399097686596068D+2,
     *  u3 = -0.1049732798234555D+1,
     *  u4 = -0.9352595329403189D-2,

     *  v0 = +0.2400674237117267D+4,
     *  v1 = +0.2971983745208492D+4,
     *  v2 = +0.9215669755265309D+3,
     *  v3 = +0.7442838974141118D+2,
     *  v4 = +1.0000000000000000D+0
     *   )

      parameter (
     *  pibyfour =  .7853981633974483D+0,   ! PI/4
     *  obyfopi  =  .7957747154594768D-1    ! 1/4/PI
     *   )
C------------------------------------------------------------
C------------------------------------------------------------
      ax = abs(x)

      if(ax.le.8.d0) then
       y = x*x
       dbes0rat = 
     *      (p0 + y*(p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 +
     *            y*(p6 + y*p7))))))) /
     *      (q0 + y*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 +
     *            y*(q6 + y*q7))))))) 

      else
       z = 8.d0/ax
       y = z*z
       xx = ax - pibyfour
       dbes0rat = sqrt(obyfopi*z) * 
     *     (
     *       cos(xx) * (r0 + y*(r1 + y*(r2 + y*(r3 + y*r4)))) /
     *                 (s0 + y*(s1 + y*(s2 + y*(s3 + y*s4))))
     *   - z*sin(xx) * (u0 + y*(u1 + y*(u2 + y*(u3 + y*u4)))) /
     *                 (v0 + y*(v1 + y*(v2 + y*(v3 + y*v4))))  
     *     )
      end if

      return
      end


C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     DHAN0RAT
C     Hankel function H01 of a REAL*8 argument.
C     J.F. Hart, et.al: 'Computer Approximations', 
C     John Wiley, 1968, p. 146-147.
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      complex*16 function dhan0rat(x)
      implicit real*8 (a-h,o-z)

c     0 <= X <= 8:
c     Coefficients of rational approximation P0(X**2)/Q0(X**2)
c     to the regular part Y0R(X) of Y0(X),  (Table 6242 p. 311).
      parameter (
     *  p0 = -0.1018609712241696D+17,
     *  p1 = +0.2438044852569504D+17,
     *  p2 = -0.1902052663082696D+16,
     *  p3 = +0.4770917646459635D+14,
     *  p4 = -0.5283488038967236D+12,
     *  p5 = +0.2869108287298390D+10,
     *  p6 = -0.7531475457816756D+7,
     *  p7 = +0.7705933648729243D+4,

     *  q0 = +0.1380149638637746D+18,
     *  q1 = +0.1786338779946877D+16,
     *  q2 = +0.1177522687238065D+14,
     *  q3 = +0.5216295494381439D+11,
     *  q4 = +0.1716517441570019D+9,
     *  q5 = +0.4334677945014705D+6,
     *  q6 = +0.8183580864178241D+3,
     *  q7 = +1.0000000000000000D+0
     *   )


c     X > 8:
c     Coefficients of rational approximations of P0(X) and Q0(X)
c     in J0(X) = SQRT(2/PI/X)*(P0(X)*COS(X0) - Q0(X)*SIN(X0)) 
c     where X0 = X - PI/4:

c     Rational approximation of P0(X) = R(Z**2)/S(Z**2) where
c     Z = 8/X,    (Table 6546, p.321).
      parameter (
     *  r0 = +0.8554822541506661D+4,
     *  r1 = +0.8894437532960619D+4,
     *  r2 = +0.2204501043965180D+4,
     *  r3 = +0.1286775857487142D+3,
     *  r4 = +0.9004793474802880D+0,

     *  s0 = +0.8554822541506663D+4,
     *  s1 = +0.8903836141709595D+4,
     *  s2 = +0.2214048851914710D+4,
     *  s3 = +0.1308849004999239D+3,
     *  s4 = +1.0000000000000000D+0
     *   )

c     Rational approximation of Q0(X) = Z*U(Z**2)/V(Z**2), where
c     Z = 8/X,    (Table 6946, p.328).
      parameter (
     *  u0 = -0.3751053495495711D+2,
     *  u1 = -0.4609382681462517D+2,
     *  u2 = -0.1399097686596068D+2,
     *  u3 = -0.1049732798234555D+1,
     *  u4 = -0.9352595329403189D-2,

     *  v0 = +0.2400674237117267D+4,
     *  v1 = +0.2971983745208492D+4,
     *  v2 = +0.9215669755265309D+3,
     *  v3 = +0.7442838974141118D+2,
     *  v4 = +1.0000000000000000D+0
     *   )

      parameter (
     *  pibyfour =  .7853981633974483D+0,   ! PI/4
     *  obyfopi  =  .7957747154594768D-1,   ! 1/4/PI
     *  twobypi  =  .6366197723675814D+0,   ! 2/PI
     *  gamma    =  .5772156649015329D+0    ! Eulers constant
     *   )
C------------------------------------------------------------
C     H0(X) = J0(X) + iY0(X) 
C     where 
C     Y0(X) = Y0R(X) + 2/PI*J0(X)*LOG(X)
C
C     This may be written
C     Re(H0(X)) = SIGN(X)*J0(X)  
C     Im(H0(X)) = Y0R(X) + 2/PI*J0(X)*LOG(ABS(X))
C     for all real nonzero X, where Y0R(X) is regular, and
C     is approximated as described above.
C------------------------------------------------------------
      ax = abs(x)


      if(ax.le.8.d0) then

       if(x.eq.0.d0) then
C	"overflow" but choose symmetry in  " J0 = 0.5d0 * ( H0(1) + H0(2) ) "
        DHAN0RAT = DCMPLX(1.0D0,0.D0)              !!! SI
        return
       end if

       besj0 = dbes0rat(x)

       y = x*x
       besy0r = 
     *      (p0 + y*(p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 +
     *            y*(p6 + y*p7))))))) /
     *      (q0 + y*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 +
     *            y*(q6 + y*q7))))))) 


       if(x.gt.0) then
        dhan0rat = dcmplx(besj0,besy0r+twobypi*besj0*log(ax))
       else
        dhan0rat = dcmplx(-besj0,besy0r+twobypi*besj0*log(ax))
       end if


      else
       z = 8.d0/ax
       y = z*z
       xx = ax - pibyfour
       pofy =  (r0 + y*(r1 + y*(r2 + y*(r3 + y*r4)))) /
     *         (s0 + y*(s1 + y*(s2 + y*(s3 + y*s4))))
       qofy =  z*(u0 + y*(u1 + y*(u2 + y*(u3 + y*u4)))) /
     *           (v0 + y*(v1 + y*(v2 + y*(v3 + y*v4))))  

       cosxx = cos(xx)     
       sinxx = sin(xx)
       sqrtz = sqrt(obyfopi*z)
       besj0 = sqrtz * (pofy*cosxx - qofy*sinxx)
       besy0 = sqrtz * (pofy*sinxx + qofy*cosxx)

       if(x.gt.0.d0) then
        dhan0rat = dcmplx(besj0,besy0)
       else
        dhan0rat = dcmplx(-besj0,besy0)
       end if

      end if

      return
      end



C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     DBES1RAT
C     Bessel function J1 of a REAL*8 argument
C     J.F. Hart, et.al: 'Computer Approximations', 
C     John Wiley, 1968, p. 146-147.
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      real*8 function dbes1rat(x)
      implicit real*8 (a-h,o-z)

c     0 <= X <= 8:
c     Coefficients of rational approximation P1(X**2)/Q1(X**2)
c     to J1(X)/X,  (Table 6046 p. 306).
      parameter (
     *  p0 = +0.2214887880421963D+18,
     *  p1 = -0.2512374214703213D+17,
     *  p2 = +0.8482420744781273D+15,
     *  p3 = -0.1249820367262025D+14,
     *  p4 = +0.9316565296724673D+11,
     *  p5 = -0.3686668987022982D+9,
     *  p6 = +0.7437023817119996D+6,
     *  p7 = -0.6079530179607414D+3,

     *  q0 = +0.4429775760843926D+18,
     *  q1 = +0.5124712716484872D+16,
     *  q2 = +0.2989836307725487D+14,
     *  q3 = +0.1158192127466889D+12,
     *  q4 = +0.3281940344534196D+9,
     *  q5 = +0.6988586184485076D+6,
     *  q6 = +0.1077741289433304D+4,
     *  q7 = +1.0000000000000000D+0
     *   )


c     X > 8:
c     Coefficients of rational approximations of P1(X) and Q1(X)
c     in J1(X) = SQRT(2/PI/X)*(P1(X)*COS(X1) - Q1(X)*SIN(X1)) 
c     where X1 = X - 3*PI/4:

c     Rational approximation of P1(X) = R(Z**2)/S(Z**2) where
c     Z = 8/X,    (Table 6747, p.324).
      parameter (
     *  r0 = +0.1290918471896188D+5,
     *  r1 = +0.1309042051103506D+5,
     *  r2 = +0.3132752956355070D+4,
     *  r3 = +0.1743137974837902D+3,
     *  r4 = +0.1228505376435904D+1,

     *  s0 = +0.1290918471896188D+5,
     *  s1 = +0.1306678308784402D+5,
     *  s2 = +0.3109281416770029D+4,
     *  s3 = +0.1690472177500861D+3,
     *  s4 = +1.0000000000000000D+0
     *   )

c     Rational approximation of Q1(X) = Z*U(Z**2)/V(Z**2), where
c     Z = 8/X,    (Table 7147, p.332).
      parameter (
     *  u0 = +0.1446528287499521D+3,
     *  u1 = +0.1744291689092426D+3,
     *  u2 = +0.5173653281836592D+2,
     *  u3 = +0.3799445379698067D+1,
     *  u4 = +0.3636346647603471D-1,

     *  v0 = +0.3085927013332317D+4,
     *  v1 = +0.3734340106016302D+4,
     *  v2 = +0.1119109852704749D+4,
     *  v3 = +0.8522392064341340D+2,
     *  v4 = +1.0000000000000000D+0
     *   )

      parameter (
     *  thrpibyf =  .2356194490192345D+1,   ! 3*PI/4
     *  obyfopi  =  .7957747154594768D-1    ! 1/4/PI
     *   )
C------------------------------------------------------------
C------------------------------------------------------------
      ax = abs(x)

      if(ax.le.8.d0) then
       y = x*x
       dbes1rat = 
     *    x * (p0 + y*(p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 +
     *              y*(p6 + y*p7))))))) /
     *        (q0 + y*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 +
     *              y*(q6 + y*q7))))))) 

      else
       z = 8.d0/ax
       y = z*z
       xx = ax - thrpibyf
       dbes1r = sqrt(obyfopi*z) * 
     *     (
     *       cos(xx) * (r0 + y*(r1 + y*(r2 + y*(r3 + y*r4)))) /
     *                 (s0 + y*(s1 + y*(s2 + y*(s3 + y*s4))))
     *   - z*sin(xx) * (u0 + y*(u1 + y*(u2 + y*(u3 + y*u4)))) /
     *                 (v0 + y*(v1 + y*(v2 + y*(v3 + y*v4))))  
     *     )

       if(x.gt.0.d0) then
        dbes1rat = dbes1r
       else
        dbes1rat = -dbes1r
       end if

      end if

      return
      end


C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     DHAN1RAT
C     Hankel function H11 of a REAL*8 argument.
C     J.F. Hart, et.al: 'Computer Approximations', 
C     John Wiley, 1968, p. 146-147.
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      complex*16 function dhan1rat(x)
      implicit real*8 (a-h,o-z)

c     0 <= X <= 8:
c     Coefficients of rational approximation P1(X**2)/Q1(X**2)
c     to Y1R(X)/X where Y1R(X) is the regular part of Y1(X),  
c     (Table 6441 p. 311).
      parameter (
     *  p0 = -0.1065255800731456D+17,
     *  p1 = +0.2801795134748275D+16,
     *  p2 = -0.1197649702061131D+15,
     *  p3 = +0.1913518877315997D+13,
     *  p4 = -0.1380063196398078D+11,
     *  p5 = +0.4574459090948714D+8,
     *  p6 = -0.5706839261364355D+5,

     *  q0 = +0.5433395942070581D+17,
     *  q1 = +0.7711237672472039D+15,
     *  q2 = +0.5616196558236932D+13,
     *  q3 = +0.2774054192456350D+11,
     *  q4 = +0.1029552862389160D+9,
     *  q5 = +0.2977641751209388D+6,
     *  q6 = +0.6583923936947126D+3,
     *  q7 = +1.0000000000000000D+0
     *   )


c     X > 8:
c     Coefficients of rational approximations of P1(X) and Q1(X)
c     in J1(X) = SQRT(2/PI/X)*(P1(X)*COS(X1) - Q1(X)*SIN(X1)) 
c     where X1 = X - 3*PI/4:

c     Rational approximation of P1(X) = R(Z**2)/S(Z**2) where
c     Z = 8/X,    (Table 6747, p.324).
      parameter (
     *  r0 = +0.1290918471896188D+5,
     *  r1 = +0.1309042051103506D+5,
     *  r2 = +0.3132752956355070D+4,
     *  r3 = +0.1743137974837902D+3,
     *  r4 = +0.1228505376435904D+1,

     *  s0 = +0.1290918471896188D+5,
     *  s1 = +0.1306678308784402D+5,
     *  s2 = +0.3109281416770029D+4,
     *  s3 = +0.1690472177500861D+3,
     *  s4 = +1.0000000000000000D+0
     *   )

c     Rational approximation of Q1(X) = Z*U(Z**2)/V(Z**2), where
c     Z = 8/X,    (Table 7147, p.332).
      parameter (
     *  u0 = +0.1446528287499521D+3,
     *  u1 = +0.1744291689092426D+3,
     *  u2 = +0.5173653281836592D+2,
     *  u3 = +0.3799445379698067D+1,
     *  u4 = +0.3636346647603471D-1,

     *  v0 = +0.3085927013332317D+4,
     *  v1 = +0.3734340106016302D+4,
     *  v2 = +0.1119109852704749D+4,
     *  v3 = +0.8522392064341340D+2,
     *  v4 = +1.0000000000000000D+0
     *   )

      parameter (
     *  obyfopi  =  .7957747154594768D-1,   ! 1/4/PI
     *  thrpibyf =  .2356194490192345D+1,   ! 3*PI/4
     *  twobypi  =  .6366197723675814D+0    ! 2/PI
     *   )
C------------------------------------------------------------
C     H1(X) = J1(X) + iY1(X) 
C     where 
C     Y1(X) = Y1R(X) + 2/PI * (J1(X)*LOG(X) - 1/X)
C
C     This may be written
C     Re(H1(X)) = SIGN(X)*J1(X)  
C     Im(H1(X)) = Y1R(X) + 2/PI*(J1(X)*LOG(ABS(X)) - 1/X)
C     for all real nonzero X, where Y1R(X) is regular, and
C     is approximated as described above.
C------------------------------------------------------------
      ax = abs(x)


      if(ax.le.8.d0) then

       if(x.eq.0.d0) then
        dhan1rat = dcmplx(0.d0,0.d0)   ! Singular point
        return
       end if

       besj1 = dbes1rat(x)

       y = x*x
       besy1r = 
     *    x * (p0 + y*(p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 +
     *              y*p6)))))) /
     *        (q0 + y*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 +
     *              y*(q6 + y*q7))))))) 

       besy1s = twobypi*(besj1*log(ax) - 1.d0/x)

       if(x.gt.0) then
        dhan1rat = dcmplx(besj1,besy1r+besy1s)
       else
        dhan1rat = dcmplx(-besj1,besy1r+besy1s)
       end if


      else
       z = 8.d0/ax
       y = z*z
       xx = ax - thrpibyf
       pofy =  (r0 + y*(r1 + y*(r2 + y*(r3 + y*r4)))) /
     *         (s0 + y*(s1 + y*(s2 + y*(s3 + y*s4))))
       qofy =  z*(u0 + y*(u1 + y*(u2 + y*(u3 + y*u4)))) /
     *           (v0 + y*(v1 + y*(v2 + y*(v3 + y*v4))))  

       cosxx = cos(xx)     
       sinxx = sin(xx)
       sqrtz = sqrt(obyfopi*z)
       besj1 = sqrtz * (pofy*cosxx - qofy*sinxx)
       besy1 = sqrtz * (pofy*sinxx + qofy*cosxx)

       if(x.gt.0.d0) then
        dhan1rat = dcmplx(besj1,besy1)
       else
        dhan1rat = dcmplx(besj1,-besy1)
       end if

      end if

      return
      end



C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     BESSJ
C     Bessel function Jn, n=2,3,.., of a REAL*8 argument
C     W.H. Press et.al: 'Numerical Recipes', 
C     Cambridge University Press, 1989, p. 175-176.
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      real*8 function bessj(n,x)
      implicit real*8 (a-h,o-z)

      parameter (iacc=225, bign0=1.0D10, bigni=1.0D-10)
C     nr of significant figures of accuracy ~ sqrt(iacc) 

      if (n.lt.2) stop 'bessj, n too small'
      ax = dabs(x)
      if (ax.eq.0.0d0) then
        bessj = 0.0d0
      else if ( sngl(ax) .gt. float(n) ) then
C       upwards recurrence
        tox = 2/ax
        bjm = dbes0rat(ax)
        bj = dbes1rat(ax)
        do j = 1, n-1
          bjp = j*tox*bj - bjm
          bjm = bj
          bj = bjp
        end do
        bessj = bj
      else
C       downwards recurrence
        tox = 2/ax
        m = 2*((n+int(sqrt(float(iacc*n))))/2)
        bessj = 0.0d0
        jsum = 0
        sum = 0.0d0
        bjp = 0.0d0
        bj = 1.0d0
        do j = m, 1, -1
          bjm = j*tox*bj - bjp
          bjp = bj
          bj = bjm
          if ( dabs(bj) .gt. bign0 ) then
            bj = bj*bigni
            bjp = bjp*bigni
            bessj = bessj*bigni
            sum = sum*bigni
          end if
          if (jsum.ne.0) sum = sum+bj
          jsum = 1-jsum
          if (j.eq.n) bessj = bjp
        end do
        sum = 2*sum - bj
        bessj = bessj/sum
      end if
      if ( x.lt.0.0d0 .and. mod(n,2).eq.1 ) bessj = -bessj

      return
      end



C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     BESSH
C     Hankel function Hn(1), n=2,3,.., of a REAL*8 argument
C     W.H. Press et.al: 'Numerical Recipes', 
C     Cambridge University Press, 1989, p. 174.
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      complex*16 function bessh(n,x)
      implicit real*8 (a-h,o-z)
      complex*16 dhan0rat,dhan1rat

      if (n.lt.2) stop 'bessh, n too small'
      if (x.eq.0.0d0) then
        bessh = dcmplx(0.0d0,0.0d0)   ! singular point
        return
      end if

      ax = dabs(x)
      tox = 2/ax
      by = dimag(dhan1rat(ax))
      bym = dimag(dhan0rat(ax))
      do j = 1, n-1
        byp = j*tox*by - bym
        bym = by
        by = byp
      end do
      bessy = by

      bessr = bessj(n,x)

      if (x.gt.0.0d0) then
        bessh = dcmplx(bessr,bessy)
      else if ( n .eq. 2*(n/2) ) then
        bessh = dcmplx(-bessr,bessy)
      else
        bessh = dcmplx(-bessr,-bessy)
      end if

      return
      end




C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     BESSJD
C     derivative Jn'(x) of Bessel function Jn, n=0,1,..
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      real*8 function bessjd(n,x)
      implicit real*8 (a-h,o-z)

      if (n.lt.0) stop 'bessjd, n negative'

      if (n.eq.0) then
        bessjd = -dbes1rat(x)
      else if (n.eq.1) then
        bessjd = 0.5d0 * ( dbes0rat(x) - bessj(2,x) )
      else if (n.eq.2) then
        bessjd = 0.5d0 * ( dbes1rat(x) - bessj(3,x) )
      else
        bessjd = 0.5d0 * ( bessj(n-1,x) - bessj(n+1,x) )
      end if

      return
      end


C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     BESSHD
C     derivative Hn(1)'(x) of Hankel function Hn(1), n=0,1,..
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      complex*16 function besshd(n,x)
      implicit real*8 (a-h,o-z)
      complex*16 dhan0rat,dhan1rat,bessh

      if (n.lt.0) stop 'besshd, n negative'

      if (n.eq.0) then
        besshd = -dhan1rat(x)
      else if (n.eq.1) then
        besshd = 0.5d0 * ( dhan0rat(x) - bessh(2,x) )
      else if (n.eq.2) then
        besshd = 0.5d0 * ( dhan1rat(x) - bessh(3,x) )
      else
        besshd = 0.5d0 * ( bessh(n-1,x) - bessh(n+1,x) )
      end if

      return
      end

