*DECK BSKIN
      SUBROUTINE BSKIN (X, N, KODE, M, Y, NZ, IERR)
C***BEGIN PROLOGUE  BSKIN
C***PURPOSE  Compute repeated integrals of the K-zero Bessel function.
C***LIBRARY   SLATEC
C***CATEGORY  C10F
C***TYPE      SINGLE PRECISION (BSKIN-S, DBSKIN-D)
C***KEYWORDS  BICKLEY FUNCTIONS, EXPONENTIAL INTEGRAL,
C             INTEGRALS OF BESSEL FUNCTIONS, K-ZERO BESSEL FUNCTION
C***AUTHOR  Amos, D. E., (SNLA)
C***DESCRIPTION
C
C         The following definitions are used in BSKIN:
C
C   Definition 1
C         KI(0,X) = K-zero Bessel function.
C
C   Definition 2
C         KI(N,X) = Bickley Function
C                 =  integral from X to infinity of KI(N-1,t)dt
C                     for X .ge. 0 and N = 1,2,...
C   ____________________________________________________________________
C      BSKIN computes sequences of Bickley functions (repeated integrals
C      of the K0 Bessel function); i.e. for fixed X and N and K=1,...,
C      BSKIN computes the M-member sequence
C
C                     Y(K) =        KI(N+K-1,X) for KODE=1
C      or
C                     Y(K) = EXP(X)*KI(N+K-1,X) for KODE=2,
C
C      for N.ge.0 and X.ge.0 (N and X cannot be zero simultaneously).
C
C      INPUT
C        X      - Argument, X .ge. 0.0E0
C        N      - Order of first member of the sequence N .ge. 0
C        KODE   - Selection parameter
C                 KODE = 1 returns Y(K)=       KI(N+K-1,X), K=1,M
C                      = 2 returns Y(K)=EXP(X)*KI(N+K-1,X), K=1,M
C        M      - Number of members in the sequence, M.ge.1
C
C      OUTPUT
C        Y      - A vector of dimension at least M containing the
C                 sequence selected by KODE.
C        NZ     - Underflow flag
C                 NZ = 0 means computation completed
C                    = M means an exponential underflow occurred on
C                        KODE=1.  Y(K)=0.0E0, K=1,...,M is returned
C        IERR   - Error flag
C                 IERR = 0, Normal return, computation completed.
C                      = 1, Input error,   no computation.
C                      = 2, Error,         no computation.  The
C                           termination condition was not met.
C
C      The nominal computational accuracy is the maximum of unit
C      roundoff (=R1MACH(4)) and 1.0e-18 since critical constants
C      are given to only 18 digits.
C
C      DBSKIN is the double precision version of BSKIN.
C
C *Long Description:
C
C         Numerical recurrence on
C
C      (L-1)*KI(L,X) = X(KI(L-3,X) - KI(L-1,X)) + (L-2)*KI(L-2,X)
C
C         is stable where recurrence is carried forward or backward
C         away from INT(X+0.5).  The power series for indices 0,1 and 2
C         on 0.le.X.le. 2 starts a stable recurrence for indices
C         greater than 2.  If N is sufficiently large (N.gt.NLIM), the
C         uniform asymptotic expansion for N to INFINITY is more
C         economical.  On X.gt.2 the recursion is started by evaluating
C         the uniform expansion for the three members whose indices are
C         closest to INT(X+0.5) within the set N,...,N+M-1.  Forward
C         recurrence, backward recurrence or both, complete the
C         sequence depending on the relation of INT(X+0.5) to the
C         indices N,...,N+M-1.
C
C***REFERENCES  D. E. Amos, Uniform asymptotic expansions for
C                 exponential integrals E(N,X) and Bickley functions
C                 KI(N,X), ACM Transactions on Mathematical Software,
C                 1983.
C               D. E. Amos, A portable Fortran subroutine for the
C                 Bickley functions KI(N,X), Algorithm 609, ACM
C                 Transactions on Mathematical Software, 1983.
C***ROUTINES CALLED  BKIAS, BKISR, EXINT, GAMRN, I1MACH, R1MACH
C***REVISION HISTORY  (YYMMDD)
C   820601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   891009  Removed unreferenced statement label.  (WRB)
C   891009  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  BSKIN
      INTEGER I, ICASE, IERR, IL, I1M, K, KK, KODE, KTRMS, M,
     * M3, N, NE, NFLG, NL, NLIM, NN, NP, NS, NT, NZ
      INTEGER I1MACH
      REAL A, ENLIM, EXI, FN, GR, H, HN, HRTPI, SS, TOL, T1, T2, W, X,
     * XLIM, XNLIM, XP, Y, YS, YSS
      REAL GAMRN, R1MACH
      DIMENSION EXI(102), A(50), YS(3), YSS(3), H(31), Y(*)
      SAVE A, HRTPI
