1*9933Ssam /* @(#)j0.c 4.1 12/25/82 */
2*9933Ssam
3*9933Ssam /*
4*9933Ssam floating point Bessel's function
5*9933Ssam of the first and second kinds
6*9933Ssam of order zero
7*9933Ssam
8*9933Ssam j0(x) returns the value of J0(x)
9*9933Ssam for all real values of x.
10*9933Ssam
11*9933Ssam There are no error returns.
12*9933Ssam Calls sin, cos, sqrt.
13*9933Ssam
14*9933Ssam There is a niggling bug in J0 which
15*9933Ssam causes errors up to 2e-16 for x in the
16*9933Ssam interval [-8,8].
17*9933Ssam The bug is caused by an inappropriate order
18*9933Ssam of summation of the series. rhm will fix it
19*9933Ssam someday.
20*9933Ssam
21*9933Ssam Coefficients are from Hart & Cheney.
22*9933Ssam #5849 (19.22D)
23*9933Ssam #6549 (19.25D)
24*9933Ssam #6949 (19.41D)
25*9933Ssam
26*9933Ssam y0(x) returns the value of Y0(x)
27*9933Ssam for positive real values of x.
28*9933Ssam For x<=0, error number EDOM is set and a
29*9933Ssam large negative value is returned.
30*9933Ssam
31*9933Ssam Calls sin, cos, sqrt, log, j0.
32*9933Ssam
33*9933Ssam The values of Y0 have not been checked
34*9933Ssam to more than ten places.
35*9933Ssam
36*9933Ssam Coefficients are from Hart & Cheney.
37*9933Ssam #6245 (18.78D)
38*9933Ssam #6549 (19.25D)
39*9933Ssam #6949 (19.41D)
40*9933Ssam */
41*9933Ssam
42*9933Ssam #include <math.h>
43*9933Ssam #include <errno.h>
44*9933Ssam
45*9933Ssam int errno;
46*9933Ssam static double pzero, qzero;
47*9933Ssam static double tpi = .6366197723675813430755350535e0;
48*9933Ssam static double pio4 = .7853981633974483096156608458e0;
49*9933Ssam static double p1[] = {
50*9933Ssam 0.4933787251794133561816813446e21,
51*9933Ssam -.1179157629107610536038440800e21,
52*9933Ssam 0.6382059341072356562289432465e19,
53*9933Ssam -.1367620353088171386865416609e18,
54*9933Ssam 0.1434354939140344111664316553e16,
55*9933Ssam -.8085222034853793871199468171e13,
56*9933Ssam 0.2507158285536881945555156435e11,
57*9933Ssam -.4050412371833132706360663322e8,
58*9933Ssam 0.2685786856980014981415848441e5,
59*9933Ssam };
60*9933Ssam static double q1[] = {
61*9933Ssam 0.4933787251794133562113278438e21,
62*9933Ssam 0.5428918384092285160200195092e19,
63*9933Ssam 0.3024635616709462698627330784e17,
64*9933Ssam 0.1127756739679798507056031594e15,
65*9933Ssam 0.3123043114941213172572469442e12,
66*9933Ssam 0.6699987672982239671814028660e9,
67*9933Ssam 0.1114636098462985378182402543e7,
68*9933Ssam 0.1363063652328970604442810507e4,
69*9933Ssam 1.0
70*9933Ssam };
71*9933Ssam static double p2[] = {
72*9933Ssam 0.5393485083869438325262122897e7,
73*9933Ssam 0.1233238476817638145232406055e8,
74*9933Ssam 0.8413041456550439208464315611e7,
75*9933Ssam 0.2016135283049983642487182349e7,
76*9933Ssam 0.1539826532623911470917825993e6,
77*9933Ssam 0.2485271928957404011288128951e4,
78*9933Ssam 0.0,
79*9933Ssam };
80*9933Ssam static double q2[] = {
81*9933Ssam 0.5393485083869438325560444960e7,
82*9933Ssam 0.1233831022786324960844856182e8,
83*9933Ssam 0.8426449050629797331554404810e7,
84*9933Ssam 0.2025066801570134013891035236e7,
85*9933Ssam 0.1560017276940030940592769933e6,
86*9933Ssam 0.2615700736920839685159081813e4,
87*9933Ssam 1.0,
88*9933Ssam };
89*9933Ssam static double p3[] = {
90*9933Ssam -.3984617357595222463506790588e4,
91*9933Ssam -.1038141698748464093880530341e5,
92*9933Ssam -.8239066313485606568803548860e4,
93*9933Ssam -.2365956170779108192723612816e4,
94*9933Ssam -.2262630641933704113967255053e3,
95*9933Ssam -.4887199395841261531199129300e1,
96*9933Ssam 0.0,
97*9933Ssam };
98*9933Ssam static double q3[] = {
99*9933Ssam 0.2550155108860942382983170882e6,
100*9933Ssam 0.6667454239319826986004038103e6,
101*9933Ssam 0.5332913634216897168722255057e6,
102*9933Ssam 0.1560213206679291652539287109e6,
103*9933Ssam 0.1570489191515395519392882766e5,
104*9933Ssam 0.4087714673983499223402830260e3,
105*9933Ssam 1.0,
106*9933Ssam };
107*9933Ssam static double p4[] = {
108*9933Ssam -.2750286678629109583701933175e20,
109*9933Ssam 0.6587473275719554925999402049e20,
110*9933Ssam -.5247065581112764941297350814e19,
111*9933Ssam 0.1375624316399344078571335453e18,
112*9933Ssam -.1648605817185729473122082537e16,
113*9933Ssam 0.1025520859686394284509167421e14,
114*9933Ssam -.3436371222979040378171030138e11,
115*9933Ssam 0.5915213465686889654273830069e8,
116*9933Ssam -.4137035497933148554125235152e5,
117*9933Ssam };
118*9933Ssam static double q4[] = {
119*9933Ssam 0.3726458838986165881989980e21,
120*9933Ssam 0.4192417043410839973904769661e19,
121*9933Ssam 0.2392883043499781857439356652e17,
122*9933Ssam 0.9162038034075185262489147968e14,
123*9933Ssam 0.2613065755041081249568482092e12,
124*9933Ssam 0.5795122640700729537480087915e9,
125*9933Ssam 0.1001702641288906265666651753e7,
126*9933Ssam 0.1282452772478993804176329391e4,
127*9933Ssam 1.0,
128*9933Ssam };
129*9933Ssam
130*9933Ssam double
j0(arg)131*9933Ssam j0(arg) double arg;{
132*9933Ssam double argsq, n, d;
133*9933Ssam double sin(), cos(), sqrt();
134*9933Ssam int i;
135*9933Ssam
136*9933Ssam if(arg < 0.) arg = -arg;
137*9933Ssam if(arg > 8.){
138*9933Ssam asympt(arg);
139*9933Ssam n = arg - pio4;
140*9933Ssam return(sqrt(tpi/arg)*(pzero*cos(n) - qzero*sin(n)));
141*9933Ssam }
142*9933Ssam argsq = arg*arg;
143*9933Ssam for(n=0,d=0,i=8;i>=0;i--){
144*9933Ssam n = n*argsq + p1[i];
145*9933Ssam d = d*argsq + q1[i];
146*9933Ssam }
147*9933Ssam return(n/d);
148*9933Ssam }
149*9933Ssam
150*9933Ssam double
y0(arg)151*9933Ssam y0(arg) double arg;{
152*9933Ssam double argsq, n, d;
153*9933Ssam double sin(), cos(), sqrt(), log(), j0();
154*9933Ssam int i;
155*9933Ssam
156*9933Ssam errno = 0;
157*9933Ssam if(arg <= 0.){
158*9933Ssam errno = EDOM;
159*9933Ssam return(-HUGE);
160*9933Ssam }
161*9933Ssam if(arg > 8.){
162*9933Ssam asympt(arg);
163*9933Ssam n = arg - pio4;
164*9933Ssam return(sqrt(tpi/arg)*(pzero*sin(n) + qzero*cos(n)));
165*9933Ssam }
166*9933Ssam argsq = arg*arg;
167*9933Ssam for(n=0,d=0,i=8;i>=0;i--){
168*9933Ssam n = n*argsq + p4[i];
169*9933Ssam d = d*argsq + q4[i];
170*9933Ssam }
171*9933Ssam return(n/d + tpi*j0(arg)*log(arg));
172*9933Ssam }
173*9933Ssam
174*9933Ssam static
asympt(arg)175*9933Ssam asympt(arg) double arg;{
176*9933Ssam double zsq, n, d;
177*9933Ssam int i;
178*9933Ssam zsq = 64./(arg*arg);
179*9933Ssam for(n=0,d=0,i=6;i>=0;i--){
180*9933Ssam n = n*zsq + p2[i];
181*9933Ssam d = d*zsq + q2[i];
182*9933Ssam }
183*9933Ssam pzero = n/d;
184*9933Ssam for(n=0,d=0,i=6;i>=0;i--){
185*9933Ssam n = n*zsq + p3[i];
186*9933Ssam d = d*zsq + q3[i];
187*9933Ssam }
188*9933Ssam qzero = (8./arg)*(n/d);
189*9933Ssam }
190