C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     Bessel and Hankel function routines, from Anders Bostrom
C     891108, a few minor changes by Ilkka Karasalo (IK)
C     also, a few minor changes by Sven Ivansson (SI)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C THE FILE CONTAINS THE MAIN SUBROUTINES 
C 
C CDHAN
C CDBES
C DHAN
C DBES
C
C WHICH COMPUTE HANKEL AND BESSEL FUNCTIONS
C
C CSUMNEUSI (adapted from csumneu) HAS BEEN ADDED BY SVEN IVANSSON
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     891212:  The functions CDBES and CDHAN now return zero
C              if ABS(Im(Z)) > ZIMMAX = 200.  (This number was
C              found by numerical experiment.)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CDHAN
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      SUBROUTINE CDHAN(MMAX1,Z,B,H)
C
C COMPUTES THE BESSEL AND HANKEL FUNTIONS OF ORDER 0,1,...,MMAX1-1
C AND COMPLEX ARGUMENT Z AND PUTS THEM IN THE ARRAYS B(1),... AND H(1),...
C MMAX1<52; AT LEAST 9 FIGURES ARE ACCURATE IN H AND 13 IN B
C USES SUBROUTINES CDBES, CSUMBES, CASBES, CSUMNEU, AND CASHAN
C
      implicit real*8 (a-h,o-z)
      complex*16 z,b(*),h(*),j0,j1,y0,y1,y2,h0,h1,
     *  twobyz


      parameter ( zimmax = 200.d0) ! Max ABS(Im(Z)) allowed
c                                  ! See comment first in file
      logical*2 ovfl

      zab = abs(z)
      y=dimag(z)
      ovfl = abs(y).gt.zimmax

      if(ovfl.or.zab.eq.0.d0) then
c     Z = 0: Return H = 0  (H has pole at 0)
      b(1) = dcmplx(1.d0,0.d0)
      if(ovfl) b(1) = dcmplx(0.d0,0.d0)
	IF (ZAB.EQ.0.0D0) THEN
C	  "overflow" but choose symmetry in  " J0 = 0.5d0 * ( H0(1) + H0(2) ) "
	  H(1) = (1.0D0,0.0D0)
	ELSE
      h(1) = dcmplx(0.d0,0.d0)
	END IF                                   !!! SI
      do 1 i = 2, mmax1
      b(i) = dcmplx(0.d0,0.d0)
 1    h(i) = dcmplx(0.d0,0.d0)
      return
      end if

     

      x=dreal(z)
      twobyz = 2.0d0/z

      call cdbes(mmax1,z,b)

      if((x/12)**2+(y/8)**2.lt.1.d0) then
      j0=b(1)
c     If MMAX1 = 1, use a dummy value for J1. Note that
c     the Y0 value returned by CSUMNEU does not depend on
c     the input J1, and is correct if J0 is correct
      if(mmax1.ge.2) then
      j1=b(2)
      else
      j1 = dcmplx(0.d0,0.d0)
      end if
 
      call csumneu(z,j0,j1,y0,y1)

      h(1) = j0 + dcmplx(-dimag(y0),dreal(y0))
      if(mmax1.ge.2)
     * h(2) = j1 + dcmplx(-dimag(y1),dreal(y1))
      do 10 m=3,mmax1
      y2 = twobyz*(m-2)*y1 - y0
      y0 = y1
      y1 = y2
   10 h(m) = b(m) + dcmplx(-dimag(y1),dreal(y1))

      else
      call cashan(z,h0,h1)
      h(1)=h0
      if(mmax1.ge.2) h(2)=h1
      do 30 m=3,mmax1
   30 h(m)=twobyz*(m-2)*h(m-1)-h(m-2)
      end if

 999  return
      end





C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CDHANSVEN     as cdhan but without "B" (only "H")
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      SUBROUTINE CDHANSVEN(MMAX1,Z,H)
C
C COMPUTES THE HANKEL FUNTIONS OF ORDER 0,1,...,MMAX1-1
C AND COMPLEX ARGUMENT Z AND PUTS THEM IN THE ARRAY H(1),...
C MMAX1<52; AT LEAST 9 FIGURES ARE ACCURATE IN H
C USES SUBROUTINES CDBES, CSUMBES, CASBES, CSUMNEU, AND CASHAN
C
      implicit real*8 (a-h,o-z)
      complex*16 z,b(11),h(*),j0,j1,y0,y1,y2,h0,h1,
     *  twobyz


      parameter ( zimmax = 200.d0) ! Max ABS(Im(Z)) allowed
c                                  ! See comment first in file
      logical*2 ovfl

	if (mmax1.gt.11) stop '%%%cdhansven b dimension'

      zab = abs(z)
      y=dimag(z)
      ovfl = y.lt.-zimmax

      if(ovfl.or.zab.eq.0.d0) then
c       Z = 0: Return H = 0  (H has pole at 0)
	IF (ZAB.EQ.0.0D0) THEN
C	  "overflow" but choose symmetry in  " J0 = 0.5d0 * ( H0(1) + H0(2) ) "
	  H(1) = (1.0D0,0.0D0)
	ELSE
          h(1) = dcmplx(0.d0,0.d0)
	END IF                                   !!! SI
        do 1 i = 2, mmax1
 1      h(i) = dcmplx(0.d0,0.d0)
        return
      end if

      x=dreal(z)
      twobyz = 2.0d0/z

      if((x/12)**2+(y/8)**2.lt.1.d0) then

        call cdbes(mmax1,z,b)
        j0=b(1)
c       If MMAX1 = 1, use a dummy value for J1. Note that
c       the Y0 value returned by CSUMNEU does not depend on
c       the input J1, and is correct if J0 is correct
        if(mmax1.ge.2) then
          j1=b(2)
        else
          j1 = dcmplx(0.d0,0.d0)
        end if
 
        call csumneu(z,j0,j1,y0,y1)

        h(1) = j0 + dcmplx(-dimag(y0),dreal(y0))
        if(mmax1.ge.2)
     *   h(2) = j1 + dcmplx(-dimag(y1),dreal(y1))
        do 10 m=3,mmax1
        y2 = twobyz*(m-2)*y1 - y0
        y0 = y1
        y1 = y2
   10   h(m) = b(m) + dcmplx(-dimag(y1),dreal(y1))

      else

        call cashan(z,h0,h1)
        h(1)=h0
        if(mmax1.ge.2) h(2)=h1
        do 30 m=3,mmax1
   30   h(m)=twobyz*(m-2)*h(m-1)-h(m-2)

      end if

 999  return
      end



