1 /* j1.c
2 *
3 * Bessel function of order one
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double x, y, j1();
10 *
11 * y = j1( x );
12 *
13 *
14 *
15 * DESCRIPTION:
16 *
17 * Returns Bessel function of order one of the argument.
18 *
19 * The domain is divided into the intervals [0, 8] and
20 * (8, infinity). In the first interval a 24 term Chebyshev
21 * expansion is used. In the second, the asymptotic
22 * trigonometric representation is employed using two
23 * rational functions of degree 5/5.
24 *
25 *
26 *
27 * ACCURACY:
28 *
29 * Absolute error:
30 * arithmetic domain # trials peak rms
31 * DEC 0, 30 10000 4.0e-17 1.1e-17
32 * IEEE 0, 30 30000 2.6e-16 1.1e-16
33 *
34 *
35 */
36 /* y1.c
37 *
38 * Bessel function of second kind of order one
39 *
40 *
41 *
42 * SYNOPSIS:
43 *
44 * double x, y, y1();
45 *
46 * y = y1( x );
47 *
48 *
49 *
50 * DESCRIPTION:
51 *
52 * Returns Bessel function of the second kind of order one
53 * of the argument.
54 *
55 * The domain is divided into the intervals [0, 8] and
56 * (8, infinity). In the first interval a 25 term Chebyshev
57 * expansion is used, and a call to j1() is required.
58 * In the second, the asymptotic trigonometric representation
59 * is employed using two rational functions of degree 5/5.
60 *
61 *
62 *
63 * ACCURACY:
64 *
65 * Absolute error:
66 * arithmetic domain # trials peak rms
67 * DEC 0, 30 10000 8.6e-17 1.3e-17
68 * IEEE 0, 30 30000 1.0e-15 1.3e-16
69 *
70 * (error criterion relative when |y1| > 1).
71 *
72 */
73
74
75 /*
76 Cephes Math Library Release 2.1: January, 1989
77 Copyright 1984, 1987, 1989 by Stephen L. Moshier
78 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
79 */
80
81 /*
82 #define PIO4 .78539816339744830962
83 #define THPIO4 2.35619449019234492885
84 #define SQ2OPI .79788456080286535588
85 */
86
87 #include "mconf.h"
88
89 static double RP[4] = {
90 -8.99971225705559398224E8,
91 4.52228297998194034323E11,
92 -7.27494245221818276015E13,
93 3.68295732863852883286E15,
94 };
95 static double RQ[8] = {
96 /* 1.00000000000000000000E0,*/
97 6.20836478118054335476E2,
98 2.56987256757748830383E5,
99 8.35146791431949253037E7,
100 2.21511595479792499675E10,
101 4.74914122079991414898E12,
102 7.84369607876235854894E14,
103 8.95222336184627338078E16,
104 5.32278620332680085395E18,
105 };
106
107 static double PP[7] = {
108 7.62125616208173112003E-4,
109 7.31397056940917570436E-2,
110 1.12719608129684925192E0,
111 5.11207951146807644818E0,
112 8.42404590141772420927E0,
113 5.21451598682361504063E0,
114 1.00000000000000000254E0,
115 };
116 static double PQ[7] = {
117 5.71323128072548699714E-4,
118 6.88455908754495404082E-2,
119 1.10514232634061696926E0,
120 5.07386386128601488557E0,
121 8.39985554327604159757E0,
122 5.20982848682361821619E0,
123 9.99999999999999997461E-1,
124 };
125
126 static double QP[8] = {
127 5.10862594750176621635E-2,
128 4.98213872951233449420E0,
129 7.58238284132545283818E1,
130 3.66779609360150777800E2,
131 7.10856304998926107277E2,
132 5.97489612400613639965E2,
133 2.11688757100572135698E2,
134 2.52070205858023719784E1,
135 };
136 static double QQ[7] = {
137 /* 1.00000000000000000000E0,*/
138 7.42373277035675149943E1,
139 1.05644886038262816351E3,
140 4.98641058337653607651E3,
141 9.56231892404756170795E3,
142 7.99704160447350683650E3,
143 2.82619278517639096600E3,
144 3.36093607810698293419E2,
145 };
146
147
148 static double YP[6] = {
149 1.26320474790178026440E9,
150 -6.47355876379160291031E11,
151 1.14509511541823727583E14,
152 -8.12770255501325109621E15,
153 2.02439475713594898196E17,
154 -7.78877196265950026825E17,
155 };
156 static double YQ[8] = {
157 /* 1.00000000000000000000E0,*/
158 5.94301592346128195359E2,
159 2.35564092943068577943E5,
160 7.34811944459721705660E7,
161 1.87601316108706159478E10,
162 3.88231277496238566008E12,
163 6.20557727146953693363E14,
164 6.87141087355300489866E16,
165 3.97270608116560655612E18,
166 };
167
168 static double Z1 = 1.46819706421238932572E1;
169 static double Z2 = 4.92184563216946036703E1;
170
j1(x)171 double j1(x)
172 double x;
173 {
174 extern double THPIO4, SQ2OPI;
175 double polevl(), p1evl();
176 double w, z, p, q, xn;
177 double sin(), cos(), sqrt();
178
179 w = x;
180 if( x < 0 )
181 w = -x;
182
183 if( w <= 5.0 )
184 {
185 z = x * x;
186 w = polevl( z, RP, 3 ) / p1evl( z, RQ, 8 );
187 w = w * x * (z - Z1) * (z - Z2);
188 return( w );
189 }
190
191 w = 5.0/x;
192 z = w * w;
193 p = polevl( z, PP, 6)/polevl( z, PQ, 6 );
194 q = polevl( z, QP, 7)/p1evl( z, QQ, 7 );
195 xn = x - THPIO4;
196 p = p * cos(xn) - w * q * sin(xn);
197 return( p * SQ2OPI / sqrt(x) );
198 }
199
200
201
202
203 extern double MAXNUM;
204
y1(x)205 double y1(x)
206 double x;
207 {
208 extern double TWOOPI, THPIO4, SQ2OPI;
209 double polevl(), p1evl();
210 double w, z, p, q, xn;
211 double j1(), log(), sin(), cos(), sqrt();
212
213
214 if( x <= 5.0 )
215 {
216 if( x <= 0.0 )
217 {
218 mtherr( "y1", DOMAIN );
219 return( -MAXNUM );
220 }
221 z = x * x;
222 w = x * (polevl( z, YP, 5 ) / p1evl( z, YQ, 8 ));
223 w += TWOOPI * ( j1(x) * log(x) - 1.0/x );
224 return( w );
225 }
226
227 w = 5.0/x;
228 z = w * w;
229 p = polevl( z, PP, 6)/polevl( z, PQ, 6 );
230 q = polevl( z, QP, 7)/p1evl( z, QQ, 7 );
231 xn = x - THPIO4;
232 p = p * sin(xn) + w * q * cos(xn);
233 return( p * SQ2OPI / sqrt(x) );
234 }
235