C-----------------------------------------------------------------------
C             COEFFICIENTS IN SERIES OF EXPONENTIAL INTEGRALS
C-----------------------------------------------------------------------
      DATA A(1), A(2), A(3), A(4), A(5), A(6), A(7), A(8), A(9), A(10),
     * A(11), A(12), A(13), A(14), A(15), A(16), A(17), A(18), A(19),
     * A(20), A(21), A(22), A(23), A(24) /1.00000000000000000E+00,
     * 5.00000000000000000E-01,3.75000000000000000E-01,
     * 3.12500000000000000E-01,2.73437500000000000E-01,
     * 2.46093750000000000E-01,2.25585937500000000E-01,
     * 2.09472656250000000E-01,1.96380615234375000E-01,
     * 1.85470581054687500E-01,1.76197052001953125E-01,
     * 1.68188095092773438E-01,1.61180257797241211E-01,
     * 1.54981017112731934E-01,1.49445980787277222E-01,
     * 1.44464448094367981E-01,1.39949934091418982E-01,
     * 1.35833759559318423E-01,1.32060599571559578E-01,
     * 1.28585320635465905E-01,1.25370687619579257E-01,
     * 1.22385671247684513E-01,1.19604178719328047E-01,
     * 1.17004087877603524E-01/
      DATA A(25), A(26), A(27), A(28), A(29), A(30), A(31), A(32),
     * A(33), A(34), A(35), A(36), A(37), A(38), A(39), A(40), A(41),
     * A(42), A(43), A(44), A(45), A(46), A(47), A(48)
     * /1.14566502713486784E-01,1.12275172659217048E-01,
     * 1.10116034723462874E-01,1.08076848895250599E-01,
     * 1.06146905164978267E-01,1.04316786110409676E-01,
     * 1.02578173008569515E-01,1.00923686347140974E-01,
     * 9.93467537479668965E-02,9.78414999033007314E-02,
     * 9.64026543164874854E-02,9.50254735405376642E-02,
     * 9.37056752969190855E-02,9.24393823875012600E-02,
     * 9.12230747245078224E-02,9.00535481254756708E-02,
     * 8.89278787739072249E-02,8.78433924473961612E-02,
     * 8.67976377754033498E-02,8.57883629175498224E-02,
     * 8.48134951571231199E-02,8.38711229887106408E-02,
     * 8.29594803475290034E-02,8.20769326842574183E-02/
      DATA A(49), A(50) /8.12219646354630702E-02,8.03931690779583449E-02
     * /
C-----------------------------------------------------------------------
C             SQRT(PI)/2
C-----------------------------------------------------------------------
      DATA HRTPI /8.86226925452758014E-01/
