1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
|
/* jnf.c
*
* Bessel function of integer order
*
*
*
* SYNOPSIS:
*
* int n;
* float x, y, jnf();
*
* y = jnf( n, x );
*
*
*
* DESCRIPTION:
*
* Returns Bessel function of order n, where n is a
* (possibly negative) integer.
*
* The ratio of jn(x) to j0(x) is computed by backward
* recurrence. First the ratio jn/jn-1 is found by a
* continued fraction expansion. Then the recurrence
* relating successive orders is applied until j0 or j1 is
* reached.
*
* If n = 0 or 1 the routine for j0 or j1 is called
* directly.
*
*
*
* ACCURACY:
*
* Absolute error:
* arithmetic range # trials peak rms
* IEEE 0, 15 30000 3.6e-7 3.6e-8
*
*
* Not suitable for large n or x. Use jvf() instead.
*
*/
/* jn.c
Cephes Math Library Release 2.2: June, 1992
Copyright 1984, 1987, 1992 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
#include <math.h>
extern float MACHEPF;
float j0f(float), j1f(float);
float jnf( int n, float xx )
{
float x, pkm2, pkm1, pk, xk, r, ans, xinv, sign;
int k;
x = xx;
sign = 1.0;
if( n < 0 )
{
n = -n;
if( (n & 1) != 0 ) /* -1**n */
sign = -1.0;
}
if( n == 0 )
return( sign * j0f(x) );
if( n == 1 )
return( sign * j1f(x) );
if( n == 2 )
return( sign * (2.0 * j1f(x) / x - j0f(x)) );
/*
if( x < MACHEPF )
return( 0.0 );
*/
/* continued fraction */
k = 24;
pk = 2 * (n + k);
ans = pk;
xk = x * x;
do
{
pk -= 2.0;
ans = pk - (xk/ans);
}
while( --k > 0 );
/*ans = x/ans;*/
/* backward recurrence */
pk = 1.0;
/*pkm1 = 1.0/ans;*/
xinv = 1.0/x;
pkm1 = ans * xinv;
k = n-1;
r = (float )(2 * k);
do
{
pkm2 = (pkm1 * r - pk * x) * xinv;
pk = pkm1;
pkm1 = pkm2;
r -= 2.0;
}
while( --k > 0 );
r = pk;
if( r < 0 )
r = -r;
ans = pkm1;
if( ans < 0 )
ans = -ans;
if( r > ans ) /* if( fabs(pk) > fabs(pkm1) ) */
ans = sign * j1f(x)/pk;
else
ans = sign * j0f(x)/pkm1;
return( ans );
}
|