C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CDHANSPEC     as cdhansven but with "common /inspec/" active
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      SUBROUTINE CDHANSPEC(MMAX1,Z,H)
C
C COMPUTES THE "inspec scaled" HANKEL FUNTIONS OF ORDER 0,1,...,MMAX1-1
C AND COMPLEX ARGUMENT Z AND PUTS THEM IN THE ARRAY H(1),...
C MMAX1<52; AT LEAST 9 FIGURES ARE ACCURATE IN H
C USES SUBROUTINES CDBES, CSUMBES, CASBES, CSUMNEU, AND CASHAN
C
      implicit real*8 (a-h,o-z)
      complex*16 z,b(11),h(*),j0,j1,y0,y1,y2,h0,h1,
     *  twobyz

	INCLUDE '../cpm/INSPEC.INC'

      parameter ( zimmax = 200.d0) ! Max ABS(Im(Z)) allowed
c                                  ! See comment first in file
      logical*2 ovfl

	if (mmax1.gt.11) stop '%%%cdhanspec b dimension'

      zab = abs(z)
      y=dimag(z)
      ovfl = y.lt.-zimmax

      if(ovfl.or.zab.eq.0.d0) then
c       Z = 0: Return H = 0  (H has pole at 0)
	IF (ZAB.EQ.0.0D0) THEN
C	  "overflow" but choose symmetry in  " J0 = 0.5d0 * ( H0(1) + H0(2) ) "
	  IF (IRKSPEC.EQ.0) THEN
	    H(1) = (1.0D0,0.0D0)
	  ELSE
	    H(1) = DCMPLX(DEXP(DLNSPEC),0.0D0)
	  END IF
	ELSE
          h(1) = dcmplx(0.d0,0.d0)
	END IF                                   !!! SI
        do 1 i = 2, mmax1
 1      h(i) = dcmplx(0.d0,0.d0)
        return
      end if

      x=dreal(z)
      twobyz = 2.0d0/z

      if((x/12)**2+(y/8)**2.lt.1.d0) then

        call cdbes(mmax1,z,b)
        j0=b(1)
c       If MMAX1 = 1, use a dummy value for J1. Note that
c       the Y0 value returned by CSUMNEU does not depend on
c       the input J1, and is correct if J0 is correct
        if(mmax1.ge.2) then
          j1=b(2)
        else
          j1 = dcmplx(0.d0,0.d0)
        end if
 
        call csumneu(z,j0,j1,y0,y1)

        h(1) = j0 + dcmplx(-dimag(y0),dreal(y0))
        if(mmax1.ge.2)
     *   h(2) = j1 + dcmplx(-dimag(y1),dreal(y1))
        do 10 m=3,mmax1
        y2 = twobyz*(m-2)*y1 - y0
        y0 = y1
        y1 = y2
   10   h(m) = b(m) + dcmplx(-dimag(y1),dreal(y1))
	IF (IRKSPEC.NE.0) THEN
	  SLASK = DEXP(DLNSPEC)
	  DO M=1,MMAX1
	    H(M) = SLASK*H(M)
	  END DO
	END IF

      else

        call cashanspec(z,h0,h1)
        h(1)=h0
        if(mmax1.ge.2) h(2)=h1
        do 30 m=3,mmax1
   30   h(m)=twobyz*(m-2)*h(m-1)-h(m-2)

      end if

 999  return
      end





C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CASHAN
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine cashan(z,h0,h1)
C
C COMPUTES H0 AND H1 USING THE ASYMPTOTIC FORMULA
C
      implicit real*8 (a-h,o-z)
      dimension ap0(13),ap1(13),aq0(13),aq1(13)
      complex*16 z,h0,h1,p0,p1,q0,q1,z2,a, onebyz, wor

      parameter ( twobypi  =  .6366197723675814D+00 ) ! 2/PI
      parameter ( pibyfour =  .7853981633974483D+00 ) ! PI/4
      parameter ( thrpibyf =  .2356194490192345D+01 ) ! 3*PI/4

      DATA AP0 /.1D1,-.703125D-1,.112152099609375D0,
     ,-.5725014209747314D0,.6074042001273483D1,
     ,-.1100171402692467D3,.3038090510922384D4,
     ,-.1188384262567832D6,.6252951493434794D7,
     ,-.4259392165047667D9,.3646840080706554D11,
     ,-.3833534661393942D13,.4854014686852896D15/

      DATA AP1 /.1D1,.1171875D0,-.144195556640625D0,
     ,.6765925884246826D0,-.6883914268109948D1,
     ,.1215978918765359D3,-.3302272294480852D4,
     ,.1276412726461746D6,-.6656367718817685D7,
     ,.4502786003050390D9,-.3833857520742787D11,
     ,.4011838599133194D13,-.5060568503314721D15/

      DATA AQ0 /-.125D0,.732421875D-1,-.2271080017089844D0,
     ,.1727727502584457D1,-.2438052969955606D2,
     ,.5513358961220205D3,-.1825775547429317D5,
     ,.8328593040162890D6,-.5006958953198890D8,
     ,.3836255180230431D10,-.3649010818849831D12,
     ,.4218971570284093D14,-.5827244631566901D16/

      DATA AQ1 /.375D0,-.1025390625D0,.2775764465332031D0,
     ,-.1993531733751297D1,.2724882731126854D2,
     ,-.6038440767050701D3,.1971837591223663D5,
     ,-.8902978767070675D6,.5310411010968520D8,
     ,-.4043620325107752D10,.3827011346598603D12,
     ,-.4406481417852275D14,.6065091351222692D16/

      IF (DIMAG(Z).LT.0.0D0) THEN
        IF (DREAL(Z).LT.0.0D0) THEN
          WRITE(*,*) '%%% besbod.f, better accuracy by symmetry; Sven I stop'
          write(*,*) 'z=',z
          STOP 5760
        END IF
      END IF

      zab = DMAX1 ( DABS(DREAL(Z)), DABS(DIMAG(Z)) ) ! SI from "cdabs(z)"  !!!
      onebyz = 1.0d0/z

      mz=zab + 1
      mm=min(mz,12)
      p0=ap0(mm+1)
      p1=ap1(mm+1)
      q0=aq0(mm+1)
      q1=aq1(mm+1)
      z2 = onebyz*onebyz
      do 10 i=1,mm
      k=mm+1-i
      p0=p0*z2+ap0(k)
      p1=p1*z2+ap1(k)
      q0=q0*z2+aq0(k)
   10 q1=q1*z2+aq1(k)
      a=cdsqrt(twobypi*onebyz)
      IF (DIMAG(Z).EQ.0.0D0) THEN
        IF (DREAL(Z).LT.0.0D0) A = -A  ! correction by Sven Ivansson dec. 1990
      END IF