C
C***FIRST EXECUTABLE STATEMENT  BSKIN
      IERR = 0
      NZ=0
      IF (X.LT.0.0E0) IERR=1
      IF (N.LT.0) IERR=1
      IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1
      IF (M.LT.1) IERR=1
      IF (X.EQ.0.0E0 .AND. N.EQ.0) IERR=1
      IF (IERR.NE.0) RETURN
      IF (X.EQ.0.0E0) GO TO 300
      I1M = -I1MACH(12)
      T1 = 2.3026E0*R1MACH(5)*I1M
      XLIM = T1 - 3.228086E0
      T2 = T1 + N + M - 1
      IF (T2.GT.1000.0E0) XLIM = T1 - 0.5E0*(LOG(T2)-0.451583E0)
      IF (X.GT.XLIM .AND. KODE.EQ.1) GO TO 320
      TOL = MAX(R1MACH(4),1.0E-18)
      I1M = I1MACH(11)
C-----------------------------------------------------------------------
C     LN(NLIM) = 0.125*LN(EPS),   NLIM = 2*KTRMS+N
C-----------------------------------------------------------------------
      XNLIM = 0.287823E0*(I1M-1)*R1MACH(5)
      ENLIM = EXP(XNLIM)
      NLIM = INT(ENLIM) + 2
      NLIM = MIN(100,NLIM)
      NLIM = MAX(20,NLIM)
      M3 = MIN(M,3)
      NL = N + M - 1
      IF (X.GT.2.0E0) GO TO 130
      IF (N.GT.NLIM) GO TO 280
C-----------------------------------------------------------------------
C     COMPUTATION BY SERIES FOR 0.LE.X.LE.2
C-----------------------------------------------------------------------
      NFLG = 0
      NN = N
      IF (NL.LE.2) GO TO 60
      M3 = 3
      NN = 0
      NFLG = 1
   60 CONTINUE
      XP = 1.0E0
      IF (KODE.EQ.2) XP = EXP(X)
      DO 80 I=1,M3
        CALL BKISR(X, NN, W, IERR)
      IF(IERR.NE.0) RETURN
        W = W*XP
        IF (NN.LT.N) GO TO 70
        KK = NN - N + 1
        Y(KK) = W
   70   CONTINUE
        YS(I) = W
        NN = NN + 1
   80 CONTINUE
      IF (NFLG.EQ.0) RETURN
      NS = NN
      XP = 1.0E0
   90 CONTINUE
C-----------------------------------------------------------------------
C     FORWARD RECURSION SCALED BY EXP(X) ON ICASE=0,1,2
C-----------------------------------------------------------------------
      FN = NS - 1
      IL = NL - NS + 1
      IF (IL.LE.0) RETURN
      DO 110 I=1,IL
        T1 = YS(2)
        T2 = YS(3)
        YS(3) = (X*(YS(1)-YS(3))+(FN-1.0E0)*YS(2))/FN
        YS(2) = T2
        YS(1) = T1
        FN = FN + 1.0E0
        IF (NS.LT.N) GO TO 100
        KK = NS - N + 1
        Y(KK) = YS(3)*XP
  100   CONTINUE
        NS = NS + 1
  110 CONTINUE
      RETURN
C-----------------------------------------------------------------------
C     COMPUTATION BY ASYMPTOTIC EXPANSION FOR X.GT.2
C-----------------------------------------------------------------------
  130 CONTINUE
      W = X + 0.5E0
      NT = INT(W)
      IF (NL.GT.NT) GO TO 270
C-----------------------------------------------------------------------
C     CASE NL.LE.NT, ICASE=0
C-----------------------------------------------------------------------
      ICASE = 0
      NN = NL
      NFLG = MIN(M-M3,1)
  140 CONTINUE
      KK = (NLIM-NN)/2
      KTRMS = MAX(0,KK)
      NS = NN + 1
      NP = NN - M3 + 1
      XP = 1.0E0
      IF (KODE.EQ.1) XP = EXP(-X)
      DO 150 I=1,M3
        KK = I
        CALL BKIAS(X, NP, KTRMS, A, W, KK, NE, GR, H, IERR)
      IF(IERR.NE.0) RETURN
        YS(I) = W
        NP = NP + 1
  150 CONTINUE
