1 #ifndef lint
2 static char sccsid[] = "@(#)j1.c 4.2 (Berkeley) 08/21/85";
3 #endif not lint
4
5 /*
6 floating point Bessel's function
7 of the first and second kinds
8 of order one
9
10 j1(x) returns the value of J1(x)
11 for all real values of x.
12
13 There are no error returns.
14 Calls sin, cos, sqrt.
15
16 There is a niggling bug in J1 which
17 causes errors up to 2e-16 for x in the
18 interval [-8,8].
19 The bug is caused by an inappropriate order
20 of summation of the series. rhm will fix it
21 someday.
22
23 Coefficients are from Hart & Cheney.
24 #6050 (20.98D)
25 #6750 (19.19D)
26 #7150 (19.35D)
27
28 y1(x) returns the value of Y1(x)
29 for positive real values of x.
30 For x<=0, if on the VAX, error number EDOM is set and
31 the reserved operand fault is generated;
32 otherwise (an IEEE machine) an invalid operation is performed.
33
34 Calls sin, cos, sqrt, log, j1.
35
36 The values of Y1 have not been checked
37 to more than ten places.
38
39 Coefficients are from Hart & Cheney.
40 #6447 (22.18D)
41 #6750 (19.19D)
42 #7150 (19.35D)
43 */
44
45 #include <math.h>
46 #ifdef VAX
47 #include <errno.h>
48 #else /* IEEE double */
49 static double zero = 0.e0;
50 #endif
51 static double pzero, qzero;
52 static double tpi = .6366197723675813430755350535e0;
53 static double pio4 = .7853981633974483096156608458e0;
54 static double p1[] = {
55 0.581199354001606143928050809e21,
56 -.6672106568924916298020941484e20,
57 0.2316433580634002297931815435e19,
58 -.3588817569910106050743641413e17,
59 0.2908795263834775409737601689e15,
60 -.1322983480332126453125473247e13,
61 0.3413234182301700539091292655e10,
62 -.4695753530642995859767162166e7,
63 0.2701122710892323414856790990e4,
64 };
65 static double q1[] = {
66 0.1162398708003212287858529400e22,
67 0.1185770712190320999837113348e20,
68 0.6092061398917521746105196863e17,
69 0.2081661221307607351240184229e15,
70 0.5243710262167649715406728642e12,
71 0.1013863514358673989967045588e10,
72 0.1501793594998585505921097578e7,
73 0.1606931573481487801970916749e4,
74 1.0,
75 };
76 static double p2[] = {
77 -.4435757816794127857114720794e7,
78 -.9942246505077641195658377899e7,
79 -.6603373248364939109255245434e7,
80 -.1523529351181137383255105722e7,
81 -.1098240554345934672737413139e6,
82 -.1611616644324610116477412898e4,
83 0.0,
84 };
85 static double q2[] = {
86 -.4435757816794127856828016962e7,
87 -.9934124389934585658967556309e7,
88 -.6585339479723087072826915069e7,
89 -.1511809506634160881644546358e7,
90 -.1072638599110382011903063867e6,
91 -.1455009440190496182453565068e4,
92 1.0,
93 };
94 static double p3[] = {
95 0.3322091340985722351859704442e5,
96 0.8514516067533570196555001171e5,
97 0.6617883658127083517939992166e5,
98 0.1849426287322386679652009819e5,
99 0.1706375429020768002061283546e4,
100 0.3526513384663603218592175580e2,
101 0.0,
102 };
103 static double q3[] = {
104 0.7087128194102874357377502472e6,
105 0.1819458042243997298924553839e7,
106 0.1419460669603720892855755253e7,
107 0.4002944358226697511708610813e6,
108 0.3789022974577220264142952256e5,
109 0.8638367769604990967475517183e3,
110 1.0,
111 };
112 static double p4[] = {
113 -.9963753424306922225996744354e23,
114 0.2655473831434854326894248968e23,
115 -.1212297555414509577913561535e22,
116 0.2193107339917797592111427556e20,
117 -.1965887462722140658820322248e18,
118 0.9569930239921683481121552788e15,
119 -.2580681702194450950541426399e13,
120 0.3639488548124002058278999428e10,
121 -.2108847540133123652824139923e7,
122 0.0,
123 };
124 static double q4[] = {
125 0.5082067366941243245314424152e24,
126 0.5435310377188854170800653097e22,
127 0.2954987935897148674290758119e20,
128 0.1082258259408819552553850180e18,
129 0.2976632125647276729292742282e15,
130 0.6465340881265275571961681500e12,
131 0.1128686837169442121732366891e10,
132 0.1563282754899580604737366452e7,
133 0.1612361029677000859332072312e4,
134 1.0,
135 };
136
137 double
j1(arg)138 j1(arg) double arg;{
139 double xsq, n, d, x;
140 double sin(), cos(), sqrt();
141 int i;
142
143 x = arg;
144 if(x < 0.) x = -x;
145 if(x > 8.){
146 asympt(x);
147 n = x - 3.*pio4;
148 n = sqrt(tpi/x)*(pzero*cos(n) - qzero*sin(n));
149 if(arg <0.) n = -n;
150 return(n);
151 }
152 xsq = x*x;
153 for(n=0,d=0,i=8;i>=0;i--){
154 n = n*xsq + p1[i];
155 d = d*xsq + q1[i];
156 }
157 return(arg*n/d);
158 }
159
160 double
y1(arg)161 y1(arg) double arg;{
162 double xsq, n, d, x;
163 double sin(), cos(), sqrt(), log(), j1();
164 int i;
165
166 x = arg;
167 if(x <= 0.){
168 #ifdef VAX
169 extern double infnan();
170 return(infnan(EDOM)); /* NaN */
171 #else /* IEEE double */
172 return(zero/zero); /* IEEE machines: invalid operation */
173 #endif
174 }
175 if(x > 8.){
176 asympt(x);
177 n = x - 3*pio4;
178 return(sqrt(tpi/x)*(pzero*sin(n) + qzero*cos(n)));
179 }
180 xsq = x*x;
181 for(n=0,d=0,i=9;i>=0;i--){
182 n = n*xsq + p4[i];
183 d = d*xsq + q4[i];
184 }
185 return(x*n/d + tpi*(j1(x)*log(x)-1./x));
186 }
187
188 static
asympt(arg)189 asympt(arg) double arg;{
190 double zsq, n, d;
191 int i;
192 zsq = 64./(arg*arg);
193 for(n=0,d=0,i=6;i>=0;i--){
194 n = n*zsq + p2[i];
195 d = d*zsq + q2[i];
196 }
197 pzero = n/d;
198 for(n=0,d=0,i=6;i>=0;i--){
199 n = n*zsq + p3[i];
200 d = d*zsq + q3[i];
201 }
202 qzero = (8./arg)*(n/d);
203 }
204