C
C     ..................................................................
C
C        SUBROUTINE BESK
C
C           COMPUTE THE K BESSEL FUNCTION FOR A GIVEN ARGUMENT AND ORDER
C
C        USAGE
C           CALL BESK(X,N,BK,IER)
C
C        DESCRIPTION OF PARAMETERS
C           X  -THE ARGUMENT OF THE K BESSEL FUNCTION DESIRED
C           N  -THE ORDER OF THE K BESSEL FUNCTION DESIRED
C           BK -THE RESULTANT K BESSEL FUNCTION
C           IER-RESULTANT ERROR CODE WHERE
C              IER=0  NO ERROR
C              IER=1  N IS NEGATIVE
C              IER=2  X IS ZERO OR NEGATIVE
C              IER=3  X .GT. 170, MACHINE RANGE EXCEEDED
C              IER=4  BK .GT. 10**70
C
C        REMARKS
C           N MUST BE GREATER THAN OR EQUAL TO ZERO
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           COMPUTES ZERO ORDER AND FIRST ORDER BESSEL FUNCTIONS USING
C           SERIES APPROXIMATIONS AND THEN COMPUTES N TH ORDER FUNCTION
C           USING RECURRENCE RELATION.
C           RECURRENCE RELATION AND POLYNOMIAL APPROXIMATION TECHNIQUE
C           AS DESCRIBED BY A.J.M.HITCHCOCK,'POLYNOMIAL APPROXIMATIONS
C           TO BESSEL FUNCTIONS OF ORDER ZERO AND ONE AND TO RELATED
C           FUNCTIONS', M.T.A.C., V.11,1957,PP.86-88, AND G.N. WATSON,
C           'A TREATISE ON THE THEORY OF BESSEL FUNCTIONS', CAMBRIDGE
C           UNIVERSITY PRESS, 1958, P. 62
C
C     ..................................................................
C
      SUBROUTINE BESK(X,N,BK,IER)
      DIMENSION T(12)
      BK=.0
      IF(N)10,11,11
   10 IER=1
      RETURN
   11 IF(X)12,12,20
   12 IER=2
      RETURN
   20 IF(X-170.0)22,22,21
   21 IER=3
      RETURN
   22 IER=0
      IF(X-1.)36,36,25
   25 A=EXP(-X)
      B=1./X
      C=SQRT(B)
      T(1)=B
      DO 26 L=2,12
   26 T(L)=T(L-1)*B
      IF(N-1)27,29,27
C
C     COMPUTE KO USING POLYNOMIAL APPROXIMATION
C
   27 G0=A*(1.2533141-.1566642*T(1)+.08811128*T(2)-.09139095*T(3)
     2+.1344596*T(4)-.2299850*T(5)+.3792410*T(6)-.5247277*T(7)
     3+.5575368*T(8)-.4262633*T(9)+.2184518*T(10)-.06680977*T(11)
     4+.009189383*T(12))*C
      IF(N)20,28,29
   28 BK=G0
      RETURN
C
C     COMPUTE K1 USING POLYNOMIAL APPROXIMATION
C
   29 G1=A*(1.2533141+.4699927*T(1)-.1468583*T(2)+.1280427*T(3)
     2-.1736432*T(4)+.2847618*T(5)-.4594342*T(6)+.6283381*T(7)
     3-.6632295*T(8)+.5050239*T(9)-.2581304*T(10)+.07880001*T(11)
     4-.01082418*T(12))*C
      IF(N-1)20,30,31
   30 BK=G1
      RETURN
C
C     FROM KO,K1 COMPUTE KN USING RECURRENCE RELATION
C
   31 DO 35 J=2,N
      GJ=2.*(FLOAT(J)-1.)*G1/X+G0
      IF(GJ-1.7E33)33,33,32
   32 IER=4
      GO TO 34
   33 G0=G1
   35 G1=GJ
   34 BK=GJ
      RETURN
   36 B=X/2.
      A=.5772157+ALOG(B)
      C=B*B
      IF(N-1)37,43,37
C
C     COMPUTE KO USING SERIES EXPANSION
C
   37 G0=-A
      X2J=1.
      FACT=1.
      HJ=.0
      DO 40 J=1,6
      RJ=1./FLOAT(J)
      X2J=X2J*C
      FACT=FACT*RJ*RJ
      HJ=HJ+RJ
   40 G0=G0+X2J*FACT*(HJ-A)
      IF(N)43,42,43
   42 BK=G0
      RETURN
C
C     COMPUTE K1 USING SERIES EXPANSION
C
   43 X2J=B
      FACT=1.
      HJ=1.
      G1=1./X+X2J*(.5+A-HJ)
      DO 50 J=2,8
      X2J=X2J*C
      RJ=1./FLOAT(J)
      FACT=FACT*RJ*RJ
      HJ=HJ+RJ
   50 G1=G1+X2J*FACT*(.5+(A-HJ)*FLOAT(J))
      IF(N-1)31,52,31
   52 BK=G1
      RETURN
      END