source: sasview/src/sas/models/c_extension/cephes/struve.c @ 79492222

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 79492222 was 79492222, checked in by krzywon, 9 years ago

Changed the file and folder names to remove all SANS references.

  • Property mode set to 100644
File size: 4.4 KB
Line 
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/*
34Cephes Math Library Release 2.81:  June, 2000
35Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
36*/
37#include "mconf.h"
38#define DEBUG 0
39#ifdef ANSIPROT
40extern double gamma ( double );
41extern double pow ( double, double );
42extern double sqrt ( double );
43extern double yn ( int, double );
44extern double jv ( double, double );
45extern double fabs ( double );
46extern double floor ( double );
47extern double sin ( double );
48extern double cos ( double );
49double yv ( double, double );
50double onef2 (double, double, double, double, double * );
51double threef0 (double, double, double, double, double * );
52#else
53double gamma(), pow(), sqrt(), yn(), yv(), jv(), fabs(), floor();
54double sin(), cos();
55double onef2(), threef0();
56#endif
57static double stop = 1.37e-17;
58extern double MACHEP;
59
60double onef2( a, b, c, x, err )
61double a, b, c, x;
62double *err;
63{
64double n, a0, sum, t;
65double an, bn, cn, max, z;
66
67an = a;
68bn = b;
69cn = c;
70a0 = 1.0;
71sum = 1.0;
72n = 1.0;
73t = 1.0;
74max = 0.0;
75
76do
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 = fabs( a0 );
93        if( z > max )
94                max = z;
95        if( sum != 0 )
96                t = fabs( a0 / sum );
97        else
98                t = z;
99        }
100while( t > stop );
101
102done:
103
104*err = fabs( MACHEP*max /sum );
105
106#if DEBUG
107        printf(" onef2 cancellation error %.5E\n", *err );
108#endif
109
110goto xit;
111
112error:
113#if DEBUG
114printf("onef2 does not converge\n");
115#endif
116*err = 1.0e38;
117
118xit:
119
120#if DEBUG
121printf("onef2( %.2E %.2E %.2E %.5E ) =  %.3E  %.6E\n", a, b, c, x, n, sum);
122#endif
123return(sum);
124}
125
126
127
128
129double threef0( a, b, c, x, err )
130double a, b, c, x;
131double *err;
132{
133double n, a0, sum, t, conv, conv1;
134double an, bn, cn, max, z;
135
136an = a;
137bn = b;
138cn = c;
139a0 = 1.0;
140sum = 1.0;
141n = 1.0;
142t = 1.0;
143max = 0.0;
144conv = 1.0e38;
145conv1 = conv;
146
147do
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 = 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 = fabs( a0 / sum );
175        else
176                t = z;
177        }
178while( t > stop );
179
180done:
181
182t = fabs( MACHEP*max/sum );
183#if DEBUG
184        printf(" threef0 cancellation error %.5E\n", t );
185#endif
186
187max = fabs( conv/sum );
188if( max > t )
189        t = max;
190#if DEBUG
191        printf(" threef0 convergence %.5E\n", max );
192#endif
193
194goto xit;
195
196error:
197#if DEBUG
198printf("threef0 does not converge\n");
199#endif
200t = 1.0e38;
201
202xit:
203
204#if DEBUG
205printf("threef0( %.2E %.2E %.2E %.5E ) =  %.3E  %.6E\n", a, b, c, x, n, sum);
206#endif
207
208*err = t;
209return(sum);
210}
211
212
213
214
215extern double PI;
216
217double struve( v, x )
218double v, x;
219{
220double y, ya, f, g, h, t;
221double onef2err, threef0err;
222
223f = floor(v);
224if( (v < 0) && ( v-f == 0.5 ) )
225        {
226        y = jv( -v, x );
227        f = 1.0 - f;
228        g =  2.0 * floor(f/2.0);
229        if( g != f )
230                y = -y;
231        return(y);
232        }
233t = 0.25*x*x;
234f = fabs(x);
235g = 1.5 * fabs(v);
236if( (f > 30.0) && (f > g) )
237        {
238        onef2err = 1.0e38;
239        y = 0.0;
240        }
241else
242        {
243        y = onef2( 1.0, 1.5, 1.5+v, -t, &onef2err );
244        }
245
246if( (f < 18.0) || (x < 0.0) )
247        {
248        threef0err = 1.0e38;
249        ya = 0.0;
250        }
251else
252        {
253        ya = threef0( 1.0, 0.5, 0.5-v, -1.0/t, &threef0err );
254        }
255
256f = sqrt( PI );
257h = pow( 0.5*x, v-1.0 );
258
259if( onef2err <= threef0err )
260        {
261        g = gamma( v + 1.5 );
262        y = y * h * t / ( 0.5 * f * g );
263        return(y);
264        }
265else
266        {
267        g = 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
280double yv( v, x )
281double v, x;
282{
283double y, t;
284int n;
285
286y = floor( v );
287if( y == v )
288        {
289        n = v;
290        y = yn( n, x );
291        return( y );
292        }
293t = PI * v;
294y = (cos(t) * jv( v, x ) - jv( -v, x ))/sin(t);
295return( 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 */
Note: See TracBrowser for help on using the repository browser.