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