![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to j1.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 43BSD / usr.lib / libom|
Return to j1.c CVS log | Up to [CSRG BSD Unix] / 43BSD / usr.lib / libom |
1.1 root 1: /* @(#)j1.c 4.1 12/25/82 */
2:
3: /*
4: floating point Bessel's function
5: of the first and second kinds
6: of order one
7:
8: j1(x) returns the value of J1(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 J1 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: #6050 (20.98D)
23: #6750 (19.19D)
24: #7150 (19.35D)
25:
26: y1(x) returns the value of Y1(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, j1.
32:
33: The values of Y1 have not been checked
34: to more than ten places.
35:
36: Coefficients are from Hart & Cheney.
37: #6447 (22.18D)
38: #6750 (19.19D)
39: #7150 (19.35D)
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.581199354001606143928050809e21,
51: -.6672106568924916298020941484e20,
52: 0.2316433580634002297931815435e19,
53: -.3588817569910106050743641413e17,
54: 0.2908795263834775409737601689e15,
55: -.1322983480332126453125473247e13,
56: 0.3413234182301700539091292655e10,
57: -.4695753530642995859767162166e7,
58: 0.2701122710892323414856790990e4,
59: };
60: static double q1[] = {
61: 0.1162398708003212287858529400e22,
62: 0.1185770712190320999837113348e20,
63: 0.6092061398917521746105196863e17,
64: 0.2081661221307607351240184229e15,
65: 0.5243710262167649715406728642e12,
66: 0.1013863514358673989967045588e10,
67: 0.1501793594998585505921097578e7,
68: 0.1606931573481487801970916749e4,
69: 1.0,
70: };
71: static double p2[] = {
72: -.4435757816794127857114720794e7,
73: -.9942246505077641195658377899e7,
74: -.6603373248364939109255245434e7,
75: -.1523529351181137383255105722e7,
76: -.1098240554345934672737413139e6,
77: -.1611616644324610116477412898e4,
78: 0.0,
79: };
80: static double q2[] = {
81: -.4435757816794127856828016962e7,
82: -.9934124389934585658967556309e7,
83: -.6585339479723087072826915069e7,
84: -.1511809506634160881644546358e7,
85: -.1072638599110382011903063867e6,
86: -.1455009440190496182453565068e4,
87: 1.0,
88: };
89: static double p3[] = {
90: 0.3322091340985722351859704442e5,
91: 0.8514516067533570196555001171e5,
92: 0.6617883658127083517939992166e5,
93: 0.1849426287322386679652009819e5,
94: 0.1706375429020768002061283546e4,
95: 0.3526513384663603218592175580e2,
96: 0.0,
97: };
98: static double q3[] = {
99: 0.7087128194102874357377502472e6,
100: 0.1819458042243997298924553839e7,
101: 0.1419460669603720892855755253e7,
102: 0.4002944358226697511708610813e6,
103: 0.3789022974577220264142952256e5,
104: 0.8638367769604990967475517183e3,
105: 1.0,
106: };
107: static double p4[] = {
108: -.9963753424306922225996744354e23,
109: 0.2655473831434854326894248968e23,
110: -.1212297555414509577913561535e22,
111: 0.2193107339917797592111427556e20,
112: -.1965887462722140658820322248e18,
113: 0.9569930239921683481121552788e15,
114: -.2580681702194450950541426399e13,
115: 0.3639488548124002058278999428e10,
116: -.2108847540133123652824139923e7,
117: 0.0,
118: };
119: static double q4[] = {
120: 0.5082067366941243245314424152e24,
121: 0.5435310377188854170800653097e22,
122: 0.2954987935897148674290758119e20,
123: 0.1082258259408819552553850180e18,
124: 0.2976632125647276729292742282e15,
125: 0.6465340881265275571961681500e12,
126: 0.1128686837169442121732366891e10,
127: 0.1563282754899580604737366452e7,
128: 0.1612361029677000859332072312e4,
129: 1.0,
130: };
131:
132: double
133: j1(arg) double arg;{
134: double xsq, n, d, x;
135: double sin(), cos(), sqrt();
136: int i;
137:
138: x = arg;
139: if(x < 0.) x = -x;
140: if(x > 8.){
141: asympt(x);
142: n = x - 3.*pio4;
143: n = sqrt(tpi/x)*(pzero*cos(n) - qzero*sin(n));
144: if(arg <0.) n = -n;
145: return(n);
146: }
147: xsq = x*x;
148: for(n=0,d=0,i=8;i>=0;i--){
149: n = n*xsq + p1[i];
150: d = d*xsq + q1[i];
151: }
152: return(arg*n/d);
153: }
154:
155: double
156: y1(arg) double arg;{
157: double xsq, n, d, x;
158: double sin(), cos(), sqrt(), log(), j1();
159: int i;
160:
161: errno = 0;
162: x = arg;
163: if(x <= 0.){
164: errno = EDOM;
165: return(-HUGE);
166: }
167: if(x > 8.){
168: asympt(x);
169: n = x - 3*pio4;
170: return(sqrt(tpi/x)*(pzero*sin(n) + qzero*cos(n)));
171: }
172: xsq = x*x;
173: for(n=0,d=0,i=9;i>=0;i--){
174: n = n*xsq + p4[i];
175: d = d*xsq + q4[i];
176: }
177: return(x*n/d + tpi*(j1(x)*log(x)-1./x));
178: }
179:
180: static
181: asympt(arg) double arg;{
182: double zsq, n, d;
183: int i;
184: zsq = 64./(arg*arg);
185: for(n=0,d=0,i=6;i>=0;i--){
186: n = n*zsq + p2[i];
187: d = d*zsq + q2[i];
188: }
189: pzero = n/d;
190: for(n=0,d=0,i=6;i>=0;i--){
191: n = n*zsq + p3[i];
192: d = d*zsq + q3[i];
193: }
194: qzero = (8./arg)*(n/d);
195: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.