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