      SUBROUTINE BESI(X, ALPHA, KODE, N, Y, NZ)
C
C     WRITTEN BY D. E. AMOS AND S. L. DANIEL, JANUARY,1975.
C
C     REFERENCE
C         SAND-75-0152
C
C         CDC 6600 SUBROUTINES IBESS AND JBESS FOR BESSEL FUNCTIONS
C         I(NU,X) AND J(NU,X), X.GE.0, NU.GE.0  BY D.E. AMOS, S.L.
C         DANIEL, M.K. WESTON. ACM TRANS MATH SOFTWARE,3,PP 76-92
C         (1977)
C
C         TABLES OF BESSEL FUNCTIONS OF MODERATE OR LARGE ORDERS,
C         NPL MATHEMATICAL TABLES, VOL. 6, BY F.W.J. OLVER, HER
C         MAJESTY-S STATIONERY OFFICE, LONDON, 1962.
C
C     ABSTRACT
C         BESI COMPUTES AN N MEMBER SEQUENCE OF I BESSEL FUNCTIONS
C         I/SUB(ALPHA+K-1)/(X), K=1,...,N OR SCALED BESSEL FUNCTIONS
C         EXP(-X)*I/SUB(ALPHA+K-1)/(X), K=1,...,N FOR NON-NEGATIVE ALPHA
C         AND X. A COMBINATION OF THE POWER SERIES, THE ASYMPTOTIC
C         EXPANSION FOR X TO INFINITY, AND THE UNIFORM ASYMPTOTIC
C         EXPANSION FOR NU TO INFINITY ARE APPLIED OVER SUBDIVISIONS OF
C         THE (NU,X) PLANE. FOR VALUES NOT COVERED BY ONE OF THESE
C         FORMULAE, THE ORDER IS INCREMENTED BY AN INTEGER SO THAT ONE
C         OF THESE FORMULAE APPLY. BACKWARD RECURSION IS USED TO REDUCE
C         ORDERS BY INTEGER VALUES. THE ASYMPTOTIC EXPANSION FOR X TO
C         INFINITY IS USED ONLY WHEN THE ENTIRE SEQUENCE (SPECIFICALLY
C         THE LAST MEMBER) LIES WITHIN THE REGION COVERED BY THE
C         EXPANSION. LEADING TERMS OF THESE EXPANSIONS ARE USED TO TEST
C         FOR OVER OR UNDERFLOW WHERE APPROPRIATE. IF A SEQUENCE IS
C         REQUESTED AND THE LAST MEMBER WOULD UNDERFLOW, THE RESULT IS
C         SET TO ZERO AND THE NEXT LOWER ORDER TRIED, ETC., UNTIL A
C         MEMBER COMES ON SCALE OR ALL ARE SET TO ZERO. AN OVERFLOW
C         CANNOT OCCUR WITH SCALING.
C
C         BESI CALLS ASIK, GAMLN, SPMPAR, AND I1MACH
C
C     DESCRIPTION OF ARGUMENTS
C
C         INPUT
C           X      - X.GE.0.0E0
C           ALPHA  - ORDER OF FIRST MEMBER OF THE SEQUENCE,
C                    ALPHA.GE.0.0E0
C           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
C                    KODE=1 RETURNS
C                           Y(K)=        I/SUB(ALPHA+K-1)/(X),
C                                K=1,...,N
C                    KODE=2 RETURNS
C                           Y(K)=EXP(-X)*I/SUB(ALPHA+K-1)/(X),
C                                K=1,...,N
C           N      - NUMBER OF MEMBERS IN THE SEQUENCE, N.GE.1
C
C         OUTPUT
C           Y      - A VECTOR WHOSE FIRST N COMPONENTS CONTAIN
C                    VALUES FOR I/SUB(ALPHA+K-1)/(X) OR SCALED
C                    VALUES FOR EXP(-X)*I/SUB(ALPHA+K-1)/(X),
C                    K=1,...,N DEPENDING ON KODE
C           NZ     - ERROR INDICATOR
C                    NZ= 0     NORMAL RETURN-COMPUTATION COMPLETED
C                    NZ=-1     X IS LESS THAN 0.0
C                    NZ=-2     ALPHA IS LESS THAN 0.0
C                    NZ=-3     N IS LESS THAN 1
C                    NZ=-4     KODE IS NOT 1 OR 2
C                    NZ=-5     X IS TOO LARGE FOR KODE=1
C                    NZ.GT.0   LAST NZ COMPONENTS OF Y SET TO 0.0
C                              BECAUSE OF UNDERFLOW
C
C     ERROR CONDITIONS
C         IMPROPER INPUT ARGUMENTS - A FATAL ERROR
C         OVERFLOW WITH KODE=1 - A FATAL ERROR
C         UNDERFLOW - A NON-FATAL ERROR(NZ.GT.0)
C
      INTEGER I, IALP, IN, INLIM, IS, I1, I2, K, KK, KM, KODE, KT,
     * N, NN, NS, NZ
      INTEGER I1MACH
      REAL AIN, AK, AKM, ALPHA, ANS, AP, ARG, ATOL, TOLLN, DFN,
     * DTM, DX, EARG, ELIM, ETX, FLGIK,FN, FNF, FNI,FNP1,FNU,GLN,RA,
     * RTTPI, S, SX, SXO2, S1, S2, T, TA, TB, TEMP, TFN, TM, TOL,
     * TRX, T2, X, XO2, XO2L, Y, Z
      REAL GAMLN, SPMPAR
      DIMENSION Y(1), TEMP(3)
      DATA RTTPI           / 3.98942280401433E-01/
      DATA INLIM           /          80         /
