      SUBROUTINE ASIK(X,FNU,KODE,FLGIK,RA,ARG,IN,TOL,Y)
C
C                  ASIK COMPUTES BESSEL FUNCTIONS I AND K
C                  FOR ARGUMENTS X.GT.0.0 AND ORDERS FNU.GE.35
C                  ON FLGIK = 1 AND FLGIK = -1 RESPECTIVELY.
C
C                                    INPUT
C
C      X    - ARGUMENT, X.GT.0.0E0
C      FNU  - ORDER OF FIRST BESSEL FUNCTION
C      KODE - A PARAMETER TO INDICATE THE SCALING OPTION
C             KODE=1 RETURNS Y(I)=        I/SUB(FNU+I-1)/(X), I=1,IN
C                    OR      Y(I)=        K/SUB(FNU+I-1)/(X), I=1,IN
C                    ON FLGIK = 1.0E0 OR FLGIK = -1.0E0
C             KODE=2 RETURNS Y(I)=EXP(-X)*I/SUB(FNU+I-1)/(X), I=1,IN
C                    OR      Y(I)=EXP( X)*K/SUB(FNU+I-1)/(X), I=1,IN
C                    ON FLGIK = 1.0E0 OR FLGIK = -1.0E0
C     FLGIK - SELECTION PARAMETER FOR I OR K FUNCTION
C             FLGIK =  1.0E0 GIVES THE I FUNCTION
C             FLGIK = -1.0E0 GIVES THE K FUNCTION
C        RA - SQRT(1.+Z*Z), Z=X/FNU
C       ARG - ARGUMENT OF THE LEADING EXPONENTIAL
C        IN - NUMBER OF FUNCTIONS DESIRED, IN=1 OR 2
C       TOL - TOLERANCE SPECIFIED BY BESI OR BESK
C
C                                    OUTPUT
C
C         Y - A VECTOR WHOSE FIRST IN COMPONENTS CONTAIN THE SEQUENCE
C
C                                 WRITTEN BY
C                                 D. E. AMOS
C
C     ABSTRACT
C         ASIK IMPLEMENTS THE UNIFORM ASYMPTOTIC EXPANSION OF
C         THE I AND K BESSEL FUNCTIONS FOR FNU.GE.35 AND REAL
C         X.GT.0.0E0. THE FORMS ARE IDENTICAL EXCEPT FOR A CHANGE
C         IN SIGN OF SOME OF THE TERMS. THIS CHANGE IN SIGN IS
C         ACCOMPLISHED BY MEANS OF THE FLAG FLGIK = 1 OR -1.
C
      INTEGER IN, J, JN, K, KK, KODE, L
      REAL AK,AP,ARG,C, COEF,CON,ETX,FLGIK,FN, FNU,GLN,RA,S1,S2,
     * T, TOL, T2, X, Y, Z
      DIMENSION Y(1), C(65), CON(2)
      DATA CON(1), CON(2)  /
     1        3.98942280401432678E-01,    1.25331413731550025E+00/
      DATA C(1), C(2), C(3), C(4), C(5), C(6), C(7), C(8), C(9), C(10),
     1     C(11), C(12), C(13), C(14), C(15), C(16), C(17), C(18),
     2     C(19), C(20), C(21), C(22), C(23), C(24)/
     3       -2.08333333333333E-01,        1.25000000000000E-01,
     4        3.34201388888889E-01,       -4.01041666666667E-01,
     5        7.03125000000000E-02,       -1.02581259645062E+00,
     6        1.84646267361111E+00,       -8.91210937500000E-01,
     7        7.32421875000000E-02,        4.66958442342625E+00,
     8       -1.12070026162230E+01,        8.78912353515625E+00,
     9       -2.36408691406250E+00,        1.12152099609375E-01,
     A       -2.82120725582002E+01,        8.46362176746007E+01,
     B       -9.18182415432400E+01,        4.25349987453885E+01,
     C       -7.36879435947963E+00,        2.27108001708984E-01,
     D        2.12570130039217E+02,       -7.65252468141182E+02,
     E        1.05999045252800E+03,       -6.99579627376133E+02/
      DATA C(25), C(26), C(27), C(28), C(29), C(30), C(31), C(32),
     1     C(33), C(34), C(35), C(36), C(37), C(38), C(39), C(40),
     2     C(41), C(42), C(43), C(44), C(45), C(46), C(47), C(48)/
     3        2.18190511744212E+02,       -2.64914304869516E+01,
     4        5.72501420974731E-01,       -1.91945766231841E+03,
     5        8.06172218173731E+03,       -1.35865500064341E+04,
     6        1.16553933368645E+04,       -5.30564697861340E+03,
     7        1.20090291321635E+03,       -1.08090919788395E+02,
     8        1.72772750258446E+00,        2.02042913309661E+04,
     9       -9.69805983886375E+04,        1.92547001232532E+05,
     A       -2.03400177280416E+05,        1.22200464983017E+05,
     B       -4.11926549688976E+04,        7.10951430248936E+03,
     C       -4.93915304773088E+02,        6.07404200127348E+00,
     D       -2.42919187900551E+05,        1.31176361466298E+06,
     E       -2.99801591853811E+06,        3.76327129765640E+06/
      DATA C(49), C(50), C(51), C(52), C(53), C(54), C(55), C(56),
     1     C(57), C(58), C(59), C(60), C(61), C(62), C(63), C(64),
     2     C(65)/
     3       -2.81356322658653E+06,        1.26836527332162E+06,
     4       -3.31645172484564E+05,        4.52187689813627E+04,
     5       -2.49983048181121E+03,        2.43805296995561E+01,
     6        3.28446985307204E+06,       -1.97068191184322E+07,
     7        5.09526024926646E+07,       -7.41051482115327E+07,
     8        6.63445122747290E+07,       -3.75671766607634E+07,
     9        1.32887671664218E+07,       -2.78561812808645E+06,
     A        3.08186404612662E+05,       -1.38860897537170E+04,
     B        1.10017140269247E+02/
C     ---------------------
      FN = FNU
      Z  = (3.0E0-FLGIK)/2.0E0
      KK = INT(Z)
      DO 50 JN=1,IN
        IF (JN.EQ.1) GO TO 10
        FN = FN - FLGIK
        Z = X/FN
        RA = SQRT(1.0E0+Z*Z)
        GLN = ALOG((1.0E0+RA)/Z)
        ETX = FLOAT(KODE-1)
        T = RA*(1.0E0-ETX) + ETX/(Z+RA)
        ARG = FN*(T-GLN)*FLGIK
   10   COEF = EXP(ARG)
        T = 1.0E0/RA
        T2 = T*T
        T = T/FN
        T = SIGN(T,FLGIK)
        S2 = 1.0E0
        AP = 1.0E0
        L = 0
        DO 30 K=2,11
          L = L + 1
          S1 = C(L)
          DO 20 J=2,K
            L = L + 1
            S1 = S1*T2 + C(L)
   20     CONTINUE
          AP = AP*T
          AK = AP*S1
          S2 = S2 + AK
          IF (AMAX1(ABS(AK),ABS(AP)) .LT. TOL) GO TO 40
   30   CONTINUE
   40   CONTINUE
      T = ABS(T)
      Y(JN) = S2*COEF*SQRT(T)*CON(KK)
   50 CONTINUE
      RETURN
      END
