1 FUNCTION BESY1(X) 2C***BEGIN PROLOGUE BESY1 3C***DATE WRITTEN 770401 (YYMMDD) 4C***REVISION DATE 820801 (YYMMDD) 5C***CATEGORY NO. C10A1 6C***KEYWORDS BESSEL FUNCTION,ORDER ONE,SECOND KIND,SPECIAL FUNCTION 7C***AUTHOR FULLERTON, W., (LANL) 8C***PURPOSE Computes the Bessel function of the second kind of order 9C one. 10C***DESCRIPTION 11C 12C BESY1(X) calculates the Bessel function of the second kind of 13C order one for real argument X. 14C 15C Series for BY1 on the interval 0. to 1.60000D+01 16C with weighted error 1.87E-18 17C log weighted error 17.73 18C significant figures required 17.83 19C decimal places required 18.30 20C 21C Series for BM1 on the interval 0. to 6.25000D-02 22C with weighted error 5.61E-17 23C log weighted error 16.25 24C significant figures required 14.97 25C decimal places required 16.91 26C 27C Series for BTH1 on the interval 0. to 6.25000D-02 28C with weighted error 4.10E-17 29C log weighted error 16.39 30C significant figures required 15.96 31C decimal places required 17.08 32C***REFERENCES (NONE) 33C***ROUTINES CALLED BESJ1,CSEVL,INITS,R1MACH,XERROR 34C***END PROLOGUE BESY1 35 DIMENSION BY1CS(14), BM1CS(21), BTH1CS(24) 36 DATA BY1 CS( 1) / .0320804710 0611908629E0 / 37 DATA BY1 CS( 2) / 1.2627078974 33500450E0 / 38 DATA BY1 CS( 3) / .0064999618 9992317500E0 / 39 DATA BY1 CS( 4) / -.0893616452 8860504117E0 / 40 DATA BY1 CS( 5) / .0132508812 2175709545E0 / 41 DATA BY1 CS( 6) / -.0008979059 1196483523E0 / 42 DATA BY1 CS( 7) / .0000364736 1487958306E0 / 43 DATA BY1 CS( 8) / -.0000010013 7438166600E0 / 44 DATA BY1 CS( 9) / .0000000199 4539657390E0 / 45 DATA BY1 CS(10) / -.0000000003 0230656018E0 / 46 DATA BY1 CS(11) / .0000000000 0360987815E0 / 47 DATA BY1 CS(12) / -.0000000000 0003487488E0 / 48 DATA BY1 CS(13) / .0000000000 0000027838E0 / 49 DATA BY1 CS(14) / -.0000000000 0000000186E0 / 50 DATA BM1 CS( 1) / .1047362510 931285E0 / 51 DATA BM1 CS( 2) / .0044244389 3702345E0 / 52 DATA BM1 CS( 3) / -.0000566163 9504035E0 / 53 DATA BM1 CS( 4) / .0000023134 9417339E0 / 54 DATA BM1 CS( 5) / -.0000001737 7182007E0 / 55 DATA BM1 CS( 6) / .0000000189 3209930E0 / 56 DATA BM1 CS( 7) / -.0000000026 5416023E0 / 57 DATA BM1 CS( 8) / .0000000004 4740209E0 / 58 DATA BM1 CS( 9) / -.0000000000 8691795E0 / 59 DATA BM1 CS(10) / .0000000000 1891492E0 / 60 DATA BM1 CS(11) / -.0000000000 0451884E0 / 61 DATA BM1 CS(12) / .0000000000 0116765E0 / 62 DATA BM1 CS(13) / -.0000000000 0032265E0 / 63 DATA BM1 CS(14) / .0000000000 0009450E0 / 64 DATA BM1 CS(15) / -.0000000000 0002913E0 / 65 DATA BM1 CS(16) / .0000000000 0000939E0 / 66 DATA BM1 CS(17) / -.0000000000 0000315E0 / 67 DATA BM1 CS(18) / .0000000000 0000109E0 / 68 DATA BM1 CS(19) / -.0000000000 0000039E0 / 69 DATA BM1 CS(20) / .0000000000 0000014E0 / 70 DATA BM1 CS(21) / -.0000000000 0000005E0 / 71 DATA BTH1CS( 1) / .7406014102 6313850E0 / 72 DATA BTH1CS( 2) / -.0045717556 59637690E0 / 73 DATA BTH1CS( 3) / .0001198185 10964326E0 / 74 DATA BTH1CS( 4) / -.0000069645 61891648E0 / 75 DATA BTH1CS( 5) / .0000006554 95621447E0 / 76 DATA BTH1CS( 6) / -.0000000840 66228945E0 / 77 DATA BTH1CS( 7) / .0000000133 76886564E0 / 78 DATA BTH1CS( 8) / -.0000000024 99565654E0 / 79 DATA BTH1CS( 9) / .0000000005 29495100E0 / 80 DATA BTH1CS(10) / -.0000000001 24135944E0 / 81 DATA BTH1CS(11) / .0000000000 31656485E0 / 82 DATA BTH1CS(12) / -.0000000000 08668640E0 / 83 DATA BTH1CS(13) / .0000000000 02523758E0 / 84 DATA BTH1CS(14) / -.0000000000 00775085E0 / 85 DATA BTH1CS(15) / .0000000000 00249527E0 / 86 DATA BTH1CS(16) / -.0000000000 00083773E0 / 87 DATA BTH1CS(17) / .0000000000 00029205E0 / 88 DATA BTH1CS(18) / -.0000000000 00010534E0 / 89 DATA BTH1CS(19) / .0000000000 00003919E0 / 90 DATA BTH1CS(20) / -.0000000000 00001500E0 / 91 DATA BTH1CS(21) / .0000000000 00000589E0 / 92 DATA BTH1CS(22) / -.0000000000 00000237E0 / 93 DATA BTH1CS(23) / .0000000000 00000097E0 / 94 DATA BTH1CS(24) / -.0000000000 00000040E0 / 95 DATA TWODPI / 0.6366197723 6758134E0 / 96 DATA PI4 / 0.7853981633 9744831E0 / 97 DATA NTY1, NTM1, NTTH1, XMIN, XSML, XMAX / 3*0, 3*0./ 98C***FIRST EXECUTABLE STATEMENT BESY1 99 IF (NTY1.NE.0) GO TO 10 100 NTY1 = INITS (BY1CS, 14, 0.1*R1MACH(3)) 101 NTM1 = INITS (BM1CS, 21, 0.1*R1MACH(3)) 102 NTTH1 = INITS (BTH1CS, 24, 0.1*R1MACH(3)) 103C 104 XMIN = 1.571*EXP ( AMAX1(ALOG(R1MACH(1)), -ALOG(R1MACH(2)))+.01) 105 XSML = SQRT (4.0*R1MACH(3)) 106 XMAX = 1.0/R1MACH(4) 107C 108 10 IF (X.LE.0.) CALL XERROR ( 'BESY1 X IS ZERO OR NEGATIVE', 29, 109 1 1, 2) 110 IF (X.GT.4.0) GO TO 20 111C 112 IF (X.LT.XMIN) CALL XERROR ( 'BESY1 X SO SMALL Y1 OVERFLOWS', 113 1 31, 3, 2) 114 Y = 0. 115 IF (X.GT.XSML) Y = X*X 116 BESY1 = TWODPI*ALOG(0.5*X)*BESJ1(X) + 117 1 (0.5 + CSEVL (.125*Y-1., BY1CS, NTY1))/X 118 RETURN 119C 120 20 IF (X.GT.XMAX) CALL XERROR ( 'BESY1 NO PRECISION BECAUSE X IS BI 121 1G', 37, 2, 2) 122C 123 Z = 32.0/X**2 - 1.0 124 AMPL = (0.75 + CSEVL (Z, BM1CS, NTM1)) / SQRT(X) 125 THETA = X - 3.0*PI4 + CSEVL (Z, BTH1CS, NTTH1) / X 126 BESY1 = AMPL * SIN (THETA) 127C 128 RETURN 129 END 130