1*DECK DBESJ1
2      DOUBLE PRECISION FUNCTION DBESJ1 (X)
3C***BEGIN PROLOGUE  DBESJ1
4C***PURPOSE  Compute the Bessel function of the first kind of order one.
5C***LIBRARY   SLATEC (FNLIB)
6C***CATEGORY  C10A1
7C***TYPE      DOUBLE PRECISION (BESJ1-S, DBESJ1-D)
8C***KEYWORDS  BESSEL FUNCTION, FIRST KIND, FNLIB, ORDER ONE,
9C             SPECIAL FUNCTIONS
10C***AUTHOR  Fullerton, W., (LANL)
11C***DESCRIPTION
12C
13C DBESJ1(X) calculates the double precision Bessel function of the
14C first kind of order one for double precision argument X.
15C
16C Series for BJ1        on the interval  0.          to  1.60000E+01
17C                                        with weighted error   1.16E-33
18C                                         log weighted error  32.93
19C                               significant figures required  32.36
20C                                    decimal places required  33.57
21C
22C***REFERENCES  (NONE)
23C***ROUTINES CALLED  D1MACH, D9B1MP, DCSEVL, INITDS, XERMSG
24C***REVISION HISTORY  (YYMMDD)
25C   780601  DATE WRITTEN
26C   890531  Changed all specific intrinsics to generic.  (WRB)
27C   890531  REVISION DATE from Version 3.2
28C   891214  Prologue converted to Version 4.0 format.  (BAB)
29C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
30C   910401  Corrected error in code which caused values to have the
31C           wrong sign for arguments less than 4.0.  (WRB)
32C***END PROLOGUE  DBESJ1
33      DOUBLE PRECISION X, BJ1CS(19), AMPL, THETA, XSML, XMIN, Y,
34     1  D1MACH, DCSEVL
35      LOGICAL FIRST
36      SAVE BJ1CS, NTJ1, XSML, XMIN, FIRST
37      DATA BJ1CS(  1) / -.1172614151 3332786560 6240574524 003 D+0    /
38      DATA BJ1CS(  2) / -.2536152183 0790639562 3030884554 698 D+0    /
39      DATA BJ1CS(  3) / +.5012708098 4469568505 3656363203 743 D-1    /
40      DATA BJ1CS(  4) / -.4631514809 6250819184 2619728789 772 D-2    /
41      DATA BJ1CS(  5) / +.2479962294 1591402453 9124064592 364 D-3    /
42      DATA BJ1CS(  6) / -.8678948686 2788258452 1246435176 416 D-5    /
43      DATA BJ1CS(  7) / +.2142939171 4379369150 2766250991 292 D-6    /
44      DATA BJ1CS(  8) / -.3936093079 1831797922 9322764073 061 D-8    /
45      DATA BJ1CS(  9) / +.5591182317 9468800401 8248059864 032 D-10   /
46      DATA BJ1CS( 10) / -.6327616404 6613930247 7695274014 880 D-12   /
47      DATA BJ1CS( 11) / +.5840991610 8572470032 6945563268 266 D-14   /
48      DATA BJ1CS( 12) / -.4482533818 7012581903 9135059199 999 D-16   /
49      DATA BJ1CS( 13) / +.2905384492 6250246630 6018688000 000 D-18   /
50      DATA BJ1CS( 14) / -.1611732197 8414416541 2118186666 666 D-20   /
51      DATA BJ1CS( 15) / +.7739478819 3927463729 8346666666 666 D-23   /
52      DATA BJ1CS( 16) / -.3248693782 1119984114 3466666666 666 D-25   /
53      DATA BJ1CS( 17) / +.1202237677 2274102272 0000000000 000 D-27   /
54      DATA BJ1CS( 18) / -.3952012212 6513493333 3333333333 333 D-30   /
55      DATA BJ1CS( 19) / +.1161678082 2664533333 3333333333 333 D-32   /
56      DATA FIRST /.TRUE./
57C***FIRST EXECUTABLE STATEMENT  DBESJ1
58      IF (FIRST) THEN
59         NTJ1 = INITDS (BJ1CS, 19, 0.1*REAL(D1MACH(3)))
60C
61         XSML = SQRT(8.0D0*D1MACH(3))
62         XMIN = 2.0D0*D1MACH(1)
63      ENDIF
64      FIRST = .FALSE.
65C
66      Y = ABS(X)
67      IF (Y.GT.4.0D0) GO TO 20
68C
69      DBESJ1 = 0.0D0
70      IF (Y.EQ.0.0D0) RETURN
71      IF (Y .LE. XMIN) CALL XERMSG ('SLATEC', 'DBESJ1',
72     +   'ABS(X) SO SMALL J1 UNDERFLOWS', 1, 1)
73      IF (Y.GT.XMIN) DBESJ1 = 0.5D0*X
74      IF (Y.GT.XSML) DBESJ1 = X*(.25D0 + DCSEVL (.125D0*Y*Y-1.D0,
75     1  BJ1CS, NTJ1) )
76      RETURN
77C
78 20   CALL D9B1MP (Y, AMPL, THETA)
79      DBESJ1 = SIGN (AMPL, X) * COS(THETA)
80C
81      RETURN
82      END
83