1 /* struve.c
2 *
3 * Struve function
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double v, x, y, struve();
10 *
11 * y = struve( v, x );
12 *
13 *
14 *
15 * DESCRIPTION:
16 *
17 * Computes the Struve function Hv(x) of order v, argument x.
18 * Negative x is rejected unless v is an integer.
19 *
20 * This module also contains the hypergeometric functions 1F2
21 * and 3F0 and a routine for the Bessel function Yv(x) with
22 * noninteger v.
23 *
24 *
25 *
26 * ACCURACY:
27 *
28 * Not accurately characterized, but spot checked against tables.
29 *
30 */
31
32
33 /*
34 Cephes Math Library Release 2.81: June, 2000
35 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
36 */
37 #include "mconf.h"
38 #define DEBUG 0
39 #ifdef ANSIPROT
40 extern double md_gamma ( double );
41 extern double md_pow ( double, double );
42 extern double sqrt ( double );
43 extern double md_yn ( int, double );
44 extern double jv ( double, double );
45 extern double md_fabs ( double );
46 extern double md_floor ( double );
47 extern double md_sin ( double );
48 extern double md_cos ( double );
49 double yv ( double, double );
50 double onef2 (double, double, double, double, double * );
51 double threef0 (double, double, double, double, double * );
52 #else
53 double md_gamma(), md_pow(), sqrt(), md_yn(), yv(), jv(), md_fabs(), md_floor();
54 double md_sin(), md_cos();
55 double onef2(), threef0();
56 #endif
57 static double stop = 1.37e-17;
58 extern double MACHEP;
59
onef2(a,b,c,x,err)60 double onef2( a, b, c, x, err )
61 double a, b, c, x;
62 double *err;
63 {
64 double n, a0, sum, t;
65 double an, bn, cn, max, z;
66
67 an = a;
68 bn = b;
69 cn = c;
70 a0 = 1.0;
71 sum = 1.0;
72 n = 1.0;
73 t = 1.0;
74 max = 0.0;
75
76 do
77 {
78 if( an == 0 )
79 goto done;
80 if( bn == 0 )
81 goto error;
82 if( cn == 0 )
83 goto error;
84 if( (a0 > 1.0e34) || (n > 200) )
85 goto error;
86 a0 *= (an * x) / (bn * cn * n);
87 sum += a0;
88 an += 1.0;
89 bn += 1.0;
90 cn += 1.0;
91 n += 1.0;
92 z = md_fabs( a0 );
93 if( z > max )
94 max = z;
95 if( sum != 0 )
96 t = md_fabs( a0 / sum );
97 else
98 t = z;
99 }
100 while( t > stop );
101
102 done:
103
104 *err = md_fabs( MACHEP*max /sum );
105
106 #if DEBUG
107 printf(" onef2 cancellation error %.5E\n", *err );
108 #endif
109
110 goto xit;
111
112 error:
113 #if DEBUG
114 printf("onef2 does not converge\n");
115 #endif
116 *err = 1.0e38;
117
118 xit:
119
120 #if DEBUG
121 printf("onef2( %.2E %.2E %.2E %.5E ) = %.3E %.6E\n", a, b, c, x, n, sum);
122 #endif
123 return(sum);
124 }
125
126
127
128
threef0(a,b,c,x,err)129 double threef0( a, b, c, x, err )
130 double a, b, c, x;
131 double *err;
132 {
133 double n, a0, sum, t, conv, conv1;
134 double an, bn, cn, max, z;
135
136 an = a;
137 bn = b;
138 cn = c;
139 a0 = 1.0;
140 sum = 1.0;
141 n = 1.0;
142 t = 1.0;
143 max = 0.0;
144 conv = 1.0e38;
145 conv1 = conv;
146
147 do
148 {
149 if( an == 0.0 )
150 goto done;
151 if( bn == 0.0 )
152 goto done;
153 if( cn == 0.0 )
154 goto done;
155 if( (a0 > 1.0e34) || (n > 200) )
156 goto error;
157 a0 *= (an * bn * cn * x) / n;
158 an += 1.0;
159 bn += 1.0;
160 cn += 1.0;
161 n += 1.0;
162 z = md_fabs( a0 );
163 if( z > max )
164 max = z;
165 if( z >= conv )
166 {
167 if( (z < max) && (z > conv1) )
168 goto done;
169 }
170 conv1 = conv;
171 conv = z;
172 sum += a0;
173 if( sum != 0 )
174 t = md_fabs( a0 / sum );
175 else
176 t = z;
177 }
178 while( t > stop );
179
180 done:
181
182 t = md_fabs( MACHEP*max/sum );
183 #if DEBUG
184 printf(" threef0 cancellation error %.5E\n", t );
185 #endif
186
187 max = md_fabs( conv/sum );
188 if( max > t )
189 t = max;
190 #if DEBUG
191 printf(" threef0 convergence %.5E\n", max );
192 #endif
193
194 goto xit;
195
196 error:
197 #if DEBUG
198 printf("threef0 does not converge\n");
199 #endif
200 t = 1.0e38;
201
202 xit:
203
204 #if DEBUG
205 printf("threef0( %.2E %.2E %.2E %.5E ) = %.3E %.6E\n", a, b, c, x, n, sum);
206 #endif
207
208 *err = t;
209 return(sum);
210 }
211
212
213
214
215 extern double PI;
216
struve(v,x)217 double struve( v, x )
218 double v, x;
219 {
220 double y, ya, f, g, h, t;
221 double onef2err, threef0err;
222
223 f = md_floor(v);
224 if( (v < 0) && ( v-f == 0.5 ) )
225 {
226 y = jv( -v, x );
227 f = 1.0 - f;
228 g = 2.0 * md_floor(f/2.0);
229 if( g != f )
230 y = -y;
231 return(y);
232 }
233 t = 0.25*x*x;
234 f = md_fabs(x);
235 g = 1.5 * md_fabs(v);
236 if( (f > 30.0) && (f > g) )
237 {
238 onef2err = 1.0e38;
239 y = 0.0;
240 }
241 else
242 {
243 y = onef2( 1.0, 1.5, 1.5+v, -t, &onef2err );
244 }
245
246 if( (f < 18.0) || (x < 0.0) )
247 {
248 threef0err = 1.0e38;
249 ya = 0.0;
250 }
251 else
252 {
253 ya = threef0( 1.0, 0.5, 0.5-v, -1.0/t, &threef0err );
254 }
255
256 f = sqrt( PI );
257 h = md_pow( 0.5*x, v-1.0 );
258
259 if( onef2err <= threef0err )
260 {
261 g = md_gamma( v + 1.5 );
262 y = y * h * t / ( 0.5 * f * g );
263 return(y);
264 }
265 else
266 {
267 g = md_gamma( v + 0.5 );
268 ya = ya * h / ( f * g );
269 ya = ya + yv( v, x );
270 return(ya);
271 }
272 }
273
274
275
276
277 /* Bessel function of noninteger order
278 */
279
yv(v,x)280 double yv( v, x )
281 double v, x;
282 {
283 double y, t;
284 int n;
285
286 y = md_floor( v );
287 if( y == v )
288 {
289 n = v;
290 y = md_yn( n, x );
291 return( y );
292 }
293 t = PI * v;
294 y = (md_cos(t) * jv( v, x ) - jv( -v, x ))/md_sin(t);
295 return( y );
296 }
297
298 /* Crossover points between ascending series and asymptotic series
299 * for Struve function
300 *
301 * v x
302 *
303 * 0 19.2
304 * 1 18.95
305 * 2 19.15
306 * 3 19.3
307 * 5 19.7
308 * 10 21.35
309 * 20 26.35
310 * 30 32.31
311 * 40 40.0
312 */
313