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