Commit | Line | Data |
---|---|---|
71e266b1 SL |
1 | /* @(#)j0.c 4.1 %G% */ |
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 | |
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 | |
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 | |
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 | } |