c     Save four complex mults, by evaluating I* explicitly
      wor = q0*onebyz
      wor = dcmplx(-dimag(wor),dreal(wor))
      h0 = a*(p0+wor)
      wor = z - pibyfour
      wor = dcmplx(-dimag(wor),dreal(wor))
      h0 = h0*cdexp(wor)

      wor = q1*onebyz
      wor = dcmplx(-dimag(wor),dreal(wor))
      h1 = a*(p1+wor)
      wor = z - thrpibyf
      wor = dcmplx(-dimag(wor),dreal(wor))
      h1 = h1*cdexp(wor)

      return
      end



C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CASHANSPEC     as cashan but with "common /inspec/" active
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine cashanspec(z,h0,h1)
C
C COMPUTES "inspec scaled" H0 AND H1 USING THE ASYMPTOTIC FORMULA
C
      implicit real*8 (a-h,o-z)
      dimension ap0(13),ap1(13),aq0(13),aq1(13)
      complex*16 z,h0,h1,p0,p1,q0,q1,z2,a, onebyz, wor

	INCLUDE '../cpm/INSPEC.INC'

      parameter ( twobypi  =  .6366197723675814D+00 ) ! 2/PI
      parameter ( pibyfour =  .7853981633974483D+00 ) ! PI/4
      parameter ( thrpibyf =  .2356194490192345D+01 ) ! 3*PI/4

      DATA AP0 /.1D1,-.703125D-1,.112152099609375D0,
     ,-.5725014209747314D0,.6074042001273483D1,
     ,-.1100171402692467D3,.3038090510922384D4,
     ,-.1188384262567832D6,.6252951493434794D7,
     ,-.4259392165047667D9,.3646840080706554D11,
     ,-.3833534661393942D13,.4854014686852896D15/

      DATA AP1 /.1D1,.1171875D0,-.144195556640625D0,
     ,.6765925884246826D0,-.6883914268109948D1,
     ,.1215978918765359D3,-.3302272294480852D4,
     ,.1276412726461746D6,-.6656367718817685D7,
     ,.4502786003050390D9,-.3833857520742787D11,
     ,.4011838599133194D13,-.5060568503314721D15/

      DATA AQ0 /-.125D0,.732421875D-1,-.2271080017089844D0,
     ,.1727727502584457D1,-.2438052969955606D2,
     ,.5513358961220205D3,-.1825775547429317D5,
     ,.8328593040162890D6,-.5006958953198890D8,
     ,.3836255180230431D10,-.3649010818849831D12,
     ,.4218971570284093D14,-.5827244631566901D16/

      DATA AQ1 /.375D0,-.1025390625D0,.2775764465332031D0,
     ,-.1993531733751297D1,.2724882731126854D2,
     ,-.6038440767050701D3,.1971837591223663D5,
     ,-.8902978767070675D6,.5310411010968520D8,
     ,-.4043620325107752D10,.3827011346598603D12,
     ,-.4406481417852275D14,.6065091351222692D16/

      IF (DIMAG(Z).LT.0.0D0) THEN
        IF (DREAL(Z).LT.0.0D0) THEN
          WRITE(*,*) '%%% besbod.f, better accuracy by symmetry; Sven I stop'
          write(*,*) 'z=',z
          STOP 5760
        END IF
      END IF

      zab = DMAX1 ( DABS(DREAL(Z)), DABS(DIMAG(Z)) ) ! SI from "cdabs(z)"  !!!
      onebyz = 1.0d0/z

      mz=zab + 1
      mm=min(mz,12)
      p0=ap0(mm+1)
      p1=ap1(mm+1)
      q0=aq0(mm+1)
      q1=aq1(mm+1)
      z2 = onebyz*onebyz
      do 10 i=1,mm
      k=mm+1-i
      p0=p0*z2+ap0(k)
      p1=p1*z2+ap1(k)
      q0=q0*z2+aq0(k)
   10 q1=q1*z2+aq1(k)
      a=cdsqrt(twobypi*onebyz)
      IF (DIMAG(Z).EQ.0.0D0) THEN
        IF (DREAL(Z).LT.0.0D0) A = -A  ! correction by Sven Ivansson dec. 1990
      END IF

