DOUBLE PRECISION FUNCTION DBESY1 (X)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
c>> 1996-03-30 DBESY1 Krogh  Added external statement.
C>> 1995-11-03 DBESY1 Krogh  Changes to simplify C conversion.
C>> 1994-11-11 DBESY1 Krogh  Declared all vars.
C>> 1994-10-20 DBESY1 Krogh  Changes to use M77CON
C>> 1991-01-14 DBESY1 CLL Changed to generic name SIN.
C>> 1990-11-29 CLL
C>> 1985-08-02 DBESY1 Lawson  Initial code.
C JULY 1977 EDITION.  W. FULLERTON, C3, LOS ALAMOS SCIENTIFIC LAB.
C C.L.LAWSON & S.CHAN, JPL, 1984 FEB ADAPTED TO JPL MATH77 LIBRARY.
c     ------------------------------------------------------------------
c--D replaces "?": ?BESY1, ?BESJ1, ?BMP1, ?INITS, ?CSEVL, ?ERM1
c     ------------------------------------------------------------------
      EXTERNAL D1MACH, DBESJ1, DCSEVL
      INTEGER NTY1
      DOUBLE PRECISION X, BY1CS(20), AMPL, THETA, TWODPI, XMIN, XSML,
     1  Y, D1MACH, DCSEVL, DBESJ1
C
C SERIES FOR BY1        ON THE INTERVAL  0.          TO  1.60000D+01
C                                        WITH WEIGHTED ERROR   8.65D-33
C                                         LOG WEIGHTED ERROR  32.06
C                               SIGNIFICANT FIGURES REQUIRED  32.17
C                                    DECIMAL PLACES REQUIRED  32.71
C
      SAVE NTY1, XMIN, XSML
C
      DATA BY1CS / +.320804710061190862932352018628015D-1,
     *  +.126270789743350044953431725999727D+1,
     *  +.649996189992317500097490637314144D-2,
     *  -.893616452886050411653144160009712D-1,
     *  +.132508812217570954512375510370043D-1,
     *  -.897905911964835237753039508298105D-3,
     *  +.364736148795830678242287368165349D-4,
     *  -.100137438166600055549075523845295D-5,
     *  +.199453965739017397031159372421243D-7,
     *  -.302306560180338167284799332520743D-9,
     *  +.360987815694781196116252914242474D-11,
     *  -.348748829728758242414552947409066D-13,
     *  +.278387897155917665813507698517333D-15,
     *  -.186787096861948768766825352533333D-17,
     *  +.106853153391168259757070336000000D-19,
     *  -.527472195668448228943872000000000D-22,
     *  +.227019940315566414370133333333333D-24,
     *  -.859539035394523108693333333333333D-27,
     *  +.288540437983379456000000000000000D-29,
     *  -.864754113893717333333333333333333D-32 /
C
      DATA TWODPI / 0.636619772367581343075535053490057D0 /
      DATA NTY1, XMIN, XSML / 0, 2*0.D0 /
C     ------------------------------------------------------------------
      IF (NTY1 .EQ. 0) then
      call DINITS (BY1CS, 20, 0.1D0*D1MACH(3), NTY1)
C
      XMIN = 1.571D0 * EXP (MAX(LOG(D1MACH(1)), -LOG(D1MACH(2))) +
     1  0.01D0)
      XSML = SQRT (4.0D0*D1MACH(3))
      endif
C
      IF (X .LE. 0.D0) THEN
        DBESY1 = 0.D0
        CALL DERM1 ('DBESY1',1,0,'X IS ZERO OR NEGATIVE','X',X,'.')
      ELSE IF (X .LE. XMIN) THEN
        DBESY1 = 0.D0
        CALL DERM1 ('DBESY1',2,0,'X SO SMALL Y1 OVERFLOWS','X',X,'.')
      ELSE IF (X .LE. 4.D0) THEN
        IF (X .LE. XSML) THEN
          Y = 0.D0
        ELSE
          Y = X * X
        END IF
        DBESY1 = TWODPI * LOG(0.5D0*X)*DBESJ1(X) + (0.5D0 +
     *           DCSEVL (.125D0*Y-1.D0, BY1CS, NTY1))/X
      ELSE
        CALL DBMP1 (X, AMPL, THETA)
        DBESY1 = AMPL * SIN(THETA)
      END IF
      RETURN
      END

.