C     -------------------
C     I1MACH(15) REPLACES I1MACH(12) IN A DOUBLE PRECISION CODE
C     I1MACH(14) REPLACES I1MACH(11) IN A DOUBLE PRECISION CODE
C
C     DEFINITION OF THE TOLERANCES TOL AND ELIM
C
      TB = I1MACH(10)
      TA = SPMPAR(1)/TB
      IF (TB.EQ.2.0E0) GO TO 1
      IF (TB.EQ.8.0E0) GO TO 2
      IF (TB.EQ.16.0E0) GO TO 3
      TB = ALOG(TB)
      GO TO 5
    1 TB = .69315E0
      GO TO 5
    2 TB = 2.07944E0
      GO TO 5
    3 TB = 2.77259E0
C
    5 TOL = AMAX1(TA,1.E-15)
      I1 = I1MACH(11)
      I2 = I1MACH(12)
C     LN(10**3) = 6.90776
      ELIM = FLOAT(-I2)*TB - 6.90776E0
C     TOLLN = -LN(TOL)
      TOLLN = FLOAT(I1)*TB
      TOLLN = AMIN1(TOLLN,34.5388E0)
C
C
C
      NZ = 0
      KT = 1
      IF (N-1) 590, 10, 20
   10 KT = 2
   20 NN = N
      IF (KODE.LT.1 .OR. KODE.GT.2) GO TO 570
      IF (X) 600, 30, 80
   30 IF (ALPHA) 580, 40, 50
   40 Y(1) = 1.0E0
      IF (N.EQ.1) RETURN
      I1 = 2
      GO TO 60
   50 I1 = 1
   60 DO 70 I=I1,N
        Y(I) = 0.0E0
   70 CONTINUE
      RETURN
   80 CONTINUE
      IF (ALPHA.LT.0.0E0) GO TO 580
C
      IALP = INT(ALPHA)
      FNI = FLOAT(IALP+N-1)
      FNF = ALPHA - FLOAT(IALP)
      DFN = FNI + FNF
      FNU = DFN
      IN = 0
      XO2 = X*0.5E0
      SXO2 = XO2*XO2
      ETX = FLOAT(KODE-1)
      SX = ETX*X
C
C     DECISION TREE FOR REGION WHERE SERIES, ASYMPTOTIC EXPANSION FOR X
C     TO INFINITY AND ASYMPTOTIC EXPANSION FOR NU TO INFINITY ARE
C     APPLIED.
C
      IF (SXO2.LE.(FNU+1.0E0)) GO TO 90
      IF (X.LE.12.0E0) GO TO 110
      FN = 0.55E0*FNU*FNU
      FN = AMAX1(17.0E0,FN)
      IF (X.GE.FN) GO TO 430
      ANS = AMAX1(36.0E0-FNU,0.0E0)
      NS = INT(ANS)
      FNI = FNI + FLOAT(NS)
      DFN = FNI + FNF
      FN = DFN
      IS = KT
      KM = N - 1 + NS
      IF (KM.GT.0) IS = 3
      GO TO 120
   90 FN = FNU
      FNP1 = FN + 1.0E0
      XO2L = ALOG(XO2)
      IS = KT
      IF (X.LE.0.5E0) GO TO 230
      NS = 0
  100 FNI = FNI + FLOAT(NS)
      DFN = FNI + FNF
      FN = DFN
      FNP1 = FN + 1.0E0
      IS = KT
      IF (N-1+NS.GT.0) IS = 3
      GO TO 230
  110 XO2L = ALOG(XO2)
      NS = INT(SXO2-FNU)
      GO TO 100
  120 CONTINUE
C
C     OVERFLOW TEST ON UNIFORM ASYMPTOTIC EXPANSION
C
      IF (KODE.EQ.2) GO TO 130
      IF (ALPHA.LT.1.0E0) GO TO 150
      Z = X/ALPHA
      RA = SQRT(1.0E0+Z*Z)
      GLN = ALOG((1.0E0+RA)/Z)
      T = RA*(1.0E0-ETX) + ETX/(Z+RA)
      ARG = ALPHA*(T-GLN)
      IF (ARG.GT.ELIM) GO TO 610
      IF (KM.EQ.0) GO TO 140
  130 CONTINUE