c     Save four complex mults, by evaluating I* explicitly
      wor = q0*onebyz
      wor = dcmplx(-dimag(wor),dreal(wor))
      h0 = a*(p0+wor)
      wor = z - pibyfour
      wor = dcmplx(-dimag(wor),dreal(wor))
	IF (IRKSPEC.EQ.0) THEN
      h0 = h0*cdexp(wor)
	ELSE
      h0 = h0*cdexp(wor+DLNSPEC)
	END IF

      wor = q1*onebyz
      wor = dcmplx(-dimag(wor),dreal(wor))
      h1 = a*(p1+wor)
      wor = z - thrpibyf
      wor = dcmplx(-dimag(wor),dreal(wor))
	IF (IRKSPEC.EQ.0) THEN
      h1 = h1*cdexp(wor)
	ELSE
      h1 = h1*cdexp(wor+DLNSPEC)
	END IF

      return
      end


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CSUMNEU
C     Computes Y0 and Y1 given J0 and J1, complex argument.
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine csumneu(z,j0,j1,y0,y1)

      implicit real*8 (a-h,o-z)
      complex*16 z,j0,j1,y0,y1,z0,z2,f0,t0,f1,t1,a, onebyz

      parameter ( twobypi  =  .6366197723675814D+00 ) ! 2/PI
      parameter ( onebypi  =  .3183098861837907D+00 ) ! 1/PI

      onebyz = 1.0d0/z

      p0 = -.5772156649015329d0 
      p1 = p0 + 1.d0
      z0 = 0.5d0*z
      a=cdlog(z0)
      y0=twobypi*(a*j0-p0)
      y1=(2*(a*j1-onebyz)-z0*(p0+p1))*onebypi

      z2=-z0*z0
      f0=twobypi*z2
      f1=0.5d0*onebypi*z0*z2
      onebyi = 0.5d0
      do 10 i=2,100
      onebyi1 = 1.d0/(i+1)
      p0=p1
      p1=p1 + onebyi
      t0=f0*p0
      t1=f1*(p0+p1)
      y0=y0-t0
      y1=y1-t1

      IF ( DMAX1(DABS(DREAL(T0)),DABS(DIMAG(T0))) .LT. 
     $1.D-13*DMAX1(DABS(DREAL(Y0)),DABS(DIMAG(Y0))) ) GOTO 20
C     SI change from  " if(cdabs(t0).lt.1.d-13*cdabs(y0)) goto 20 "    !!!

      f0 = f0*z2*onebyi*onebyi
      f1 = f1*z2*onebyi*onebyi1
 10   onebyi = onebyi1

      stop 451

   20 continue
      return
      end





C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CSUMNEUSI        <adapted from csumneu by sven ivansson>
C     Computes "Y0 and Y1 except for the logarithmic terms and the 
C     1/z term", returned in y0reg and y1reg which are even entire 
C     functions whose power-series have real coefficients, complex 
C     argument.
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine csumneusi(z,y0reg,y1reg)

      implicit real*8 (a-h,o-z)
      complex*16 z,y0reg,y1reg,z0,z2,f0,t0,f1,t1

      parameter ( twobypi  =  .6366197723675814D+00 ) ! 2/PI
      parameter ( onebypi  =  .3183098861837907D+00 ) ! 1/PI

      p0 = -.5772156649015329d0 
      p1 = p0 + 1.d0
      z0 = 0.5d0*z
      y0reg=twobypi*(-p0)
      y1reg=(-z0*(p0+p1))*onebypi

      z2=-z0*z0
      f0=twobypi*z2
      f1=0.5d0*onebypi*z0*z2
      onebyi = 0.5d0
      do 10 i=2,100
      onebyi1 = 1.d0/(i+1)
      p0=p1
      p1=p1 + onebyi
      t0=f0*p0
      t1=f1*(p0+p1)
      y0reg=y0reg-t0
      y1reg=y1reg-t1
      IF ( DMAX1(DABS(DREAL(T0)),DABS(DIMAG(T0))) .LT. 
     $1.D-13*DMAX1(DABS(DREAL(Y0REG)),DABS(DIMAG(Y0REG))) ) GOTO 20
C     SI change from  " if(cdabs(t0).lt.1.d-13*cdabs(y0reg)) goto 20 "    !!!
      f0 = f0*z2*onebyi*onebyi
      f1 = f1*z2*onebyi*onebyi1
 10   onebyi = onebyi1

      stop 451

   20 continue
      return
      end


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CDBES
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine cdbes(mmax1,za,b)
C
C COMPUTES THE BESSEL FUNTIONS OF ORDER 0,1,...,MMAX1-1
C AND COMPLEX ARGUMENT Z AND PUTS THEM IN THE ARRAY B(1),...
C MMAX1<52; ABOUT 13 FIGURES ARE ACCURATE
C USES SUBROUTINES CSUMBES CASBES
C
c     891209 IK:  Argument ZA is now input parameter only.
c                 Note that Z is changed below.

      implicit real*8 (a-h,o-z)
      complex*16 za,z,b(*),j0,j1,j2, twobyz


      parameter ( zimmax = 200.d0) ! Max ABS(Im(Z)) allowed
c                                  ! See comment first in file
      logical*2 ovfl, negaz

      z = za
      x = DMAX1 ( DABS(DREAL(Z)), DABS(DIMAG(Z)) ) ! SI from "cdabs(z)"  !!!
      mx=x
      ovfl = abs(dimag(z)).gt.zimmax

      if(x.eq.0.d0.or.ovfl) then
      b(1) = dcmplx(1.d0,0.d0)
      if(ovfl) b(1) = dcmplx(0.d0,0.d0)
      do 1 i = 2, mmax1
 1    b(i) = dcmplx(0.d0,0.d0)
      return
      end if

      twobyz = 2.0d0/z

      if(mx.gt.mmax1.and.mx.ge.20) then
      negaz = dreal(z).lt.0.d0
      if(negaz) z = -z
      call casbes(z,j0,j1)
      if(negaz) z = -z
      if(negaz) j1 = -j1
      b(1)=j0

c     891209 IK: Do not set B(2) when MMAX1 = 1 !!
      if(mmax1.ge.2) b(2)=j1
      do 10 m = 3,mmax1
 10   b(m) = (m-2)*twobyz*b(m-1) - b(m-2)
      return

      else
      mx=2*x-10
      mm=max(mmax1,mx,3)
      call csumbes(mm-1,z,j2)   ! J2 = Jmm-1(Z)
      call csumbes(mm-2,z,j1)   ! J1 = Jmm-2(Z)

      if(mm.eq.mmax1) b(mm) = j2
      if(mm-1.le.mmax1) b(mm-1) = j1

      do 20 i=3,mm
      m=mm-i+1
      j0 = twobyz*m*j1 - j2
      j2 = j1
      j1 = j0
   20 if(m.le.mmax1) b(m) = j0 

      end if

      return
      end


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CSUMBES
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      SUBROUTINE CSUMBES(M,Z,JM)
C
C COMPUTES JM(Z) FOR COMPLEX Z USING THE SUM FORMULA 
C
      COMPLEX*16 Z,JM,TERM,Z2
      JM=1.D0
      TERM=1.D0
      Z2=-.25D0*Z*Z

      do 10 i=1,100
      kden = i*(i+m)             
      term=term*z2
      term=term/kden
      jm = jm + term
   10 IF ( DMAX1(DABS(DREAL(TERM)),DABS(DIMAG(TERM))) .LT. 
     $1.D-13 ) GOTO 20
