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 zero
7
8 j0(x) returns the value of J0(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 J0 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 #5849 (19.22D)
23 #6549 (19.25D)
24 #6949 (19.41D)
25
26 y0(x) returns the value of Y0(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, j0.
32
33 The values of Y0 have not been checked
34 to more than ten places.
35
36 Coefficients are from Hart & Cheney.
37 #6245 (18.78D)
38 #6549 (19.25D)
39 #6949 (19.41D)
40 */
41
42 static double pzero, qzero;
43 static double tpi = .6366197723675813430755350535e0;
44 static double pio4 = .7853981633974483096156608458e0;
45 static double p1[] = {
46 0.4933787251794133561816813446e21,
47 -.1179157629107610536038440800e21,
48 0.6382059341072356562289432465e19,
49 -.1367620353088171386865416609e18,
50 0.1434354939140344111664316553e16,
51 -.8085222034853793871199468171e13,
52 0.2507158285536881945555156435e11,
53 -.4050412371833132706360663322e8,
54 0.2685786856980014981415848441e5,
55 };
56 static double q1[] = {
57 0.4933787251794133562113278438e21,
58 0.5428918384092285160200195092e19,
59 0.3024635616709462698627330784e17,
60 0.1127756739679798507056031594e15,
61 0.3123043114941213172572469442e12,
62 0.6699987672982239671814028660e9,
63 0.1114636098462985378182402543e7,
64 0.1363063652328970604442810507e4,
65 1.0
66 };
67 static double p2[] = {
68 0.5393485083869438325262122897e7,
69 0.1233238476817638145232406055e8,
70 0.8413041456550439208464315611e7,
71 0.2016135283049983642487182349e7,
72 0.1539826532623911470917825993e6,
73 0.2485271928957404011288128951e4,
74 0.0,
75 };
76 static double q2[] = {
77 0.5393485083869438325560444960e7,
78 0.1233831022786324960844856182e8,
79 0.8426449050629797331554404810e7,
80 0.2025066801570134013891035236e7,
81 0.1560017276940030940592769933e6,
82 0.2615700736920839685159081813e4,
83 1.0,
84 };
85 static double p3[] = {
86 -.3984617357595222463506790588e4,
87 -.1038141698748464093880530341e5,
88 -.8239066313485606568803548860e4,
89 -.2365956170779108192723612816e4,
90 -.2262630641933704113967255053e3,
91 -.4887199395841261531199129300e1,
92 0.0,
93 };
94 static double q3[] = {
95 0.2550155108860942382983170882e6,
96 0.6667454239319826986004038103e6,
97 0.5332913634216897168722255057e6,
98 0.1560213206679291652539287109e6,
99 0.1570489191515395519392882766e5,
100 0.4087714673983499223402830260e3,
101 1.0,
102 };
103 static double p4[] = {
104 -.2750286678629109583701933175e20,
105 0.6587473275719554925999402049e20,
106 -.5247065581112764941297350814e19,
107 0.1375624316399344078571335453e18,
108 -.1648605817185729473122082537e16,
109 0.1025520859686394284509167421e14,
110 -.3436371222979040378171030138e11,
111 0.5915213465686889654273830069e8,
112 -.4137035497933148554125235152e5,
113 };
114 static double q4[] = {
115 0.3726458838986165881989980e21,
116 0.4192417043410839973904769661e19,
117 0.2392883043499781857439356652e17,
118 0.9162038034075185262489147968e14,
119 0.2613065755041081249568482092e12,
120 0.5795122640700729537480087915e9,
121 0.1001702641288906265666651753e7,
122 0.1282452772478993804176329391e4,
123 1.0,
124 };
125
126 static void
asympt(double arg)127 asympt(double arg)
128 {
129 double zsq, n, d;
130 int i;
131
132 zsq = 64 / (arg*arg);
133 for(n=0,d=0,i=6;i>=0;i--) {
134 n = n*zsq + p2[i];
135 d = d*zsq + q2[i];
136 }
137 pzero = n/d;
138 for(n=0,d=0,i=6;i>=0;i--) {
139 n = n*zsq + p3[i];
140 d = d*zsq + q3[i];
141 }
142 qzero = (8/arg)*(n/d);
143 }
144
145 double
j0(double arg)146 j0(double arg)
147 {
148 double argsq, n, d;
149 int i;
150
151 if(arg < 0)
152 arg = -arg;
153 if(arg > 8) {
154 asympt(arg);
155 n = arg - pio4;
156 return sqrt(tpi/arg)*(pzero*cos(n) - qzero*sin(n));
157 }
158 argsq = arg*arg;
159 for(n=0,d=0,i=8;i>=0;i--) {
160 n = n*argsq + p1[i];
161 d = d*argsq + q1[i];
162 }
163 return n/d;
164 }
165
166 double
y0(double arg)167 y0(double arg)
168 {
169 double argsq, n, d;
170 int i;
171
172 errno = 0;
173 if(arg <= 0) {
174 errno = EDOM;
175 return(-HUGE_VAL);
176 }
177 if(arg > 8) {
178 asympt(arg);
179 n = arg - pio4;
180 return sqrt(tpi/arg)*(pzero*sin(n) + qzero*cos(n));
181 }
182 argsq = arg*arg;
183 for(n=0,d=0,i=8;i>=0;i--) {
184 n = n*argsq + p4[i];
185 d = d*argsq + q4[i];
186 }
187 return n/d + tpi*j0(arg)*log(arg);
188 }
189