1      DOUBLE PRECISION FUNCTION DBESY1(X)
2C***BEGIN PROLOGUE  DBESY1
3C***DATE WRITTEN   770701   (YYMMDD)
4C***REVISION DATE  820801   (YYMMDD)
5C***CATEGORY NO.  C10A1
6C***KEYWORDS  BESSEL FUNCTION,DOUBLE PRECISION,ORDER ONE,SECOND KIND,
7C             SPECIAL FUNCTION
8C***AUTHOR  FULLERTON, W., (LANL)
9C***PURPOSE  Computes the d.p. Bessel function of the second kind of
10C            order one.
11C***DESCRIPTION
12C
13C DBESY1(X) calculates the double precision Bessel function of the
14C second kind of order for double precision argument X.
15C
16C Series for BY1        on the interval  0.          to  1.60000E+01
17C                                        with weighted error   8.65E-33
18C                                         log weighted error  32.06
19C                               significant figures required  32.17
20C                                    decimal places required  32.71
21C***REFERENCES  (NONE)
22C***ROUTINES CALLED  D1MACH,D9B1MP,DBESJ1,DCSEVL,INITDS,XERROR
23C***END PROLOGUE  DBESY1
24      DOUBLE PRECISION X, BY1CS(20), AMPL, THETA, TWODPI, XMIN, XSML,
25     1  Y, Z, D1MACH, DCSEVL, DBESJ1
26      DATA BY1 CS(  1) / +.3208047100 6119086293 2352018628 015 D-1    /
27      DATA BY1 CS(  2) / +.1262707897 4335004495 3431725999 727 D+1    /
28      DATA BY1 CS(  3) / +.6499961899 9231750009 7490637314 144 D-2    /
29      DATA BY1 CS(  4) / -.8936164528 8605041165 3144160009 712 D-1    /
30      DATA BY1 CS(  5) / +.1325088122 1757095451 2375510370 043 D-1    /
31      DATA BY1 CS(  6) / -.8979059119 6483523775 3039508298 105 D-3    /
32      DATA BY1 CS(  7) / +.3647361487 9583067824 2287368165 349 D-4    /
33      DATA BY1 CS(  8) / -.1001374381 6660005554 9075523845 295 D-5    /
34      DATA BY1 CS(  9) / +.1994539657 3901739703 1159372421 243 D-7    /
35      DATA BY1 CS( 10) / -.3023065601 8033816728 4799332520 743 D-9    /
36      DATA BY1 CS( 11) / +.3609878156 9478119611 6252914242 474 D-11   /
37      DATA BY1 CS( 12) / -.3487488297 2875824241 4552947409 066 D-13   /
38      DATA BY1 CS( 13) / +.2783878971 5591766581 3507698517 333 D-15   /
39      DATA BY1 CS( 14) / -.1867870968 6194876876 6825352533 333 D-17   /
40      DATA BY1 CS( 15) / +.1068531533 9116825975 7070336000 000 D-19   /
41      DATA BY1 CS( 16) / -.5274721956 6844822894 3872000000 000 D-22   /
42      DATA BY1 CS( 17) / +.2270199403 1556641437 0133333333 333 D-24   /
43      DATA BY1 CS( 18) / -.8595390353 9452310869 3333333333 333 D-27   /
44      DATA BY1 CS( 19) / +.2885404379 8337945600 0000000000 000 D-29   /
45      DATA BY1 CS( 20) / -.8647541138 9371733333 3333333333 333 D-32   /
46      DATA TWODPI / 0.6366197723 6758134307 5535053490 057 D0 /
47      DATA NTY1, XMIN, XSML / 0, 2*0.D0 /
48C***FIRST EXECUTABLE STATEMENT  DBESY1
49      IF (NTY1.NE.0) GO TO 10
50      NTY1 = INITDS (BY1CS, 20, 0.1*SNGL(D1MACH(3)))
51C
52      XMIN = 1.571D0 * DEXP (DMAX1(DLOG(D1MACH(1)), -DLOG(D1MACH(2))) +
53     1  0.01D0)
54      XSML = DSQRT (4.0D0*D1MACH(3))
55C
56 10   IF (X.LE.0.D0) CALL XERROR ( 'DBESY1  X IS ZERO OR NEGATIVE', 29,
57     1  1, 2)
58      IF (X.GT.4.0D0) GO TO 20
59C
60      IF (X.LT.XMIN) CALL XERROR ( 'DBESY1  X SO SMALL Y1 OVERFLOWS',
61     1  31, 3, 2)
62      Y = 0.D0
63      IF (X.GT.XSML) Y = X*X
64      DBESY1 = TWODPI * DLOG(0.5D0*X)*DBESJ1(X) + (0.5D0 +
65     1  DCSEVL (.125D0*Y-1.D0, BY1CS, NTY1))/X
66      RETURN
67C
68 20   CALL D9B1MP (X, AMPL, THETA)
69      DBESY1 = AMPL * DSIN(THETA)
70      RETURN
71C
72      END
73