C
C     UNDERFLOW TEST ON UNIFORM ASYMPTOTIC EXPANSION
C
      Z = X/FN
      RA = SQRT(1.0E0+Z*Z)
      GLN = ALOG((1.0E0+RA)/Z)
      T = RA*(1.0E0-ETX) + ETX/(Z+RA)
      ARG = FN*(T-GLN)
  140 IF (ARG.LT.(-ELIM)) GO TO 280
      GO TO 190
  150 IF (X.GT.ELIM) GO TO 610
      GO TO 130
C
C     UNIFORM ASYMPTOTIC EXPANSION FOR NU TO INFINITY
C
  160 IF (KM.NE.0) GO TO 170
      Y(1) = TEMP(3)
      RETURN
  170 TEMP(1) = TEMP(3)
      IN = NS
      KT = 1
      I1 = 0
  180 CONTINUE
      IS = 2
      FNI = FNI - 1.0E0
      DFN = FNI + FNF
      FN = DFN
      IF(I1.EQ.2) GO TO 350
      Z = X/FN
      RA = SQRT(1.0E0+Z*Z)
      GLN = ALOG((1.0E0+RA)/Z)
      T = RA*(1.0E0-ETX) + ETX/(Z+RA)
      ARG = FN*(T-GLN)
  190 CONTINUE
      I1 = IABS(3-IS)
      I1 = MAX0(I1,1)
      FLGIK = 1.0E0
      CALL ASIK(X,FN,KODE,FLGIK,RA,ARG,I1,TOL,TEMP(IS))
      GO TO (180, 350, 510), IS
C
C     SERIES FOR (X/2)**2.LE.NU+1
C
  230 CONTINUE
      GLN = GAMLN(FNP1)
      ARG = FN*XO2L - GLN - SX
      IF (ARG.LT.(-ELIM)) GO TO 300
      EARG = EXP(ARG)
  240 CONTINUE
      S = 1.0E0
      IF (X.LT.TOL) GO TO 260
      AK = 3.0E0
      T2 = 1.0E0
      T = 1.0E0
      S1 = FN
      DO 250 K=1,17
        S2 = T2 + S1
        T = T*SXO2/S2
        S = S + T
        IF (ABS(T).LT.TOL) GO TO 260
        T2 = T2 + AK
        AK = AK + 2.0E0
        S1 = S1 + FN
  250 CONTINUE
  260 CONTINUE
      TEMP(IS) = S*EARG
      GO TO (270, 350, 500), IS
  270 EARG = EARG*FN/XO2
      FNI = FNI - 1.0E0
      DFN = FNI + FNF
      FN = DFN
      IS = 2
      GO TO 240
C
C     SET UNDERFLOW VALUE AND UPDATE PARAMETERS
C
  280 Y(NN) = 0.0E0
      NN = NN - 1
      FNI = FNI - 1.0E0
      DFN = FNI + FNF
      FN = DFN
      IF (NN-1) 340, 290, 130
  290 KT = 2
      IS = 2
      GO TO 130
  300 Y(NN) = 0.0E0
      NN = NN - 1
      FNP1 = FN
      FNI = FNI - 1.0E0
      DFN = FNI + FNF
      FN = DFN
      IF (NN-1) 340, 310, 320
  310 KT = 2
      IS = 2
  320 IF (SXO2.LE.FNP1) GO TO 330
      GO TO 130
  330 ARG = ARG - XO2L + ALOG(FNP1)
      IF (ARG.LT.(-ELIM)) GO TO 300
      GO TO 230
  340 NZ = N - NN
      RETURN
C
C     BACKWARD RECURSION SECTION
C
  350 CONTINUE
      NZ = N - NN
  360 CONTINUE
      IF(KT.EQ.2) GO TO 420
      S1 = TEMP(1)
      S2 = TEMP(2)
      TRX = 2.0E0/X
      DTM = FNI
      TM = (DTM+FNF)*TRX
      IF (IN.EQ.0) GO TO 390
C     BACKWARD RECUR TO INDEX ALPHA+NN-1
      DO 380 I=1,IN
        S = S2
        S2 = TM*S2 + S1
        S1 = S
        DTM = DTM - 1.0E0
        TM = (DTM+FNF)*TRX
  380 CONTINUE
      Y(NN) = S1
      IF (NN.EQ.1) RETURN
      Y(NN-1) = S2
      IF (NN.EQ.2) RETURN
      GO TO 400
  390 CONTINUE
