      Subroutine spbesl (lmaxa,xx,fj)
c *** see RH Landau, Program lpott, 1981
c     spbesl calculates spherical bessel fuctions of 1st      kind for a
c     for all order from zero to lmaxa, at argumrnt xx
      Implicit Real*8  (a-h,o-z)
      Dimension fj(50)
      lmax = lmaxa
      r = xx
      x = xx
      lx = 2.*r
      If (49..lt.lx) lx = 49.
c     calculate small argument limit
      If (r.gt..01) GoTo 10
      fj0 = 1.-x**2/6.+x**4/120.-x**6/42./120.
      If (lmax.gt.1) GoTo 90
      fj(1) = fj0
      fj(2) = x*(1.-x**2/10.+x**4/280.-x**6/30./56./9.)/3.
      GoTo 130
c     regular value of x, start with large l,iterate Down,renormalize
 10   fj0 = sin(x)/x
 20   l2 = max0(lmax+5,lx)
      xl2 = 2.*l2
      fj(l2+1) = 1.e-10
      fj(l2) = 1.e-10*((xl2+1.)/x-x/(xl2+3.))
      l3 = l2-1
      Do 40 ll=1,l3
         l1 = l2-ll
         fl1 = l1
         fj(l1) = (2.*fl1+1.)*fj(l1+1)/x-fj(l1+2)
         rhl = abs(fj(l1))
         If (rhl.le.1.e30) GoTo 40
c     fj too large, scale Down and/or set too small fj*s to 0
         Do 30 l4=l1,l2
            If ((abs(fj(l4)).lt.1.e-26).and.(l4.ne.1)) fj(l4) = 0.
c     nb changed 50 to 26 for vax
            fj(l4) = 1.e-10*fj(l4)
 30      Continue
 40   Continue
c     applay renormaliz factor
      zzj = fj0/fj(1)
      Do 60 l1=1,2
         If (fj(l1).eq.0.) GoTo 50
         If (log10(abs(zzj))+log10(abs(fj(l1))).le.-60.) fj(l1) = 0.
 50      fj(l1) = zzj*fj(l1)
 60   Continue
      If (lmax.le.1) GoTo 130
      l2 = l2-4
      Do 80 l1=3,l2
         If (fj(l1).eq.0.) GoTo 70
         If (log10(abs(zzj))+log10(abs(fj(l1))).le.-60.) fj(l1) = 0.
 70      fj(l1) = zzj*fj(l1)
 80   Continue
      GoTo 130
 90   If (r) 120,100,20
c     x=0 values
 100  Write (6,140)
      l2 = lmax+1
      Do 110 ll=2,l2
         fj(ll) = 0.
 110  Continue
      fj(1) = fj0
      GoTo 130
 120  Write (6,150)
 130  Continue
      Return
c
 140  format (26h spbesl(x) called for x=0,)
 150  format (31h error in spbesl, neg. argument)
      End
