1      REAL             FUNCTION SBESY1 (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 SBESY1 Krogh  Added external statement.
6C>> 1995-11-03 SBESY1 Krogh  Changes to simplify C conversion.
7C>> 1994-11-11 SBESY1 Krogh  Declared all vars.
8C>> 1994-10-20 SBESY1 Krogh  Changes to use M77CON
9C>> 1991-01-14 SBESY1 CLL Changed to generic name SIN.
10C>> 1990-11-29 CLL
11C>> 1985-08-02 SBESY1 Lawson  Initial code.
12C JULY 1977 EDITION.  W. FULLERTON, C3, LOS ALAMOS SCIENTIFIC LAB.
13C C.L.LAWSON & S.CHAN, JPL, 1984 FEB ADAPTED TO JPL MATH77 LIBRARY.
14c     ------------------------------------------------------------------
15c--S replaces "?": ?BESY1, ?BESJ1, ?BMP1, ?INITS, ?CSEVL, ?ERM1
16c     ------------------------------------------------------------------
17      EXTERNAL R1MACH, SBESJ1, SCSEVL
18      INTEGER NTY1
19      REAL             X, BY1CS(20), AMPL, THETA, TWODPI, XMIN, XSML,
20     1  Y, R1MACH, SCSEVL, SBESJ1
21C
22C SERIES FOR BY1        ON THE INTERVAL  0.          TO  1.60000E+01
23C                                        WITH WEIGHTED ERROR   8.65E-33
24C                                         LOG WEIGHTED ERROR  32.06
25C                               SIGNIFICANT FIGURES REQUIRED  32.17
26C                                    DECIMAL PLACES REQUIRED  32.71
27C
28      SAVE NTY1, XMIN, XSML
29C
30      DATA BY1CS / +.320804710061190862932352018628015E-1,
31     *  +.126270789743350044953431725999727E+1,
32     *  +.649996189992317500097490637314144E-2,
33     *  -.893616452886050411653144160009712E-1,
34     *  +.132508812217570954512375510370043E-1,
35     *  -.897905911964835237753039508298105E-3,
36     *  +.364736148795830678242287368165349E-4,
37     *  -.100137438166600055549075523845295E-5,
38     *  +.199453965739017397031159372421243E-7,
39     *  -.302306560180338167284799332520743E-9,
40     *  +.360987815694781196116252914242474E-11,
41     *  -.348748829728758242414552947409066E-13,
42     *  +.278387897155917665813507698517333E-15,
43     *  -.186787096861948768766825352533333E-17,
44     *  +.106853153391168259757070336000000E-19,
45     *  -.527472195668448228943872000000000E-22,
46     *  +.227019940315566414370133333333333E-24,
47     *  -.859539035394523108693333333333333E-27,
48     *  +.288540437983379456000000000000000E-29,
49     *  -.864754113893717333333333333333333E-32 /
50C
51      DATA TWODPI / 0.636619772367581343075535053490057E0 /
52      DATA NTY1, XMIN, XSML / 0, 2*0.E0 /
53C     ------------------------------------------------------------------
54      IF (NTY1 .EQ. 0) then
55      call SINITS (BY1CS, 20, 0.1E0*R1MACH(3), NTY1)
56C
57      XMIN = 1.571E0 * EXP (MAX(LOG(R1MACH(1)), -LOG(R1MACH(2))) +
58     1  0.01E0)
59      XSML = SQRT (4.0E0*R1MACH(3))
60      endif
61C
62      IF (X .LE. 0.E0) THEN
63        SBESY1 = 0.E0
64        CALL SERM1 ('SBESY1',1,0,'X IS ZERO OR NEGATIVE','X',X,'.')
65      ELSE IF (X .LE. XMIN) THEN
66        SBESY1 = 0.E0
67        CALL SERM1 ('SBESY1',2,0,'X SO SMALL Y1 OVERFLOWS','X',X,'.')
68      ELSE IF (X .LE. 4.E0) THEN
69        IF (X .LE. XSML) THEN
70          Y = 0.E0
71        ELSE
72          Y = X * X
73        END IF
74        SBESY1 = TWODPI * LOG(0.5E0*X)*SBESJ1(X) + (0.5E0 +
75     *           SCSEVL (.125E0*Y-1.E0, BY1CS, NTY1))/X
76      ELSE
77        CALL SBMP1 (X, AMPL, THETA)
78        SBESY1 = AMPL * SIN(THETA)
79      END IF
80      RETURN
81      END
82