C     BACKWARD RECUR FROM INDEX ALPHA+NN-1 TO ALPHA
      Y(NN) = S1
      Y(NN-1) = S2
      IF (NN.EQ.2) RETURN
  400 K = NN + 1
      DO 410 I=3,NN
        K = K - 1
        Y(K-2) = TM*Y(K-1) + Y(K)
        DTM = DTM - 1.0E0
        TM = (DTM+FNF)*TRX
  410 CONTINUE
      RETURN
  420 Y(1) = TEMP(2)
      RETURN
C
C     ASYMPTOTIC EXPANSION FOR X TO INFINITY
C
  430 CONTINUE
      EARG = RTTPI/SQRT(X)
      IF (KODE.EQ.2) GO TO 440
      IF (X.GT.ELIM) GO TO 610
      EARG = EARG*EXP(X)
  440 ETX = 8.0E0*X
      IS = KT
      IN = 0
      FN = FNU
  450 DX = FNI + FNI
      TM = 0.0E0
      IF (FNI.EQ.0.0E0 .AND. ABS(FNF).LT.TOL) GO TO 460
      TM = 4.0E0*FNF*(FNI+FNI+FNF)
  460 CONTINUE
      DTM = DX*DX
      S1 = ETX
      TRX = DTM - 1.0E0
      DX = -(TRX+TM)/ETX
      T = DX
      S = 1.0E0 + DX
      ATOL = TOL*ABS(S)
      S2 = 1.0E0
      AK = 8.0E0
      DO 470 K=1,25
        S1 = S1 + ETX
        S2 = S2 + AK
        DX = DTM - S2
        AP = DX + TM
        T = -T*AP/S1
        S = S + T
        IF (ABS(T).LE.ATOL) GO TO 480
        AK = AK + 8.0E0
  470 CONTINUE
  480 TEMP(IS) = S*EARG
      IF(IS.EQ.2) GO TO 360
      IS = 2
      FNI = FNI - 1.0E0
      DFN = FNI + FNF
      FN = DFN
      GO TO 450
C
C     BACKWARD RECURSION WITH NORMALIZATION BY
C     ASYMPTOTIC EXPANSION FOR NU TO INFINITY OR POWER SERIES.
C
  500 CONTINUE
C     COMPUTATION OF LAST ORDER FOR SERIES NORMALIZATION
      AKM = AMAX1(3.0E0-FN,0.0E0)
      KM = INT(AKM)
      TFN = FN + FLOAT(KM)
      TA = (GLN+TFN-0.9189385332E0-0.0833333333E0/TFN)/(TFN+0.5E0)
      TA = XO2L - TA
      TB = -(1.0E0-1.0E0/TFN)/TFN
      AIN = TOLLN/(-TA+SQRT(TA*TA-TOLLN*TB)) + 1.5E0
      IN = INT(AIN)
      IN = IN + KM
      GO TO 520
  510 CONTINUE
C     COMPUTATION OF LAST ORDER FOR ASYMPTOTIC EXPANSION NORMALIZATION
      T = 1.0E0/(FN*RA)
      AIN = TOLLN/(GLN+SQRT(GLN*GLN+T*TOLLN)) + 1.5E0
      IN = INT(AIN)
      IF (IN.GT.INLIM) GO TO 160
  520 CONTINUE
      TRX = 2.0E0/X
      DTM = FNI + FLOAT(IN)
      TM = (DTM+FNF)*TRX
      TA = 0.0E0
      TB = TOL
      KK = 1
  530 CONTINUE
C
C     BACKWARD RECUR UNINDEXED
C
      DO 540 I=1,IN
        S = TB
        TB = TM*TB + TA
        TA = S
        DTM = DTM - 1.0E0
        TM = (DTM+FNF)*TRX
  540 CONTINUE
C     NORMALIZATION
      IF (KK.NE.1) GO TO 550
      TA = (TA/TB)*TEMP(3)
      TB = TEMP(3)
      KK = 2
      IN = NS
      IF (NS.NE.0) GO TO 530
  550 Y(NN) = TB
      NZ = N - NN
      IF (NN.EQ.1) RETURN
      TB = TM*TB + TA
      K = NN - 1
      Y(K) = TB
      IF (NN.EQ.2) RETURN
      DTM = DTM - 1.0E0
      TM = (DTM+FNF)*TRX
      KM = K - 1
C
C     BACKWARD RECUR INDEXED
C
      DO 560 I=1,KM
        Y(K-1) = TM*Y(K) + Y(K+1)
        DTM = DTM - 1.0E0
        TM = (DTM+FNF)*TRX
        K = K - 1
  560 CONTINUE
      RETURN
C
C
C
  570 CONTINUE
      NZ = -4
      RETURN
  580 CONTINUE
      NZ = -2
      RETURN
  590 CONTINUE
      NZ = -3
      RETURN
  600 CONTINUE
      NZ = -1
      RETURN
  610 CONTINUE
      NZ = -5
      RETURN
      END