C-----------------------------------------------------------------------
C     SUM SERIES OF EXPONENTIAL INTEGRALS BACKWARD
C-----------------------------------------------------------------------
      IF (KTRMS.EQ.0) GO TO 160
      NE = KTRMS + KTRMS + 1
      NP = NN - M3 + 2
      CALL EXINT(X, NP, 2, NE, TOL, EXI, NZ, IERR)
      IF(NZ.NE.0) GO TO 320
      IF(IERR.EQ.2) RETURN
  160 CONTINUE
      DO 190 I=1,M3
        SS = 0.0E0
        IF (KTRMS.EQ.0) GO TO 180
        KK = I + KTRMS + KTRMS - 2
        IL = KTRMS
        DO 170 K=1,KTRMS
          SS = SS + A(IL)*EXI(KK)
          KK = KK - 2
          IL = IL - 1
  170   CONTINUE
  180   CONTINUE
        YS(I) = YS(I) + SS
  190 CONTINUE
      IF (ICASE.EQ.1) GO TO 200
      IF (NFLG.NE.0) GO TO 220
  200 CONTINUE
      DO 210 I=1,M3
        Y(I) = YS(I)*XP
  210 CONTINUE
      IF (ICASE.EQ.1 .AND. NFLG.EQ.1) GO TO 90
      RETURN
  220 CONTINUE
C-----------------------------------------------------------------------
C     BACKWARD RECURSION SCALED BY EXP(X) ICASE=0,2
C-----------------------------------------------------------------------
      KK = NN - N + 1
      K = M3
      DO 230 I=1,M3
        Y(KK) = YS(K)*XP
        YSS(I) = YS(I)
        KK = KK - 1
        K = K - 1
  230 CONTINUE
      IL = KK
      IF (IL.LE.0) GO TO 250
      FN = NN - 3
      DO 240 I=1,IL
        T1 = YS(2)
        T2 = YS(1)
        YS(1) = YS(2) + ((FN+2.0E0)*YS(3)-(FN+1.0E0)*YS(1))/X
        YS(2) = T2
        YS(3) = T1
        Y(KK) = YS(1)*XP
        KK = KK - 1
        FN = FN - 1.0E0
  240 CONTINUE
  250 CONTINUE
      IF (ICASE.NE.2) RETURN
      DO 260 I=1,M3
        YS(I) = YSS(I)
  260 CONTINUE
      GO TO 90
  270 CONTINUE
      IF (N.LT.NT) GO TO 290
C-----------------------------------------------------------------------
C     ICASE=1, NT.LE.N.LE.NL WITH FORWARD RECURSION
C-----------------------------------------------------------------------
  280 CONTINUE
      NN = N + M3 - 1
      NFLG = MIN(M-M3,1)
      ICASE = 1
      GO TO 140
C-----------------------------------------------------------------------
C     ICASE=2, N.LT.NT.LT.NL WITH BOTH FORWARD AND BACKWARD RECURSION
C-----------------------------------------------------------------------
  290 CONTINUE
      NN = NT + 1
      NFLG = MIN(M-M3,1)
      ICASE = 2
      GO TO 140
C-----------------------------------------------------------------------
C     X=0 CASE
C-----------------------------------------------------------------------
  300 CONTINUE
      FN = N
      HN = 0.5E0*FN
      GR = GAMRN(HN)
      Y(1) = HRTPI*GR
      IF (M.EQ.1) RETURN
      Y(2) = HRTPI/(HN*GR)
      IF (M.EQ.2) RETURN
      DO 310 K=3,M
        Y(K) = FN*Y(K-2)/(FN+1.0E0)
        FN = FN + 1.0E0
  310 CONTINUE
      RETURN
C-----------------------------------------------------------------------
C     UNDERFLOW ON KODE=1, X.GT.XLIM
C-----------------------------------------------------------------------
  320 CONTINUE
      NZ=M
      DO 330 I=1,M
        Y(I) = 0.0E0
  330 CONTINUE
      RETURN
      END