       Subroutine beslia (zz,bi)
c *** see RH Landau, Program lpott, 1981
c     ------------------------------------------------------------------
c     calculates  exp(-z)*(bessel  functions  of  imaginary  argument)
c      z=argument  of  spherical  bessel  function
c      for  z.lt.1.   series  expansion  in  powers  of  z
c     for100.szs1.
c       for  100.>z>1.      recursion  Down   is  used   for  z>1..
c       for  100.>z>1.      recursion  Down   is  used   for  z>100.
c     recursion   up.
c     ------------------------------------------------------------------
      Implicit  Real*8   (a-h,o-z)
      Real*8 zz
      Dimension bi(100)
      Dimension r1b(100),r2b(100)
      z = zz
      If (z.gt.1.) GoTo 60
      a = exp(-z)
      Do 30 i=1,51
         b = 1.
         c = 1.
         Do 10 i1=1,30
            b = b*z**2*.5/i1/(2.*(i+i1)-1)
            c = c+b
            If (b.lt.1.e-14) GoTo 20
 10      Continue
 20      bi(i) = a*c
         If (a.le.1.e-35) GoTo 40
c     nb changed 60 to 35 for vax
         a = a*z/(2.*i+1)
 30   Continue
      Return
 40   Do 50 i2=i,51
         bi(i2) = 0.
 50   Continue
      Return
 60   If (z.lt.100.) GoTo 100
      bi(1) = 0.5/z
      bi(2) = 0.5/z*(1.-1./z)
      Do 70 i=3,71
         bi(i) = bi(i-2)-(2.*i-3)/z*bi(i-1)
 70   Continue
      Do 80 i=1,69
         j = 70-i
         bi(j) = bi(j+2)+(2.*j+1)/z*bi(j+1)
 80   Continue
      c = 0.5/z/bi(1)
      Do 90 i=1,71
         bi(i) = bi(i)*c
 90   Continue
      Return
 100  Continue
c     bi(71) = 1.e-60
c     bi(70) = 1.e-60
c     Do 110 i=1,69
c        j = 70-i
c        bi(j) = bi(j+2)+(2.*j+1)/z*bi(j+1)
c110  Continue
c     c = 0.5*(1.-exp(-2.*z))/z/bi(1)
c     Do 120 i=1,71
c        bi(i) = bi(i)*c
c        If (bi(i).le.1.e-35) GoTo 130
c120  Continue
c   change algorithm to prevent overflow/underflow
c   problems on vax....b.m.b
      r1b(69) = 10
      r2b(69) = 1+1390/z
      Do 110 i=1,68
      j = 69-i
      r1b(j) = r2b(j+1)/r1b(j+1)
      r2b(j) = 1+(2*j+1)*r1b(j)/z
 110  Continue
      r1b0 = r2b(1)/r1b(1)
      bi(1) = 0.5*(1.-exp(-2.*z))/z
      bi(2) = bi(1)/r1b0
      Do 120 i=3,71
      isv=i
      bi(i) = bi(i-1)/r1b(i-2)
      If (bi(i).le.1.e-36) GoTo 130
 120  Continue
      Return
 130  Do 140 i2=isv,71
         bi(i2) = 0.
 140  Continue
      Return
c
      End
