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