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