```   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: }
```

#### Defined functions

asympt defined in line 180; used 2 times
j1 defined in line 132; used 11 times

#### Defined variables

errno defined in line 45; used 2 times
p1 defined in line 49; used 1 times
p2 defined in line 71; used 1 times
p3 defined in line 89; used 1 times
p4 defined in line 107; used 1 times
pio4 defined in line 48; used 2 times
pzero defined in line 46; used 3 times
q1 defined in line 60; used 1 times
q2 defined in line 80; used 1 times
q3 defined in line 98; used 1 times
q4 defined in line 119; used 1 times
qzero defined in line 46; used 3 times
tpi defined in line 47; used 3 times
 Last modified: 1985-06-05 Generated: 2016-12-26 Generated by src2html V0.67 page hit count: 2648