C
      REAL FUNCTION BESL21*8 (X)
C
C   FUNCTION            - BY CHEBYSHEV INTERPOLATION, WE OBTAIN THE
C                           VALUE OF THE BESSEL FUNCTION OF THE SECOND
C                           KIND OF ORDER ONE.
C
C   NOTICE              - WHEN !X! DOES NOT SATISFY 8.797D-77<=!X!<
C                           1.0D12, ERROR MESSAGE IS PRINTED OUT AND
C                           $DMY0M IS SET TO 0.0 FOR !X!>=1.0D12 OR
C                           -1.0D75 FOR !X!<8.797D-77.
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION NAME(2)
      DATA C0,C1,C2,C3,C4,C5,C6,C7,C8,C9
     *   /  1.2249 4472           D-8, -3.9299 65711          D-7,
     *      8.9514 96872 0        D-6, -1.6124 40490 740      D-4,
     *      2.2574 92590 3842     D-3, -2.3703 70267 63209    D-2,
     *      1.7777 77775 58267 6  D-1, -8.8888 88888 61512 4  D-1,
     *      2.6666 66666 66491 50 D 0, -3.9999 99999 99995 63 D 0 /
      DATA D0,D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11
     *   /  2.3353 113            D-9, -8.7149 8135           D-8,
     *      2.3316 55657 5        D-6, -5.0097 83252 95       D-5,
     *      8.5332 26991 090      D-4, -1.1176 62851 97065    D-2,
     *      1.0801 30262 83770 8  D-1, -7.2710 13435 63856 6  D-1,
     *      3.1262 10900 22178 23 D 0, -7.3980 37853 32564 88 D 0,
     *      6.8529 24964 20573 60 D 0,  3.9325 61521 02602 95 D-1 /
      DATA E0,E1,E2,E3,E4,E5,E6,E7,E8,E9,E10,E11,E12,E13,E14,E15
     *   /  4.8468 2              D-12, -2.2800 90             D-11,
     *      6.3132 56             D-11, -5.2574 110            D-10,
     *      7.8189 4200           D- 9, -5.9240 4185           D- 9,
     *     -5.9222 98307 9        D- 7, -6.1872 86946 2        D- 7,
     *      5.2495 83490 136      D- 5, -4.1370 37951 76       D- 6,
     *     -2.3535 74442 27453    D- 3,  7.6458 54114 1893     D- 4,
     *      5.3923 88944 88126 4  D- 2, -3.7165 28343 50646 6  D- 2,
     *     -3.3809 02539 27279 20 D- 1,  1.4786 31433 91226 85 D- 1 /
      DATA F0,F1,F2,F3,F4,F5,F6,F7,F8,F9,F10,F11,F12,F13,F14
     *   /  1.7873 7              D-12,  1.6011 35            D-11,
     *     -3.8677 585            D-10, -2.3735 6163          D- 9,
     *      5.4322 38759          D- 8,  2.6412 21207 5       D- 7,
     *     -5.3968 07752 06       D- 6, -1.7970 21894 410     D- 5,
     *      3.3358 85263 4720     D- 4,  6.7436 67691 6589    D- 4,
     *     -1.1206 70091 19909 2  D- 2, -9.4701 49780 74827   D- 3,
     *      1.4701 02892 64358 67 D- 1,  1.7288 43275 05313 8 D- 2,
     *     -3.0266 72370 24184 87 D- 1                             /
      DATA A0,A1,A2,A3,A4,A5,A6,A7,A8,A9,A10
     *   / -3.0207 732            D-10,  2.0806 2912        D-9,
     *     -7.0003 5699           D- 9,  1.6557 68548       D-8,
     *     -3.5582 32648          D- 8,  8.9550 40706       D-8,
     *     -3.2722 78326 5        D- 7,  2.0593 18935 93    D-6,
     *     -2.8088 72167 333      D- 5,  1.4609 70265 11633 D-3,
     *      7.9788 45608 02865 46 D- 1                          /
      DATA B0,B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12
     *   /  9.4650 569            D-10, -7.0536 9184        D-9,
     *      2.4646 12368          D- 8, -5.5142 32274       D-8,
     *      9.3432 85420          D- 8, -1.4024 48404 3     D-7,
     *      2.2297 66956 2        D- 7, -4.4752 03358 4     D-7,
     *      1.2957 31867 99       D- 6, -6.0676 75972 38    D-6,
     *      5.4070 79074 732      D- 5, -1.2783 48981 98997 D-3,
     *      2.9920 67103 01074 48 D- 1                          /
      DATA PI2,PI4
     *   /  0.63661 97723 67581 34D0, 0.78539 81633 97448 31D0 /
      DATA NAME / '$DMY', '1M  ' /
      AX=ABS(X)
      IF (AX .LT. 8.797D-77) GO TO 60
      IF (AX .GE. 1.0D12) GO TO 50
      IF (AX .GE. 8.0D0) GO TO 30
      IF (AX .GE. 6.0D0) GO TO 20
      IF (AX .GE. 4.0D0) GO TO 10
C
C     0 < !X! < 4
C
      Y=X*0.25D0
      Z=Y*Y
      BJ1=((((((((((C0*Z +C1)*Z +C2)*Z +C3)*Z +C4)*Z +C5)*Z +C6)*Z
     *       +C7)*Z +C8)*Z +C9)*Z +2.0D0) *Y
      ADD=((((((((((((D0*Z +D1)*Z +D2)*Z +D3)*Z +D4)*Z +D5)*Z +D6)*Z
     *       +D7)*Z +D8)*Z +D9)*Z +D10)*Z +D11)*Z -PI2) /X
      BESL21=PI2*BJ1*LOG(AX*0.5D0) +ADD
      GO TO 9999
C
C     4 <= !X! < 6
C
   10 Z=AX-5.0D0
      BESL21=((((((((((((((E0*Z +E1)*Z +E2)*Z +E3)*Z +E4)*Z +E5)*Z
     *       +E6)*Z +E7)*Z +E8)*Z +E9)*Z +E10)*Z +E11)*Z +E12)*Z
     *       +E13)*Z +E14)*Z +E15
      GO TO 40
C
C     6 <= !X! < 8
C
   20 Z=AX-7.0D0
      BESL21=(((((((((((((F0*Z +F1)*Z +F2)*Z +F3)*Z +F4)*Z +F5)*Z
     *       +F6)*Z +F7)*Z +F8)*Z +F9)*Z +F10)*Z +F11)*Z +F12)*Z
     *       +F13)*Z +F14
      GO TO 40
C
C     8 <= !X! < 1.0D12
C
   30 Z=(8.0D0/AX)**2
      A=(((((((((A0*Z +A1)*Z +A2)*Z +A3)*Z +A4)*Z +A5)*Z +A6)*Z +A7)*Z
     *       +A8)*Z +A9)*Z +A10
      B=((((((((((((B0*Z +B1)*Z +B2)*Z +B3)*Z +B4)*Z +B5)*Z +B6)*Z
     *       +B7)*Z +B8)*Z +B9)*Z +B10)*Z +B11)*Z +B12) /AX
      Q=AX-PI4
      BESL21=(B*SIN(Q) - A*COS(Q)) /SQRT(AX)
   40 IF (X .LT. 0.0D0) BESL21=-BESL21
      GO TO 9999
C
C     PARAMETER ERROR
C
   50 BESL21=0.0D0
      GO TO 70
C
   60 BESL21=-1.0D75
   70 IER=2000
C
 9999 RETURN
      END
