C
      REAL FUNCTION BESJ1*8 (X)
C
C   NAME                - BESJ1 : DOUBLE PRECISION
C
C   USAGE               - Y=BESJ1(X)
C
C   FUNCTION            - BY CHEBYSHEV INTERPOLATION, WE OBTAIN THE
C                           VALUE OF THE BESSEL FUNCTION OF THE FIRST
C                           KIND OF ORDER ONE.
C
C   ARGUMENTS  X        - INPUT. VARIABLE OF THE BESSEL FUNCTION
C                           (!X!<1.0D12).
C              (BESJ1) - OUTPUT. THE VALUE OF THE BESSEL FUNCTION.
C
C   NOTICE              - WHEN !X!>=1.0D12, ERROR MESSAGE IS PRINTED
C                           OUT AND BESJ1 IS SET TO 0.0.
C
C   REQD. MSLII SUBPROG.- $SUERM
C
C   STATUS              - S-1511-1 05-02
C
C   HISTORY             - DATE.    1981.10
C                                  1986. 6
C
C-----------------------------------------------------------------------
C
      IMPLICIT REAL*8 (A-H,O-Z)
C      GENERIC
      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,D12,D13,D14
     *   /  2.2193 1              D-12,  4.1163 4              D-12,
     *     -4.3816 515            D-10, -5.7599 652            D-10,
     *      6.2306 13774          D- 8,  4.8857 13666          D- 8,
     *     -6.1029 92681 89       D- 6, -1.7066 40506 78       D- 6,
     *      3.7656 50660 9458     D- 4, -6.7213 28811 703      D- 5,
     *     -1.2601 05985 12045 4  D- 2,  6.8295 50387 40438    D- 3,
     *      1.6844 60804 23507 80 D- 1, -1.1208 09437 96045 28 D- 1,
     *     -3.2757 91375 91465 22 D- 1                              /
      DATA E0,E1,E2,E3,E4,E5,E6,E7,E8,E9,E10,E11,E12,E13,E14
     *   / -1.0746 3             D-12,  2.7780 03             D-11,
     *      2.0080 677           D-10, -4.7578 8198           D- 9,
     *     -2.6460 77937         D- 8,  5.7071 42412 9        D- 7,
     *      2.3218 13149 97      D- 6, -4.5285 96911 743      D- 5,
     *     -1.2058 28825 8765    D- 4,  2.1170 94601 37864    D- 3,
     *      2.9657 44869 54234   D- 3, -4.7160 51742 84081 7  D- 2,
     *     -1.9188 38969 35370 8 D- 2,  3.0074 82453 02747 75 D- 1,
     *     -4.6828 23482 34583   D- 3                              /
      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 PI4  / 0.78539 81633 97448 31D0 /
C      DATA NAME / '$DMJ', '1M  ' /
      AX=ABS(X)
      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     !X! < 4
C
      Y=X*0.25D0
      Z=Y*Y
      BESJ1=((((((((((C0*Z +C1)*Z +C2)*Z +C3)*Z +C4)*Z +C5)*Z +C6)*Z
     *       +C7)*Z +C8)*Z +C9)*Z +2.0D0) *Y
      GO TO 9999
C
C     4 <= !X! < 6
C
   10 Z=AX-5.0D0
      BESJ1=(((((((((((((D0*Z +D1)*Z +D2)*Z +D3)*Z +D4)*Z +D5)*Z
     *       +D6)*Z +D7)*Z +D8)*Z +D9)*Z +D10)*Z +D11)*Z +D12)*Z
     *       +D13)*Z +D14
      GO TO 40
C
C     6 <= !X! < 8
C
   20 Z=AX-7.0D0
      BESJ1=(((((((((((((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
      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
      BESJ1=(A*SIN(Q) + B*COS(Q)) /SQRT(AX)
   40 IF (X .LT. 0.0D0) BESJ1=-BESJ1
      GO TO 9999
C
C     PARAMETER ERROR
C
   50 BESJ1=0.0D0
C      CALL $SUERM (NAME,2000,IER)
C
 9999 RETURN
      END