C     SI change from  " 10 if(cdabs(term).lt.1.d-13) goto 20 "    !!!

   20 if(m.eq.0) return

      z2 = 0.5d0*z
      do 30 j=1,m
      jm = jm*z2
   30 jm = dcmplx(dreal(jm)/j,dimag(jm)/j)

      return
      end


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CASBES
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      SUBROUTINE CASBES(Z,J0,J1)
C
C COMPUTES J0 AND J1 FOR COMPLEX Z USING THE ASYMPTOTIC FORMULA
C

      implicit real*8 (a-h,o-z)
      dimension  ap0(13),ap1(13),aq0(13),aq1(13)

      complex*16 z,j0,j1,p0,p1,q0,q1,z2,y,s,c,a,onebyz

      parameter ( twobypi  =  .6366197723675814D+00 ) ! 2/PI
      parameter ( pibyfour =  .7853981633974483D+00 ) ! PI/4

      DATA AP0 /.1D1,-.703125D-1,.112152099609375D0,
     ,-.5725014209747314D0,.6074042001273483D1,
     ,-.1100171402692467D3,.3038090510922384D4,
     ,-.1188384262567832D6,.6252951493434794D7,
     ,-.4259392165047667D9,.3646840080706554D11,
     ,-.3833534661393942D13,.4854014686852896D15/
      DATA AP1 /.1D1,.1171875D0,-.144195556640625D0,
     ,.6765925884246826D0,-.6883914268109948D1,
     ,.1215978918765359D3,-.3302272294480852D4,
     ,.1276412726461746D6,-.6656367718817685D7,
     ,.4502786003050390D9,-.3833857520742787D11,
     ,.4011838599133194D13,-.5060568503314721D15/
      DATA AQ0 /-.125D0,.732421875D-1,-.2271080017089844D0,
     ,.1727727502584457D1,-.2438052969955606D2,
     ,.5513358961220205D3,-.1825775547429317D5,
     ,.8328593040162890D6,-.5006958953198890D8,
     ,.3836255180230431D10,-.3649010818849831D12,
     ,.4218971570284093D14,-.5827244631566901D16/
      DATA AQ1 /.375D0,-.1025390625D0,.2775764465332031D0,
     ,-.1993531733751297D1,.2724882731126854D2,
     ,-.6038440767050701D3,.1971837591223663D5,
     ,-.8902978767070675D6,.5310411010968520D8,
     ,-.4043620325107752D10,.3827011346598603D12,
     ,-.4406481417852275D14,.6065091351222692D16/
      DATA PI /3.141592653589793D0/

      IF (DREAL(Z).LT.0.0D0) THEN
        WRITE(*,*) '%%% besbod.f, more accurate for "-z", J0 even; Sven I stop'
        STOP 5760
      END IF

      P0=AP0(13)
      P1=AP1(13)
      Q0=AQ0(13)
      Q1=AQ1(13)

      onebyz = 1.0d0/z
      z2 = onebyz*onebyz

      DO 10 I=1,12
      K=13-I
      P0=P0*Z2+AP0(K)
      P1=P1*Z2+AP1(K)
      Q0=Q0*Z2+AQ0(K)
   10 Q1=Q1*Z2+AQ1(K)
      Y=Z-pibyfour
      S=CDSIN(Y)
      C=CDCOS(Y)
      a = cdsqrt(twobypi*onebyz)
      j0=a*(p0*c-q0*s*onebyz)
      j1=a*(p1*s+q1*c*onebyz)
      return
      end


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     DHAN
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      SUBROUTINE DHAN(MMAX1,X,B,H)
C
C COMPUTES THE BESSEL AND HANKEL FUNTIONS OF ORDER 0,1,...,MMAX1-1
C AND REAL ARGUMENT X>0 AND PUTS THEM IN THE ARRAYS B(1),... AND H(1),...
C MMAX1<52; ABOUT 12 FIGURES ARE ACCURATE
C USES SUBROUTINES SUMBES, SUMNEU, AND ASBN
C
      implicit real*8 (a-h,o-z)
      real*8 j0, j1, j2, b(*)
      complex*16 h(*)

      mx=x
      if(x.ne.0.d0) twobyx = 2.d0/x

      if(mx.gt.mmax1.and.mx.ge.20) then
      call asbn(x,j0,j1,y0,y1)
      b(1)=j0
      if(mmax1.gt.1) b(2)=j1
      do 10 m=3,mmax1
   10 b(m)=twobyx*(m-2)*b(m-1) - b(m-2)

      else

c     891209 IK: Avoid division by zero
      if(x.ne.0.d0) then
c     X is nonzero:
c     Use Taylor's formula and backward recursion w.r.t. order
      mx=2*x-10
      mm=max(mmax1,mx,3)
      call sumbes(mm-1,x,j2)
      call sumbes(mm-2,x,j1)
      
      if(mm.le.mmax1) b(mm) = j2
      if(mm-1.le.mmax1) b(mm-1) = j1

      do 20 i=3,mm
      m=mm-i+1
      j0 = twobyx*m*j1 - j2
      j2 = j1
      j1 = j0
   20 if(m.le.mmax1) b(m) = j0

      mx=x
      if(mx.lt.12) then
c     J0 and J1 are input data to SUMNEU
      call sumneu(x,j0,j1,y0,y1)
      else
      call asbn(x,j0,j1,y0,y1)
      endif

      else   
c     X is zero 
c     The bessel fcns of the second kind Y0, Y1,... and
c     the Hankel fcns have a pole at X = 0, and are thus
c     udefined. Their value is set to zero.
      b(1) = 1.d0
