      SUBROUTINE BESLIA (ZZ,BI)
C *** SEE R H 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.) GO TO 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) GO TO 20
 10      CONTINUE
 20      BI(I) = A*C
         IF (A.LE.1.E-35) GO TO 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.) GO TO 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) GO TO 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) GO TO 130
 120  CONTINUE
      RETURN
 130  DO 140 I2=ISV,71
         BI(I2) = 0.
 140  CONTINUE
      RETURN
C
      END
