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