function Csbj(l, z)
c***********************************************************************
c    Spherical Bessel function for complex z and integer l = 0, 1, 2
c***********************************************************************
c
      implicit none
      integer l
      double precision a1, a3, a5, a7, a9, a11
      double precision a2, a4, a6, a8, a10, a12
      double precision tiny, quarter, one, three, x, y
      double complex z, z2, zero, Csbj
c
      parameter (zero = (0.0d0, 0.0d0))
      parameter (tiny = 0.01d0, quarter = 0.25d0)
      parameter (one = 1.0d0, three = 3.0d0)
      parameter (a1=  0.3333333333333333333333333333d00)
      parameter (a3= -0.3333333333333333333333333333d-1)
      parameter (a5=  0.1190476190476190476190476190d-2)
      parameter (a7= -0.2204585537918871252204585538d-4)
      parameter (a9=  0.2505210838544171877505210839d-6)
      parameter (a11=-0.1927085260418593751927085260d-8)
      parameter (a2 = 0.666666666666666666666666666667d-1)
      parameter (a4 =-0.476190476190476190476190476190d-2)
      parameter (a6 = 0.132275132275132275132275132275d-3)
      parameter (a8 =-0.200416867083533750200416867084d-5)
      parameter (a10= 0.192708526041859375192708526042d-7)
      parameter (a12=-0.128472350694572916795139017361d-9)
c
      if (l .eq. 0) then
         if (z .eq. zero) then
            Csbj = one
         else
            Csbj = sin(z)/z
         endif
      else if (l .eq. 1) then
         x = dble(z)
         y = imag(z)
         z2 = z*z
         if ((x*x+y*y) .lt. tiny) then
            Csbj = z*(a1+z2*(a3+z2*(a5+z2*(a7+z2*(a9+a11*z2)))))
         else
            Csbj = sin(z)/z2-cos(z)/z
         endif
      else if (l .eq. 2) then
         z2=z*z
         x = dble(z)
         y = imag(z)
         if ((x*x+y*y) .lt. quarter) then
            Csbj=z2*(a2+z2*(a4+z2*(a6+z2*(a8+z2*(a10+a12*z2)))))
         else
            Csbj=sin(z)/z*(three/z2-one)-three*cos(z)/z2
         endif
      else
         write(8, 900) l
         stop 'Csbj: l value out of range'
      endif
c
  900 format('0', 'l value = ', i2, ', only 0, 1, 2 allowed')
c
      return
      end