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