function SBj(l, x)
c
      implicit none
      integer l
      double precision x, x2, SBj
      double precision a1, a3, a5, a7, a9, a11
      double precision a2, a4, a6, a8, a10, a12
      double precision zero, tiny, quarter, one, three
c
      parameter (zero = 0.0d0, one = 1.0d0, three = 3.0d0)
      parameter (tiny=0.01d0, quarter = 0.25d0)
      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 (x .eq. zero) then
            SBj = one
         else
            SBj = sin(x)/x
         endif
      else if (l .eq. 1) then
         x2 = x*x
         if (x2 .lt. tiny) then
            SBj = x*(a1+x2*(a3+x2*(a5+x2*(a7+x2*(a9+a11*x2)))))
         else
            SBj = sin(x)/x2-cos(x)/x
         endif
      else if (l .eq. 2) then
         x2=x*x
         if (x2 .lt. quarter) then
            SBj=x2*(a2+x2*(a4+x2*(a6+x2*(a8+x2*(a10+a12*x2)))))
         else
            SBj=sin(x)/x*(three/x2-one)-three*cos(x)/x2
         endif
      else
         write(8, 900) l
         stop 'SBj: l value out of range'
      endif
c
  900 format('0', 'l = ', i2, ', out of range. Only 0, 1, 2 allowed')
c
      return
      end