← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/python-scipy/oneiric-proposed

« back to all changes in this revision

Viewing changes to Lib/special/cephes/struve.c

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/special/cephes/struve.c
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 gamma ( double );
41
 
extern double pow ( double, double );
42
 
extern double sqrt ( double );
43
 
extern double yn ( int, double );
44
 
extern double jv ( double, double );
45
 
extern double fabs ( double );
46
 
extern double floor ( double );
47
 
extern double sin ( double );
48
 
extern double 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 gamma(), pow(), sqrt(), yn(), yv(), jv(), fabs(), floor();
54
 
double sin(), cos();
55
 
double onef2(), threef0();
56
 
#endif
57
 
static double stop = 1.37e-17;
58
 
extern double MACHEP, INFINITY;
59
 
 
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 = fabs( a0 );
93
 
        if( z > max )
94
 
                max = z;
95
 
        if( sum != 0 )
96
 
                t = fabs( a0 / sum );
97
 
        else
98
 
                t = z;
99
 
        }
100
 
while( t > stop );
101
 
 
102
 
done:
103
 
 
104
 
*err = 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
 
 
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 = 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
 
        }
178
 
while( t > stop );
179
 
 
180
 
done:
181
 
 
182
 
t = fabs( MACHEP*max/sum );
183
 
#if DEBUG
184
 
        printf(" threef0 cancellation error %.5E\n", t );
185
 
#endif
186
 
 
187
 
max = 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
 
 
217
 
double struve( v, x )
218
 
double v, x;
219
 
{
220
 
double y, ya, f, g, h, t;
221
 
double onef2err, threef0err;
222
 
 
223
 
if (x == 0.0) {
224
 
    if (v > -1) {
225
 
        return 0.0;
226
 
    } else if (v < -1) {
227
 
        if ((int)(floor(0.5-v)-1) % 2) return -INFINITY;
228
 
        else return INFINITY;
229
 
    } else {
230
 
        return 2.0/PI;
231
 
    }
232
 
}
233
 
 
234
 
f = floor(v);
235
 
if( (v < 0) && ( v-f == 0.5 ) )
236
 
        {
237
 
        y = jv( -v, x );
238
 
        f = 1.0 - f;
239
 
        g =  2.0 * floor(f/2.0);
240
 
        if( g != f )
241
 
                y = -y;
242
 
        return(y);
243
 
        }
244
 
t = 0.25*x*x;
245
 
f = fabs(x);
246
 
g = 1.5 * fabs(v);
247
 
if( (f > 30.0) && (f > g) )
248
 
        {
249
 
        onef2err = 1.0e38;
250
 
        y = 0.0;
251
 
        }
252
 
else
253
 
        {
254
 
        y = onef2( 1.0, 1.5, 1.5+v, -t, &onef2err );
255
 
        }
256
 
 
257
 
if( (f < 18.0) || (x < 0.0) )
258
 
        {
259
 
        threef0err = 1.0e38;
260
 
        ya = 0.0;
261
 
        }
262
 
else
263
 
        {
264
 
        ya = threef0( 1.0, 0.5, 0.5-v, -1.0/t, &threef0err );
265
 
        }
266
 
 
267
 
f = sqrt( PI );
268
 
h = pow( 0.5*x, v-1.0 );
269
 
 
270
 
if( onef2err <= threef0err )
271
 
        {
272
 
        g = gamma( v + 1.5 );
273
 
        y = y * h * t / ( 0.5 * f * g );
274
 
        return(y);
275
 
        }
276
 
else
277
 
        {
278
 
        g = gamma( v + 0.5 );
279
 
        ya = ya * h / ( f * g );
280
 
        ya = ya + yv( v, x );
281
 
        return(ya);
282
 
        }
283
 
}
284
 
 
285
 
 
286
 
 
287
 
 
288
 
/* Bessel function of noninteger order
289
 
 */
290
 
 
291
 
double yv( v, x )
292
 
double v, x;
293
 
{
294
 
double y, t;
295
 
int n;
296
 
 
297
 
y = floor( v );
298
 
if( y == v )
299
 
        {
300
 
        n = v;
301
 
        y = yn( n, x );
302
 
        return( y );
303
 
        }
304
 
t = PI * v;
305
 
y = (cos(t) * jv( v, x ) - jv( -v, x ))/sin(t);
306
 
return( y );
307
 
}
308
 
 
309
 
/* Crossover points between ascending series and asymptotic series
310
 
 * for Struve function
311
 
 *
312
 
 *       v       x
313
 
 * 
314
 
 *       0      19.2
315
 
 *       1      18.95
316
 
 *       2      19.15
317
 
 *       3      19.3
318
 
 *       5      19.7
319
 
 *      10      21.35
320
 
 *      20      26.35
321
 
 *      30      32.31
322
 
 *      40      40.0
323
 
 */