![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to besj0.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10dc / libc / math|
Return to besj0.c CVS log | Up to [Research Unix] / researchv10dc / libc / math |
1.1 root 1: /*
2: floating point Bessel's function
3: of the first and second kinds
4: of order zero
5:
6: besj0(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: besy0(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, besj0.
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: besj0(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: besy0(arg) double arg;{
150: double argsq, n, d;
151: double sin(), cos(), sqrt(), log(), besj0();
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*besj0(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: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.