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