C
      REAL FUNCTION BESL20*8 (X)
C
C   FUNCTION            - BY CHEBYSHEV INTERPOLATION, WE OBTAIN THE
C                           VALUE OF THE BESSEL FUNCTION OF THE SECOND
C                           KIND OF ORDER ZERO.
C
C   NOTICE              - WHEN !X! DOES NOT SATISFY 0<!X!<1.0D12, ERROR
C                           MESSAGE IS PRINTED OUT AND BESL20 IS SET TO
C                           0.0 FOR !X!>=1.0D12 OR -1.0D75 FOR X=0.0.
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION NAME(2)
      DATA C0,C1,C2,C3,C4,C5,C6,C7,C8,C9,C10
     *   / -2.2270 614            D-9,  7.8594 0046           D-8,
     *     -1.9891 96578 6        D-6,  4.0310 96038 59       D-5,
     *     -6.4499 78225 950      D-4,  7.9012 34182 8039     D-3,
     *     -7.1111 11100 46304    D-2,  4.4444 44444 25769 5  D-1,
     *     -1.7777 77777 77584 46 D 0,  3.9999 99999 99989 66 D 0,
     *     -3.9999 99999 99999 78 D 0                             /
      DATA D0,D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11
     *   /  3.4432 403            D-9, -1.1755 92825          D-7,
     *      2.8512 73775 9        D-6, -5.4934 40736 91       D-5,
     *      8.2765 91200 190      D-4, -9.4202 58441 9331     D-3,
     *      7.7237 20812 88392    D-2, -4.2614 41276 30476 6  D-1,
     *      1.4216 34389 58614 94 D 0, -2.3498 51013 41885 87 D 0,
     *      1.0766 11468 68385 81 D 0,  3.6746 69051 96615 97 D-1 /
      DATA E0,E1,E2,E3,E4,E5,E6,E7,E8,E9,E10,E11,E12,E13,E14,E15
     *   /  1.5041 0              D-12, -5.7211 7              D-12,
     *      4.0501 46             D-11, -6.4960 948            D-10,
     *      5.3845 950            D-10,  5.9221 31699          D- 8,
     *      6.8747 70130          D- 8, -6.5619 78581 69       D- 6,
     *      5.9100 53936 3        D- 7,  3.9226 24068 4696     D- 4,
     *     -1.5291 70822 7790     D- 4, -1.3480 97236 21783 1  D- 2,
     *      1.2388 42781 16876 7  D- 2,  1.6904 51269 63638 42 D- 1,
     *     -1.4786 31433 91226 86 D- 1, -3.0851 76252 49033 78 D- 1 /
      DATA F0,F1,F2,F3,F4,F5,F6,F7,F8,F9,F10,F11,F12,F13,F14
     *   / -1.1446 0              D-12,  2.9305 15             D-11,
     *      1.9780 007            D-10, -4.9377 2861           D- 9,
     *     -2.6412 21656          D- 8,  5.9964 47937 8        D- 7,
     *      2.2462 77371 07       D- 6, -4.7655 50355 443      D- 5,
     *     -1.1239 44615 2872     D- 4,  2.2413 40182 35420    D- 3,
     *      2.3675 37445 18725    D- 3, -4.9003 42975 47821 5  D- 2,
     *     -8.6442 16375 26570    D- 3,  3.0266 72370 24184 76 D- 1,
     *     -2.5949 74396 72092 65 D- 2                              /
      DATA A0,A1,A2,A3,A4,A5,A6,A7,A8,A9,A10
     *   /  2.8228 605            D-10, -1.9415 5134       D-9,
     *      6.5160 2613           D- 9, -1.5339 82923      D-8,
     *      3.2676 79187          D- 8, -8.0992 84883      D-8,
     *      2.8872 15007 2        D- 7, -1.7424 98968 02   D-6,
     *      2.1846 78334 973      D- 5, -8.7658 21590 6110 D-4,
     *      7.9788 45608 02865 26 D- 1                         /
      DATA B0,B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12
     *   / -8.9331 052           D-10,  6.6532 7442       D-9,
     *     -2.3224 79945         D- 8,  5.1878 39807      D-8,
     *     -8.7651 63373         D- 8,  1.3089 75220 0    D-7,
     *     -2.0632 29204 9       D- 7,  4.0856 34631 8    D-7,
     *     -1.1593 30628 42      D- 6,  5.2586 51475 08   D-6,
     *     -4.4239 73781 374     D- 5,  9.1310 64157 0497 D-4,
     *     -9.9735 57010 03581 5 D- 2                         /
      DATA PI2,PI4
     *   /  0.63661 97723 67581 34D0, 0.78539 81633 97448 31D0 /
      DATA NAME / '$DMY', '0M  ' /
      AX=ABS(X)
      IF (AX .EQ. 0.0D0) GO TO  50
      IF (AX .GE. 1.0D12) GO TO 40
      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
      Z=(AX*0.25D0)**2
      BJ0=((((((((((C0*Z +C1)*Z +C2)*Z +C3)*Z +C4)*Z +C5)*Z +C6)*Z
     *    +C7)*Z +C8)*Z +C9)*Z +C10)*Z +1.0D0
      ADD=((((((((((D0*Z +D1)*Z +D2)*Z +D3)*Z +D4)*Z +D5)*Z +D6)*Z
     *    +D7)*Z +D8)*Z +D9)*Z +D10)*Z +D11
      BESL20=PI2*BJ0*LOG(AX*0.5D0) +ADD
      GO TO 9999
C
C     4 <= !X! < 6
C
   10 Z=AX-5.0D0
      BESL20=((((((((((((((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 9999
C
C     6 <= !X! < 8
C
   20 Z=AX-7.0D0
      BESL20=(((((((((((((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 9999
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
      BESL20=(A*SIN(Q) + B*COS(Q)) /SQRT(AX)
      GO TO 9999
C
C     PARAMETER ERROR
C
   40 BESL20=0.0D0
      GO TO 60
C
   50 BESL20=-1.0D75
   60 IER=2000
C
 9999 RETURN
      END
