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