C	"overflow" but choose symmetry in  " J0 = 0.5d0 * ( H0(1) + H0(2) ) "
	H(1) = (1.0D0,0.0D0)                   !!! SI
      do 30 m = 3, mmax1
      h(m) = dcmplx(0.d0,0.d0)
 30   b(m) = 0.d0
      return
      end if   ! End X = 0

      end if  ! End use Taylor's formel and backw recursion


c     Here B(k) = Jk-1(X), k = 1,...,MMAX1 and
c     Y0 = Y0(X), Y1 = Y1(X)
      h(1) = dcmplx(b(1),y0)
      if(mmax1.ge.2) h(2) = dcmplx(b(2),y1)

c     When is the forward recursion below stable?
      do 40 m=3,mmax1
      y2 = twobyx*(m-2)*y1 - y0
      y0 = y1
      y1 = y2
   40 h(m) = dcmplx(b(m),y1)

      return
      end


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     SUMBES
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      SUBROUTINE SUMBES(M,X,JM)
C
C COMPUTES JM(X) USING THE SUM FORMULA
C
      implicit real*8 (a-h,o-z)
      real*8 jm

      JM=1.D0
      TERM=1.D0
      X2=-.25D0*X*X
      DO 10 I=1,100
      TERM=TERM*X2/((M+I)*I)
      JM=JM+TERM
   10 IF(DABS(TERM).LT.1.D-13) GOTO 20
   20 IF(M.EQ.0) RETURN
      DO 30 J=1,M
   30 JM=JM*X/(2*(M+1-J))
      RETURN
      END


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     ASBN
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      SUBROUTINE ASBN(X,J0,J1,Y0,Y1)
C
C COMPUTES J0,J1,Y0, AND Y1 USING THE ASYMPTOTIC FORMULA
C
      implicit real*8 (a-h,o-z)
      real*8 j0, j1
      dimension ap0(13),ap1(13),aq0(13),aq1(13)
      parameter ( twobypi  =  .6366197723675814D+00 ) ! 2/PI

      DATA AP0 /.1D1,-.703125D-1,.112152099609375D0,
     ,-.5725014209747314D0,.6074042001273483D1,
     ,-.1100171402692467D3,.3038090510922384D4,
     ,-.1188384262567832D6,.6252951493434794D7,
     ,-.4259392165047667D9,.3646840080706554D11,
     ,-.3833534661393942D13,.4854014686852896D15/
      DATA AP1 /.1D1,.1171875D0,-.144195556640625D0,
     ,.6765925884246826D0,-.6883914268109948D1,
     ,.1215978918765359D3,-.3302272294480852D4,
     ,.1276412726461746D6,-.6656367718817685D7,
     ,.4502786003050390D9,-.3833857520742787D11,
     ,.4011838599133194D13,-.5060568503314721D15/
      DATA AQ0 /-.125D0,.732421875D-1,-.2271080017089844D0,
     ,.1727727502584457D1,-.2438052969955606D2,
     ,.5513358961220205D3,-.1825775547429317D5,
     ,.8328593040162890D6,-.5006958953198890D8,
     ,.3836255180230431D10,-.3649010818849831D12,
     ,.4218971570284093D14,-.5827244631566901D16/
      DATA AQ1 /.375D0,-.1025390625D0,.2775764465332031D0,
     ,-.1993531733751297D1,.2724882731126854D2,
     ,-.6038440767050701D3,.1971837591223663D5,
     ,-.8902978767070675D6,.5310411010968520D8,
     ,-.4043620325107752D10,.3827011346598603D12,
     ,-.4406481417852275D14,.6065091351222692D16/
      DATA PI /3.141592653589793D0/

      IF (X.LT.0.0D0) THEN
        WRITE(*,*) '%%% besbod.f, better accuracy by symmetry; Sven I stop'
        write(*,*) 'z=',z
        STOP 5760
      END IF

      P0=AP0(13)
      P1=AP1(13)
      Q0=AQ0(13)
      Q1=AQ1(13)
      X2=1/(X*X)
      DO 10 I=1,12
      K=13-I
      P0=P0*X2+AP0(K)
      P1=P1*X2+AP1(K)
      Q0=Q0*X2+AQ0(K)
   10 Q1=Q1*X2+AQ1(K)
      Y=X-PI/4
      S=DSIN(Y)
      C=DCOS(Y)
      a=dsqrt(twobypi/x)
      J0=A*(P0*C-Q0*S/X)
      J1=A*(P1*S+Q1*C/X)
      Y0=A*(P0*S+Q0*C/X)
      Y1=A*(-P1*C+Q1*S/X)
      RETURN
      END


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     SUMNEU
C     Computes Y0 and Y1 given J0 and J1
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine sumneu(x,j0,j1,y0,y1)
      implicit real*8 (a-h,o-z)
      real*8 j0, j1
      parameter ( twobypi  =  .6366197723675814D+00 ) ! 2/PI
      parameter ( onebypi  =  .3183098861837907D+00 ) ! 1/PI

      P0=-.5772156649015329D0
      P1=P0+1
      X0=X/2
      A=DLOG(X0)
      Y0 = twobypi*(A*J0-P0)
      Y1 = onebypi*(2*A*J1-1/X0-X0*(P0+P1))
      X2=-X0*X0
      F0 = twobypi*x2
      F1 = 0.5*onebypi*x0*x2
      DO 10 I=2,100
      P0=P1
      P1=P1+1.D0/I
      T0=F0*P0
      T1=F1*(P0+P1)
      Y0=Y0-T0
      Y1=Y1-T1
      IF(DABS(T0).LT.1.D-13*DABS(Y0)) GOTO 20
      F0=F0*X2/(I*I)
   10 F1=F1*X2/(I*(I+1))

      STOP 451

   20 CONTINUE
      RETURN
      END


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     DBES
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      SUBROUTINE DBES(MMAX1,X,B)
C
C COMPUTES THE BESSEL FUNTIONS OF ORDER 0,1,...,MMAX1-1
C AND REAL ARGUMENT X AND PUTS THEM IN THE ARRAY B(1),...
C MMAX1<52; ABOUT 13 FIGURES ARE ACCURATE
C USES SUBROUTINES SUMBES AND ASBES
C
      implicit real*8 (a-h,o-z)
      real*8 b(*),j0,j1,j2

      X1=DABS(X)
      MX=X1

      IF(MX.GT.MMAX1.AND.MX.GE.20) THEN
      CALL ASBES(X1,J0,J1)
      IF(X.LT.0.D0) J1=-J1
      B(1)=J0
      B(2)=J1
      IF(MMAX1.LE.2) RETURN
      DO 10 M=3,MMAX1
   10 B(M)=2*(M-2)*B(M-1)/X-B(M-2)

      ELSE

