1*DECK DBESY1
2      DOUBLE PRECISION FUNCTION DBESY1 (X)
3C***BEGIN PROLOGUE  DBESY1
4C***PURPOSE  Compute the Bessel function of the second kind of order
5C            one.
6C***LIBRARY   SLATEC (FNLIB)
7C***CATEGORY  C10A1
8C***TYPE      DOUBLE PRECISION (BESY1-S, DBESY1-D)
9C***KEYWORDS  BESSEL FUNCTION, FNLIB, ORDER ONE, SECOND KIND,
10C             SPECIAL FUNCTIONS
11C***AUTHOR  Fullerton, W., (LANL)
12C***DESCRIPTION
13C
14C DBESY1(X) calculates the double precision Bessel function of the
15C second kind of order for double precision argument X.
16C
17C Series for BY1        on the interval  0.          to  1.60000E+01
18C                                        with weighted error   8.65E-33
19C                                         log weighted error  32.06
20C                               significant figures required  32.17
21C                                    decimal places required  32.71
22C
23C***REFERENCES  (NONE)
24C***ROUTINES CALLED  D1MACH, D9B1MP, DBESJ1, DCSEVL, INITDS, XERMSG
25C***REVISION HISTORY  (YYMMDD)
26C   770701  DATE WRITTEN
27C   890531  Changed all specific intrinsics to generic.  (WRB)
28C   890531  REVISION DATE from Version 3.2
29C   891214  Prologue converted to Version 4.0 format.  (BAB)
30C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
31C***END PROLOGUE  DBESY1
32      DOUBLE PRECISION X, BY1CS(20), AMPL, THETA, TWODPI, XMIN, XSML,
33     1  Y, D1MACH, DCSEVL, DBESJ1
34      LOGICAL FIRST
35      SAVE BY1CS, TWODPI, NTY1, XMIN, XSML, FIRST
36      DATA BY1CS(  1) / +.3208047100 6119086293 2352018628 015 D-1    /
37      DATA BY1CS(  2) / +.1262707897 4335004495 3431725999 727 D+1    /
38      DATA BY1CS(  3) / +.6499961899 9231750009 7490637314 144 D-2    /
39      DATA BY1CS(  4) / -.8936164528 8605041165 3144160009 712 D-1    /
40      DATA BY1CS(  5) / +.1325088122 1757095451 2375510370 043 D-1    /
41      DATA BY1CS(  6) / -.8979059119 6483523775 3039508298 105 D-3    /
42      DATA BY1CS(  7) / +.3647361487 9583067824 2287368165 349 D-4    /
43      DATA BY1CS(  8) / -.1001374381 6660005554 9075523845 295 D-5    /
44      DATA BY1CS(  9) / +.1994539657 3901739703 1159372421 243 D-7    /
45      DATA BY1CS( 10) / -.3023065601 8033816728 4799332520 743 D-9    /
46      DATA BY1CS( 11) / +.3609878156 9478119611 6252914242 474 D-11   /
47      DATA BY1CS( 12) / -.3487488297 2875824241 4552947409 066 D-13   /
48      DATA BY1CS( 13) / +.2783878971 5591766581 3507698517 333 D-15   /
49      DATA BY1CS( 14) / -.1867870968 6194876876 6825352533 333 D-17   /
50      DATA BY1CS( 15) / +.1068531533 9116825975 7070336000 000 D-19   /
51      DATA BY1CS( 16) / -.5274721956 6844822894 3872000000 000 D-22   /
52      DATA BY1CS( 17) / +.2270199403 1556641437 0133333333 333 D-24   /
53      DATA BY1CS( 18) / -.8595390353 9452310869 3333333333 333 D-27   /
54      DATA BY1CS( 19) / +.2885404379 8337945600 0000000000 000 D-29   /
55      DATA BY1CS( 20) / -.8647541138 9371733333 3333333333 333 D-32   /
56      DATA TWODPI / 0.6366197723 6758134307 5535053490 057 D0 /
57      DATA FIRST /.TRUE./
58C***FIRST EXECUTABLE STATEMENT  DBESY1
59      IF (FIRST) THEN
60         NTY1 = INITDS (BY1CS, 20, 0.1*REAL(D1MACH(3)))
61C
62         XMIN = 1.571D0 * EXP (MAX(LOG(D1MACH(1)), -LOG(D1MACH(2))) +
63     1     0.01D0)
64         XSML = SQRT(4.0D0*D1MACH(3))
65      ENDIF
66      FIRST = .FALSE.
67C
68      IF (X .LE. 0.D0) CALL XERMSG ('SLATEC', 'DBESY1',
69     +   'X IS ZERO OR NEGATIVE', 1, 2)
70      IF (X.GT.4.0D0) GO TO 20
71C
72      IF (X .LT. XMIN) CALL XERMSG ('SLATEC', 'DBESY1',
73     +   'X SO SMALL Y1 OVERFLOWS', 3, 2)
74      Y = 0.D0
75      IF (X.GT.XSML) Y = X*X
76      DBESY1 = TWODPI * LOG(0.5D0*X)*DBESJ1(X) + (0.5D0 +
77     1  DCSEVL (.125D0*Y-1.D0, BY1CS, NTY1))/X
78      RETURN
79C
80 20   CALL D9B1MP (X, AMPL, THETA)
81      DBESY1 = AMPL * SIN(THETA)
82      RETURN
83C
84      END
85