1      DOUBLE PRECISION FUNCTION DBESY0 (X)
2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
3c ALL RIGHTS RESERVED.
4c Based on Government Sponsored Research NAS7-03001.
5c>> 1996-03-30 DBESY0 Krogh  Added external statement.
6C>> 1995-11-13 DBESY0 Krogh  Changes to simplify conversion to C.
7C>> 1994-11-11 DBESY0 Krogh   Declared all vars.
8C>> 1994-10-20 DBESY0 Krogh  Changes to use M77CON
9C>> 1990-11-29 DBESY0 CLL
10C>> 1985-08-02 DBESY0 Lawson  Initial code.
11C JULY 1977 EDITION.  W. FULLERTON, C3, LOS ALAMOS SCIENTIFIC LAB.
12C C.L.LAWSON & S.CHAN, JPL, 1984 FEB ADAPTED TO JPL MATH77 LIBRARY.
13c     ------------------------------------------------------------------
14c--D replaces "?": ?BESY0, ?BESJ0, ?BMP0, ?INITS, ?CSEVL, ?ERM1
15c     ------------------------------------------------------------------
16      EXTERNAL D1MACH, DBESJ0, DCSEVL
17      INTEGER NTY0
18      DOUBLE PRECISION X, BY0CS(19), AMPL, THETA, TWODPI, XSML,
19     1  Y, D1MACH, DCSEVL, DBESJ0
20C
21C SERIES FOR BY0        ON THE INTERVAL  0.          TO  1.60000D+01
22C                                        WITH WEIGHTED ERROR   8.14D-32
23C                                         LOG WEIGHTED ERROR  31.09
24C                               SIGNIFICANT FIGURES REQUIRED  30.31
25C                                    DECIMAL PLACES REQUIRED  31.73
26C
27      SAVE NTY0, XSML
28C
29      DATA BY0CS / -.1127783939286557321793980546028D-1,
30     *  -.1283452375604203460480884531838D+0,
31     *  -.1043788479979424936581762276618D+0,
32     *  +.2366274918396969540924159264613D-1,
33     *  -.2090391647700486239196223950342D-2,
34     *  +.1039754539390572520999246576381D-3,
35     *  -.3369747162423972096718775345037D-5,
36     *  +.7729384267670667158521367216371D-7,
37     *  -.1324976772664259591443476068964D-8,
38     *  +.1764823261540452792100389363158D-10,
39     *  -.1881055071580196200602823012069D-12,
40     *  +.1641865485366149502792237185749D-14,
41     *  -.1195659438604606085745991006720D-16,
42     *  +.7377296297440185842494112426666D-19,
43     *  -.3906843476710437330740906666666D-21,
44     *  +.1795503664436157949829120000000D-23,
45     *  -.7229627125448010478933333333333D-26,
46     *  +.2571727931635168597333333333333D-28,
47     *  -.8141268814163694933333333333333D-31     /
48C
49      DATA TWODPI / 0.636619772367581343075535053490057D0 /
50      DATA NTY0, XSML / 0, 0.D0 /
51C     ------------------------------------------------------------------
52      IF (NTY0 .eq. 0) then
53      call DINITS (BY0CS, 19, 0.1D0*D1MACH(3), NTY0)
54      XSML = SQRT (4.0D0*D1MACH(3))
55      endif
56C
57      IF (X .LE. 0.D0) THEN
58        DBESY0 = 0.D0
59        CALL DERM1 ('DBESY0',1,0,'X IS ZERO OR NEGATIVE','X',X,'.')
60      ELSE IF (X .LE. 4.D0) THEN
61        IF (X .LE. XSML) THEN
62          Y = 0.D0
63        ELSE
64          Y = X * X
65        END IF
66        DBESY0 = TWODPI*LOG(0.5D0*X)*DBESJ0(X) + .375D0 +
67     *            DCSEVL (.125D0*Y-1.D0, BY0CS, NTY0)
68      ELSE
69        CALL DBMP0 (X, AMPL, THETA)
70        DBESY0 = AMPL * SIN(THETA)
71      END IF
72C
73      RETURN
74C
75      END
76