1 #include <math.h>
2 #include <errno.h>
3 /*
4 floating point Bessel's function
5 of the first and second kinds
6 of order one
7
8 j1(x) returns the value of J1(x)
9 for all real values of x.
10
11 There are no error returns.
12 Calls sin, cos, sqrt.
13
14 There is a niggling bug in J1 which
15 causes errors up to 2e-16 for x in the
16 interval [-8,8].
17 The bug is caused by an inappropriate order
18 of summation of the series. rhm will fix it
19 someday.
20
21 Coefficients are from Hart & Cheney.
22 #6050 (20.98D)
23 #6750 (19.19D)
24 #7150 (19.35D)
25
26 y1(x) returns the value of Y1(x)
27 for positive real values of x.
28 For x<=0, error number EDOM is set and a
29 large negative value is returned.
30
31 Calls sin, cos, sqrt, log, j1.
32
33 The values of Y1 have not been checked
34 to more than ten places.
35
36 Coefficients are from Hart & Cheney.
37 #6447 (22.18D)
38 #6750 (19.19D)
39 #7150 (19.35D)
40 */
41
42 static double pzero, qzero;
43 static double tpi = .6366197723675813430755350535e0;
44 static double pio4 = .7853981633974483096156608458e0;
45 static double p1[] = {
46 0.581199354001606143928050809e21,
47 -.6672106568924916298020941484e20,
48 0.2316433580634002297931815435e19,
49 -.3588817569910106050743641413e17,
50 0.2908795263834775409737601689e15,
51 -.1322983480332126453125473247e13,
52 0.3413234182301700539091292655e10,
53 -.4695753530642995859767162166e7,
54 0.2701122710892323414856790990e4,
55 };
56 static double q1[] = {
57 0.1162398708003212287858529400e22,
58 0.1185770712190320999837113348e20,
59 0.6092061398917521746105196863e17,
60 0.2081661221307607351240184229e15,
61 0.5243710262167649715406728642e12,
62 0.1013863514358673989967045588e10,
63 0.1501793594998585505921097578e7,
64 0.1606931573481487801970916749e4,
65 1.0,
66 };
67 static double p2[] = {
68 -.4435757816794127857114720794e7,
69 -.9942246505077641195658377899e7,
70 -.6603373248364939109255245434e7,
71 -.1523529351181137383255105722e7,
72 -.1098240554345934672737413139e6,
73 -.1611616644324610116477412898e4,
74 0.0,
75 };
76 static double q2[] = {
77 -.4435757816794127856828016962e7,
78 -.9934124389934585658967556309e7,
79 -.6585339479723087072826915069e7,
80 -.1511809506634160881644546358e7,
81 -.1072638599110382011903063867e6,
82 -.1455009440190496182453565068e4,
83 1.0,
84 };
85 static double p3[] = {
86 0.3322091340985722351859704442e5,
87 0.8514516067533570196555001171e5,
88 0.6617883658127083517939992166e5,
89 0.1849426287322386679652009819e5,
90 0.1706375429020768002061283546e4,
91 0.3526513384663603218592175580e2,
92 0.0,
93 };
94 static double q3[] = {
95 0.7087128194102874357377502472e6,
96 0.1819458042243997298924553839e7,
97 0.1419460669603720892855755253e7,
98 0.4002944358226697511708610813e6,
99 0.3789022974577220264142952256e5,
100 0.8638367769604990967475517183e3,
101 1.0,
102 };
103 static double p4[] = {
104 -.9963753424306922225996744354e23,
105 0.2655473831434854326894248968e23,
106 -.1212297555414509577913561535e22,
107 0.2193107339917797592111427556e20,
108 -.1965887462722140658820322248e18,
109 0.9569930239921683481121552788e15,
110 -.2580681702194450950541426399e13,
111 0.3639488548124002058278999428e10,
112 -.2108847540133123652824139923e7,
113 0.0,
114 };
115 static double q4[] = {
116 0.5082067366941243245314424152e24,
117 0.5435310377188854170800653097e22,
118 0.2954987935897148674290758119e20,
119 0.1082258259408819552553850180e18,
120 0.2976632125647276729292742282e15,
121 0.6465340881265275571961681500e12,
122 0.1128686837169442121732366891e10,
123 0.1563282754899580604737366452e7,
124 0.1612361029677000859332072312e4,
125 1.0,
126 };
127
128 static void
asympt(double arg)129 asympt(double arg)
130 {
131 double zsq, n, d;
132 int i;
133
134 zsq = 64/(arg*arg);
135 for(n=0,d=0,i=6;i>=0;i--) {
136 n = n*zsq + p2[i];
137 d = d*zsq + q2[i];
138 }
139 pzero = n/d;
140 for(n=0,d=0,i=6;i>=0;i--) {
141 n = n*zsq + p3[i];
142 d = d*zsq + q3[i];
143 }
144 qzero = (8/arg)*(n/d);
145 }
146
147 double
j1(double arg)148 j1(double arg)
149 {
150 double xsq, n, d, x;
151 int i;
152
153 x = arg;
154 if(x < 0)
155 x = -x;
156 if(x > 8) {
157 asympt(x);
158 n = x - 3*pio4;
159 n = sqrt(tpi/x)*(pzero*cos(n) - qzero*sin(n));
160 if(arg < 0)
161 n = -n;
162 return n;
163 }
164 xsq = x*x;
165 for(n=0,d=0,i=8;i>=0;i--) {
166 n = n*xsq + p1[i];
167 d = d*xsq + q1[i];
168 }
169 return arg*n/d;
170 }
171
172 double
y1(double arg)173 y1(double arg)
174 {
175 double xsq, n, d, x;
176 int i;
177
178 errno = 0;
179 x = arg;
180 if(x <= 0) {
181 errno = EDOM;
182 return -HUGE_VAL;
183 }
184 if(x > 8) {
185 asympt(x);
186 n = x - 3*pio4;
187 return sqrt(tpi/x)*(pzero*sin(n) + qzero*cos(n));
188 }
189 xsq = x*x;
190 for(n=0,d=0,i=9;i>=0;i--) {
191 n = n*xsq + p4[i];
192 d = d*xsq + q4[i];
193 }
194 return x*n/d + tpi*(j1(x)*log(x)-1/x);
195 }
196