c     891207, IK:  Avoid division by zero.
      if(x.ne.0.d0) then
c     X is nonzero:
	WRITE(*,*) ' %%% ERROR IN DBES IN BESBOD HERE, NOT ',
     $	'CORRECTED %%%'  ! all components of b(1:mmax1) are not always updated
	STOP 100  ! Sven Ivansson 18/2 91, Klas Oestlund got error; USE BESRAT
      mx=2*x1-10
      mm=max(mmax1,mx,3)
      call sumbes(mm-1,x,j2)
      call sumbes(mm-2,x,j1)

      if(mm.le.mmax1) b(mm) = j2
      if(mm-1.le.mmax1) b(mm) = j1

      twobyx = 2.d0/x
      do 20 i=3,mm
      m=mm-i+1
      j0 = twobyx*m*j1 - j2
      j2 = j1
      j1 = j0
 20   if(m.le.mmax1) b(m) = j0

      else
c     X is zero:
      b(1) = 1.d0
      do 30 m = 2,mmax1
 30   b(m) = 0.d0
      end if   ! End X = 0

      ENDIF

      RETURN
      END


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     ASBES
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      SUBROUTINE ASBES(X,J0,J1)
C
C COMPUTES J0 AND J1 USING THE ASYMPTOTIC FORMULA
C
      implicit real*8 (a-h,o-z)
      real*8 j0, j1
      dimension ap0(13),ap1(13),aq0(13),aq1(13)

      parameter ( twobypi  =  .6366197723675814D+00 ) ! 2/PI
      parameter ( pibyfour =  .7853981633974483D+00 ) ! PI/4

      DATA AP0 /.1D1,-.703125D-1,.112152099609375D0,
     ,-.5725014209747314D0,.6074042001273483D1,
     ,-.1100171402692467D3,.3038090510922384D4,
     ,-.1188384262567832D6,.6252951493434794D7,
     ,-.4259392165047667D9,.3646840080706554D11,
     ,-.3833534661393942D13,.4854014686852896D15/
      DATA AP1 /.1D1,.1171875D0,-.144195556640625D0,
     ,.6765925884246826D0,-.6883914268109948D1,
     ,.1215978918765359D3,-.3302272294480852D4,
     ,.1276412726461746D6,-.6656367718817685D7,
     ,.4502786003050390D9,-.3833857520742787D11,
     ,.4011838599133194D13,-.5060568503314721D15/
      DATA AQ0 /-.125D0,.732421875D-1,-.2271080017089844D0,
     ,.1727727502584457D1,-.2438052969955606D2,
     ,.5513358961220205D3,-.1825775547429317D5,
     ,.8328593040162890D6,-.5006958953198890D8,
     ,.3836255180230431D10,-.3649010818849831D12,
     ,.4218971570284093D14,-.5827244631566901D16/
      DATA AQ1 /.375D0,-.1025390625D0,.2775764465332031D0,
     ,-.1993531733751297D1,.2724882731126854D2,
     ,-.6038440767050701D3,.1971837591223663D5,
     ,-.8902978767070675D6,.5310411010968520D8,
     ,-.4043620325107752D10,.3827011346598603D12,
     ,-.4406481417852275D14,.6065091351222692D16/

      IF (X.LT.0.0D0) THEN
        WRITE(*,*) '%%% besbod.f, better accuracy by symmetry; Sven I stop'
        write(*,*) 'z=',z
        STOP 5760
      END IF

      P0=AP0(13)
      P1=AP1(13)
      Q0=AQ0(13)
      Q1=AQ1(13)
      X2=1/(X*X)
      DO 10 I=1,12
      K=13-I
      P0=P0*X2+AP0(K)
      P1=P1*X2+AP1(K)
      Q0=Q0*X2+AQ0(K)
   10 Q1=Q1*X2+AQ1(K)
      y=x-pibyfour
      S=DSIN(Y)
      C=DCOS(Y)
      a=dsqrt(twobypi/x)
      J0=A*(P0*C-Q0*S/X)
      J1=A*(P1*S+Q1*C/X)
      RETURN
      END


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     DBES0BOD
C     Computes J0(Z) for real Z, by calling DBES
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      real*8 function dbes0bod(z)
      real*8 z,j0(2)
C-----------------------------------------------------------------
C     About 13 figures are accurate
C-----------------------------------------------------------------
      call dbes(1,z,j0)
      dbes0bod = j0(1)
      return
      end


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CBES0BOD
C     Computes J0(Z) for complex Z, by calling CDBES
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      complex*16 function cbes0bod(z)
      complex*16 z,j0(2)
C-----------------------------------------------------------------
C     About 13 figures are accurate
C-----------------------------------------------------------------
      call cdbes(1,z,j0)
      cbes0bod = j0(1)
      return
      end


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     DHAN0BOD
C     Computes H0(Z) for real Z, by calling DHAN
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      complex*16 function dhan0bod(z)
      real*8 z, j0(2)
      complex*16 h0(2)
C-----------------------------------------------------------------
C-----------------------------------------------------------------
      call dhan(1,z,j0,h0)
      dhan0bod = h0(1)
      return
      end


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CHAN0BOD
C     Computes H0(Z) for complex Z, by calling CDHANSVEN (changed from cdhan)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      complex*16 function chan0bod(z)
      complex*16 z, h0(2)
C-----------------------------------------------------------------
C-----------------------------------------------------------------
      call cdhansven(1,z,h0)   ! changed from  "call cdhan(1,z,j0,h0)"
      chan0bod = h0(1)
      return
      end



