1 | /* j0.c |
---|
2 | * |
---|
3 | * Bessel function of order zero |
---|
4 | * |
---|
5 | * |
---|
6 | * |
---|
7 | * SYNOPSIS: |
---|
8 | * |
---|
9 | * double x, y, j0(); |
---|
10 | * |
---|
11 | * y = j0( x ); |
---|
12 | * |
---|
13 | * |
---|
14 | * |
---|
15 | * DESCRIPTION: |
---|
16 | * |
---|
17 | * Returns Bessel function of order zero of the argument. |
---|
18 | * |
---|
19 | * The domain is divided into the intervals [0, 5] and |
---|
20 | * (5, infinity). In the first interval the following rational |
---|
21 | * approximation is used: |
---|
22 | * |
---|
23 | * |
---|
24 | * 2 2 |
---|
25 | * (w - r ) (w - r ) P (w) / Q (w) |
---|
26 | * 1 2 3 8 |
---|
27 | * |
---|
28 | * 2 |
---|
29 | * where w = x and the two r's are zeros of the function. |
---|
30 | * |
---|
31 | * In the second interval, the Hankel asymptotic expansion |
---|
32 | * is employed with two rational functions of degree 6/6 |
---|
33 | * and 7/7. |
---|
34 | * |
---|
35 | * |
---|
36 | * |
---|
37 | * ACCURACY: |
---|
38 | * |
---|
39 | * Absolute error: |
---|
40 | * arithmetic domain # trials peak rms |
---|
41 | * DEC 0, 30 10000 4.4e-17 6.3e-18 |
---|
42 | * IEEE 0, 30 60000 4.2e-16 1.1e-16 |
---|
43 | * |
---|
44 | */ |
---|
45 | |
---|
46 | /* |
---|
47 | Cephes Math Library Release 2.8: June, 2000 |
---|
48 | Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier |
---|
49 | */ |
---|
50 | |
---|
51 | /* Note: all coefficients satisfy the relative error criterion |
---|
52 | * except YP, YQ which are designed for absolute error. */ |
---|
53 | |
---|
54 | double J0(double x) { |
---|
55 | |
---|
56 | //Cephes single precission |
---|
57 | #if FLOAT_SIZE>4 |
---|
58 | double w, z, p, q, xn; |
---|
59 | |
---|
60 | //const double TWOOPI = 6.36619772367581343075535E-1; |
---|
61 | const double SQ2OPI = 7.9788456080286535587989E-1; |
---|
62 | const double PIO4 = 7.85398163397448309616E-1; |
---|
63 | |
---|
64 | const double DR1 = 5.78318596294678452118E0; |
---|
65 | const double DR2 = 3.04712623436620863991E1; |
---|
66 | |
---|
67 | |
---|
68 | if( x < 0 ) |
---|
69 | x = -x; |
---|
70 | |
---|
71 | if( x <= 5.0 ) { |
---|
72 | z = x * x; |
---|
73 | if( x < 1.0e-5 ) |
---|
74 | return( 1.0 - z/4.0 ); |
---|
75 | |
---|
76 | p = (z - DR1) * (z - DR2); |
---|
77 | p = p * polevl( z, RPJ0, 3)/p1evl( z, RQJ0, 8 ); |
---|
78 | return( p ); |
---|
79 | } |
---|
80 | |
---|
81 | w = 5.0/x; |
---|
82 | q = 25.0/(x*x); |
---|
83 | p = polevl( q, PPJ0, 6)/polevl( q, PQJ0, 6 ); |
---|
84 | q = polevl( q, QPJ0, 7)/p1evl( q, QQJ0, 7 ); |
---|
85 | xn = x - PIO4; |
---|
86 | |
---|
87 | double sn, cn; |
---|
88 | SINCOS(xn, sn, cn); |
---|
89 | p = p * cn - w * q * sn; |
---|
90 | |
---|
91 | return( p * SQ2OPI / sqrt(x) ); |
---|
92 | //Cephes single precission |
---|
93 | #else |
---|
94 | double xx, w, z, p, q, xn; |
---|
95 | |
---|
96 | //const double YZ1 = 0.43221455686510834878; |
---|
97 | //const double YZ2 = 22.401876406482861405; |
---|
98 | //const double YZ3 = 64.130620282338755553; |
---|
99 | const double DR1 = 5.78318596294678452118; |
---|
100 | const double PIO4F = 0.7853981633974483096; |
---|
101 | |
---|
102 | if( x < 0 ) |
---|
103 | xx = -x; |
---|
104 | else |
---|
105 | xx = x; |
---|
106 | |
---|
107 | if( x <= 2.0 ) { |
---|
108 | z = xx * xx; |
---|
109 | if( x < 1.0e-3 ) |
---|
110 | return( 1.0 - 0.25*z ); |
---|
111 | |
---|
112 | p = (z-DR1) * polevl( z, JPJ0, 4); |
---|
113 | return( p ); |
---|
114 | } |
---|
115 | |
---|
116 | q = 1.0/x; |
---|
117 | w = sqrt(q); |
---|
118 | |
---|
119 | p = w * polevl( q, MOJ0, 7); |
---|
120 | w = q*q; |
---|
121 | xn = q * polevl( w, PHJ0, 7) - PIO4F; |
---|
122 | p = p * cos(xn + xx); |
---|
123 | return(p); |
---|
124 | #endif |
---|
125 | |
---|
126 | } |
---|
127 | |
---|