C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CHAN0SPEC     as chan0bod but with "common /inspec/" active
C     Computes "inspec scaled" H0(Z) for complex Z, by calling CDHANSPEC
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      complex*16 function chan0spec(z)
      complex*16 z, h0(2)
C-----------------------------------------------------------------
C-----------------------------------------------------------------
      call cdhanspec(1,z,h0)
      chan0spec = h0(1)
      return
      end



C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     DBES1BOD
C     Computes J1(Z) for real Z, by calling DBES
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      real*8 function dbes1bod(z)
      real*8 z,j1(2)
C-----------------------------------------------------------------
C     About 13 figures are accurate
C-----------------------------------------------------------------
      call dbes(2,z,j1)
      dbes1bod = j1(2)
      return
      end


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CBES1BOD
C     Computes J1(Z) for complex Z, by calling CDBES
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      complex*16 function cbes1bod(z)
      complex*16 z,j1(2)
C-----------------------------------------------------------------
C     About 13 figures are accurate
C-----------------------------------------------------------------
      call cdbes(2,z,j1)
      cbes1bod = j1(2)
      return
      end


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     DHAN1BOD
C     Computes H1(Z) for real Z, by calling DHAN
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      complex*16 function dhan1bod(z)
      real*8 z, j1(2)
      complex*16 h1(2)
C-----------------------------------------------------------------
C-----------------------------------------------------------------
      call dhan(2,z,j1,h1)
      dhan1bod = h1(2)
      return
      end


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CHAN1BOD
C     Computes H1(Z) for complex Z, by calling CDHANSVEN (changed from cdhan)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      complex*16 function chan1bod(z)
      complex*16 z, h1(2)
C-----------------------------------------------------------------
C-----------------------------------------------------------------
      call cdhansven(2,z,h1)
      chan1bod = h1(2)
      return
      end



C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CHAN1SPEC
C     Computes "inspec scaled" H1(Z) for complex Z, by calling CDHANSPEC
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      complex*16 function chan1spec(z)
      complex*16 z, h1(2)
C-----------------------------------------------------------------
C-----------------------------------------------------------------
      call cdhanspec(2,z,h1)
      chan1spec = h1(2)
      return
      end




C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     DBESBOD
C     Computes Jn(Z) for real Z, by calling DBES
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      real*8 function dbesbod(n,z)
      real*8 z,j1(11)
C-----------------------------------------------------------------
C     About 13 figures are accurate
C-----------------------------------------------------------------
      if ( n.lt.0 .or. n.gt.10 ) stop 'dbesbod'
      call dbes(n+1,z,j1)
      dbesbod = j1(n+1)
      return
      end


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CBESBOD
C     Computes Jn(Z) for complex Z, by calling CDBES
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      complex*16 function cbesbod(n,z)
      complex*16 z,j1(11)
C-----------------------------------------------------------------
C     About 13 figures are accurate
C-----------------------------------------------------------------
      if ( n.lt.0 .or. n.gt.10 ) stop 'cbesbod'
      call cdbes(n+1,z,j1)
      cbesbod = j1(n+1)
      return
      end


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     DHANBOD
C     Computes Hn(Z) for real Z, by calling DHAN
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      complex*16 function dhanbod(n,z)
      real*8 z, j1(11)
      complex*16 h1(11)
C-----------------------------------------------------------------
C-----------------------------------------------------------------
      if ( n.lt.0 .or. n.gt.10 ) stop 'dhanbod'
      call dhan(n+1,z,j1,h1)
      dhanbod = h1(n+1)
      return
      end


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CHANBOD
C     Computes Hn(Z) for complex Z, by calling CDHANSVEN (changed from cdhan)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      complex*16 function chanbod(n,z)
      complex*16 z, h1(11)
C-----------------------------------------------------------------
C-----------------------------------------------------------------
      if ( n.lt.0 .or. n.gt.10 ) stop 'chanbod'
      call cdhansven(n+1,z,h1)
      chanbod = h1(n+1)
      return
      end



C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CHANSPEC
C     Computes "inspec scaled" Hn(Z) for complex Z, by calling CDHANSPEC
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      complex*16 function chanspec(n,z)
      complex*16 z, h1(11)
C-----------------------------------------------------------------
C-----------------------------------------------------------------
      if ( n.lt.0 .or. n.gt.10 ) stop 'chanspec'
      call cdhanspec(n+1,z,h1)
      chanspec = h1(n+1)
      return
      end



C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CBESBODD
C     Computes Jn'(Z) for complex Z, by calling CBESBOD
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      complex*16 function cbesbodd(n,z)
      complex*16 z, cbesbod
C-----------------------------------------------------------------
C     About 13 figures are accurate
C-----------------------------------------------------------------
      if ( n.lt.0 ) stop 'cbesbodd'
      if (n.eq.0) then
        cbesbodd = -cbesbod(1,z)
      else
        cbesbodd = 0.5d0 * ( cbesbod(n-1,z) - cbesbod(n+1,z) )
      end if
      return
      end



C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CHANBODD
C     Computes Hn'(Z) for complex Z, by calling CHANBOD
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      complex*16 function chanbodd(n,z)
      complex*16 z, chanbod
C-----------------------------------------------------------------
C     About 13 figures are accurate
C-----------------------------------------------------------------
      if ( n.lt.0 ) stop 'chanbodd'
      if (n.eq.0) then
        chanbodd = -chanbod(1,z)
      else
        chanbodd = 0.5d0 * ( chanbod(n-1,z) - chanbod(n+1,z) )
      end if
      return
      end



C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CHANSPECD
C     Computes  "inspec scaled" Hn'(Z) for complex Z, by calling CHANSPEC
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      complex*16 function chanspecd(n,z)
      complex*16 z, chanspec
C-----------------------------------------------------------------
C     About 13 figures are accurate
C-----------------------------------------------------------------

      if ( n.lt.0 ) stop 'chanspecd'
      if (n.eq.0) then
        chanspecd = -chanspec(1,z)
      else
        chanspecd = 0.5d0 * ( chanspec(n-1,z) - chanspec(n+1,z) )
      end if
      return
      end

