bzr branch
/lh/gsl/trunk
|
2193
by bjg
added gpl headers |
1 |
/* specfunc/bessel_j.c
|
2 |
*
|
|
|
3432
by jungman
probably fixed BUG#3 |
3 |
* Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003 Gerard Jungman
|
|
2193
by bjg
added gpl headers |
4 |
*
|
5 |
* This program is free software; you can redistribute it and/or modify
|
|
6 |
* it under the terms of the GNU General Public License as published by
|
|
|
4212
by bjg
updated to gpl v3 |
7 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
2193
by bjg
added gpl headers |
8 |
* your option) any later version.
|
9 |
*
|
|
10 |
* This program is distributed in the hope that it will be useful, but
|
|
11 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
12 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
13 |
* General Public License for more details.
|
|
14 |
*
|
|
15 |
* You should have received a copy of the GNU General Public License
|
|
16 |
* along with this program; if not, write to the Free Software
|
|
|
3825
by bjg
updated FSF address in copyright notices |
17 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
2193
by bjg
added gpl headers |
18 |
*/
|
19 |
||
|
2952
by bjg
removed auto-expanding RCS tokens from comments as they |
20 |
/* Author: G. Jungman */
|
21 |
||
|
1484
by bjg
added #include <config.h> to all top-level source files |
22 |
#include <config.h> |
|
1997
by jungman
added gsl/ and fixed build to use it |
23 |
#include <gsl/gsl_math.h> |
24 |
#include <gsl/gsl_errno.h> |
|
|
3009
by bjg
made #includes consistent to simplify build |
25 |
#include <gsl/gsl_sf_pow_int.h> |
26 |
#include <gsl/gsl_sf_trig.h> |
|
27 |
#include <gsl/gsl_sf_bessel.h> |
|
|
438
by jungman
added bessel_j.c |
28 |
|
|
2769
by bjg
moved domain, overflow and underflow errors into separate macros, |
29 |
#include "error.h" |
30 |
||
31 |
#include "bessel.h" |
|
32 |
#include "bessel_olver.h" |
|
|
438
by jungman
added bessel_j.c |
33 |
|
|
2530
by bjg
unified error handling across library |
34 |
/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
|
|
438
by jungman
added bessel_j.c |
35 |
|
|
2530
by bjg
unified error handling across library |
36 |
int gsl_sf_bessel_j0_e(const double x, gsl_sf_result * result) |
|
480
by jungman
work in progress |
37 |
{
|
|
1370
by jungman
work in progress |
38 |
double ax = fabs(x); |
39 |
||
|
2530
by bjg
unified error handling across library |
40 |
/* CHECK_POINTER(result) */
|
41 |
||
42 |
if(ax < 0.5) { |
|
|
1520
by jungman
work in progress |
43 |
const double y = x*x; |
44 |
const double c1 = -1.0/6.0; |
|
45 |
const double c2 = 1.0/120.0; |
|
46 |
const double c3 = -1.0/5040.0; |
|
47 |
const double c4 = 1.0/362880.0; |
|
48 |
const double c5 = -1.0/39916800.0; |
|
49 |
const double c6 = 1.0/6227020800.0; |
|
|
1535
by jungman
specfunc broken again -- major changes |
50 |
result->val = 1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*(c5 + y*c6))))); |
|
1562
by jungman
work in progress |
51 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
1370
by jungman
work in progress |
52 |
return GSL_SUCCESS; |
|
480
by jungman
work in progress |
53 |
}
|
54 |
else { |
|
|
4844
by Patrick Alken
partial fix for bug #36152 (bessel j0/j1) |
55 |
result->val = sin(x) / x; |
56 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
57 |
return GSL_SUCCESS; |
|
|
480
by jungman
work in progress |
58 |
}
|
59 |
}
|
|
60 |
||
|
1535
by jungman
specfunc broken again -- major changes |
61 |
|
|
2530
by bjg
unified error handling across library |
62 |
int gsl_sf_bessel_j1_e(const double x, gsl_sf_result * result) |
|
469
by jungman
work in progress |
63 |
{
|
|
1370
by jungman
work in progress |
64 |
double ax = fabs(x); |
65 |
||
|
2530
by bjg
unified error handling across library |
66 |
/* CHECK_POINTER(result) */
|
67 |
||
68 |
if(x == 0.0) { |
|
|
1535
by jungman
specfunc broken again -- major changes |
69 |
result->val = 0.0; |
70 |
result->err = 0.0; |
|
|
1520
by jungman
work in progress |
71 |
return GSL_SUCCESS; |
72 |
}
|
|
|
1858
by jungman
fixed some foo with DBL_MAX, etc |
73 |
else if(ax < 3.1*GSL_DBL_MIN) { |
|
2769
by bjg
moved domain, overflow and underflow errors into separate macros, |
74 |
UNDERFLOW_ERROR(result); |
|
469
by jungman
work in progress |
75 |
}
|
|
1520
by jungman
work in progress |
76 |
else if(ax < 0.25) { |
77 |
const double y = x*x; |
|
78 |
const double c1 = -1.0/10.0; |
|
79 |
const double c2 = 1.0/280.0; |
|
80 |
const double c3 = -1.0/15120.0; |
|
81 |
const double c4 = 1.0/1330560.0; |
|
82 |
const double c5 = -1.0/172972800.0; |
|
83 |
const double sum = 1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*c5)))); |
|
|
1535
by jungman
specfunc broken again -- major changes |
84 |
result->val = x/3.0 * sum; |
|
1562
by jungman
work in progress |
85 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
469
by jungman
work in progress |
86 |
return GSL_SUCCESS; |
87 |
}
|
|
88 |
else { |
|
|
4844
by Patrick Alken
partial fix for bug #36152 (bessel j0/j1) |
89 |
const double cos_x = cos(x); |
90 |
const double sin_x = sin(x); |
|
|
1558
by jungman
work in progress |
91 |
result->val = (sin_x/x - cos_x)/x; |
|
4844
by Patrick Alken
partial fix for bug #36152 (bessel j0/j1) |
92 |
result->err = 2.0 * GSL_DBL_EPSILON * (fabs(sin_x/(x*x)) + fabs(cos_x/x)); |
|
1560
by jungman
work in progress |
93 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
4844
by Patrick Alken
partial fix for bug #36152 (bessel j0/j1) |
94 |
return GSL_SUCCESS; |
|
469
by jungman
work in progress |
95 |
}
|
96 |
}
|
|
97 |
||
|
1535
by jungman
specfunc broken again -- major changes |
98 |
|
|
2530
by bjg
unified error handling across library |
99 |
int gsl_sf_bessel_j2_e(const double x, gsl_sf_result * result) |
|
469
by jungman
work in progress |
100 |
{
|
|
1370
by jungman
work in progress |
101 |
double ax = fabs(x); |
102 |
||
|
2530
by bjg
unified error handling across library |
103 |
/* CHECK_POINTER(result) */
|
104 |
||
105 |
if(x == 0.0) { |
|
|
1535
by jungman
specfunc broken again -- major changes |
106 |
result->val = 0.0; |
107 |
result->err = 0.0; |
|
|
1520
by jungman
work in progress |
108 |
return GSL_SUCCESS; |
109 |
}
|
|
110 |
else if(ax < 4.0*GSL_SQRT_DBL_MIN) { |
|
|
2769
by bjg
moved domain, overflow and underflow errors into separate macros, |
111 |
UNDERFLOW_ERROR(result); |
|
469
by jungman
work in progress |
112 |
}
|
|
1492
by jungman
work in progress |
113 |
else if(ax < 1.3) { |
114 |
const double y = x*x; |
|
115 |
const double c1 = -1.0/14.0; |
|
116 |
const double c2 = 1.0/504.0; |
|
117 |
const double c3 = -1.0/33264.0; |
|
118 |
const double c4 = 1.0/3459456.0; |
|
119 |
const double c5 = -1.0/518918400; |
|
120 |
const double c6 = 1.0/105859353600.0; |
|
121 |
const double c7 = -1.0/28158588057600.0; |
|
122 |
const double c8 = 1.0/9461285587353600.0; |
|
123 |
const double c9 = -1.0/3916972233164390400.0; |
|
|
1520
by jungman
work in progress |
124 |
const double sum = 1.0+y*(c1+y*(c2+y*(c3+y*(c4+y*(c5+y*(c6+y*(c7+y*(c8+y*c9)))))))); |
|
1535
by jungman
specfunc broken again -- major changes |
125 |
result->val = y/15.0 * sum; |
126 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
|
469
by jungman
work in progress |
127 |
return GSL_SUCCESS; |
128 |
}
|
|
129 |
else { |
|
|
1560
by jungman
work in progress |
130 |
gsl_sf_result cos_result; |
131 |
gsl_sf_result sin_result; |
|
|
2530
by bjg
unified error handling across library |
132 |
const int stat_cos = gsl_sf_cos_e(x, &cos_result); |
133 |
const int stat_sin = gsl_sf_sin_e(x, &sin_result); |
|
|
1560
by jungman
work in progress |
134 |
const double cos_x = cos_result.val; |
135 |
const double sin_x = sin_result.val; |
|
136 |
const double f = (3.0/(x*x) - 1.0); |
|
137 |
result->val = (f * sin_x - 3.0*cos_x/x)/x; |
|
138 |
result->err = fabs(f * sin_result.err/x) + fabs((3.0*cos_result.err/x)/x); |
|
|
1744
by jungman
fixed problems exposed by optimization |
139 |
result->err += 2.0 * GSL_DBL_EPSILON * (fabs(f*sin_x/x) + 3.0*fabs(cos_x/(x*x))); |
|
1560
by jungman
work in progress |
140 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
141 |
return GSL_ERROR_SELECT_2(stat_cos, stat_sin); |
|
|
469
by jungman
work in progress |
142 |
}
|
|
1370
by jungman
work in progress |
143 |
}
|
|
469
by jungman
work in progress |
144 |
|
|
1535
by jungman
specfunc broken again -- major changes |
145 |
|
|
1574
by jungman
oh, the painful torture... |
146 |
int
|
|
2530
by bjg
unified error handling across library |
147 |
gsl_sf_bessel_jl_e(const int l, const double x, gsl_sf_result * result) |
|
607
by jungman
cleaned up spherical bessels |
148 |
{
|
149 |
if(l < 0 || x < 0.0) { |
|
|
2769
by bjg
moved domain, overflow and underflow errors into separate macros, |
150 |
DOMAIN_ERROR(result); |
|
607
by jungman
cleaned up spherical bessels |
151 |
}
|
152 |
else if(x == 0.0) { |
|
|
2171
by jungman
fixed jl(0,0) |
153 |
result->val = ( l > 0 ? 0.0 : 1.0 ); |
|
1535
by jungman
specfunc broken again -- major changes |
154 |
result->err = 0.0; |
|
607
by jungman
cleaned up spherical bessels |
155 |
return GSL_SUCCESS; |
156 |
}
|
|
157 |
else if(l == 0) { |
|
|
2530
by bjg
unified error handling across library |
158 |
return gsl_sf_bessel_j0_e(x, result); |
|
607
by jungman
cleaned up spherical bessels |
159 |
}
|
160 |
else if(l == 1) { |
|
|
2530
by bjg
unified error handling across library |
161 |
return gsl_sf_bessel_j1_e(x, result); |
|
607
by jungman
cleaned up spherical bessels |
162 |
}
|
163 |
else if(l == 2) { |
|
|
2530
by bjg
unified error handling across library |
164 |
return gsl_sf_bessel_j2_e(x, result); |
|
607
by jungman
cleaned up spherical bessels |
165 |
}
|
|
1520
by jungman
work in progress |
166 |
else if(x*x < 10.0*(l+0.5)/M_E) { |
|
1576
by jungman
still mighty broken... |
167 |
gsl_sf_result b; |
|
2530
by bjg
unified error handling across library |
168 |
int status = gsl_sf_bessel_IJ_taylor_e(l+0.5, x, -1, 50, GSL_DBL_EPSILON, &b); |
|
1576
by jungman
still mighty broken... |
169 |
double pre = sqrt((0.5*M_PI)/x); |
170 |
result->val = pre * b.val; |
|
171 |
result->err = pre * b.err; |
|
172 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
|
1520
by jungman
work in progress |
173 |
return status; |
174 |
}
|
|
|
3432
by jungman
probably fixed BUG#3 |
175 |
else if(GSL_ROOT4_DBL_EPSILON * x > (l*l + l + 1.0)) { |
|
1536
by jungman
still broken... |
176 |
gsl_sf_result b; |
|
2530
by bjg
unified error handling across library |
177 |
int status = gsl_sf_bessel_Jnu_asympx_e(l + 0.5, x, &b); |
|
1536
by jungman
still broken... |
178 |
double pre = sqrt((0.5*M_PI)/x); |
179 |
result->val = pre * b.val; |
|
|
1562
by jungman
work in progress |
180 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + pre * b.err; |
|
1520
by jungman
work in progress |
181 |
return status; |
182 |
}
|
|
183 |
else if(l > 1.0/GSL_ROOT6_DBL_EPSILON) { |
|
|
1536
by jungman
still broken... |
184 |
gsl_sf_result b; |
|
2530
by bjg
unified error handling across library |
185 |
int status = gsl_sf_bessel_Jnu_asymp_Olver_e(l + 0.5, x, &b); |
|
1536
by jungman
still broken... |
186 |
double pre = sqrt((0.5*M_PI)/x); |
187 |
result->val = pre * b.val; |
|
|
1562
by jungman
work in progress |
188 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + pre * b.err; |
|
1520
by jungman
work in progress |
189 |
return status; |
190 |
}
|
|
|
4384
by Brian Gough
use full asymptotic expansion for gsl_sf_bessel_Jnu_asympx_e, fixes bug #21819 |
191 |
else if(x > 1000.0 && x > l*l) |
|
3432
by jungman
probably fixed BUG#3 |
192 |
{
|
|
4384
by Brian Gough
use full asymptotic expansion for gsl_sf_bessel_Jnu_asympx_e, fixes bug #21819 |
193 |
/* We need this path to avoid feeding large x to CF1 below; */
|
|
3432
by jungman
probably fixed BUG#3 |
194 |
gsl_sf_result b; |
195 |
int status = gsl_sf_bessel_Jnu_asympx_e(l + 0.5, x, &b); |
|
196 |
double pre = sqrt((0.5*M_PI)/x); |
|
197 |
result->val = pre * b.val; |
|
198 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + pre * b.err; |
|
199 |
return status; |
|
200 |
}
|
|
|
607
by jungman
cleaned up spherical bessels |
201 |
else { |
|
1576
by jungman
still mighty broken... |
202 |
double sgn; |
|
1574
by jungman
oh, the painful torture... |
203 |
double ratio; |
|
4384
by Brian Gough
use full asymptotic expansion for gsl_sf_bessel_Jnu_asympx_e, fixes bug #21819 |
204 |
/* The CF1 call will hit 10000 iterations for x > 10000 + l */
|
|
1576
by jungman
still mighty broken... |
205 |
int stat_CF1 = gsl_sf_bessel_J_CF1(l+0.5, x, &ratio, &sgn); |
|
4850
by Patrick Alken
partial fix for bug #37209 + test case |
206 |
const double BESSEL_J_SMALL = GSL_DBL_MIN / GSL_DBL_EPSILON; |
207 |
double jellp1 = BESSEL_J_SMALL * ratio; |
|
208 |
double jell = BESSEL_J_SMALL; |
|
|
1536
by jungman
still broken... |
209 |
double jellm1; |
|
607
by jungman
cleaned up spherical bessels |
210 |
int ell; |
|
1574
by jungman
oh, the painful torture... |
211 |
for(ell = l; ell > 0; ell--) { |
|
607
by jungman
cleaned up spherical bessels |
212 |
jellm1 = -jellp1 + (2*ell + 1)/x * jell; |
213 |
jellp1 = jell; |
|
214 |
jell = jellm1; |
|
215 |
}
|
|
|
1520
by jungman
work in progress |
216 |
|
|
1574
by jungman
oh, the painful torture... |
217 |
if(fabs(jell) > fabs(jellp1)) { |
218 |
gsl_sf_result j0_result; |
|
|
2530
by bjg
unified error handling across library |
219 |
int stat_j0 = gsl_sf_bessel_j0_e(x, &j0_result); |
|
4850
by Patrick Alken
partial fix for bug #37209 + test case |
220 |
double pre = BESSEL_J_SMALL / jell; |
|
1574
by jungman
oh, the painful torture... |
221 |
result->val = j0_result.val * pre; |
222 |
result->err = j0_result.err * fabs(pre); |
|
|
4382
by Brian Gough
increase error estimate |
223 |
result->err += 4.0 * GSL_DBL_EPSILON * (0.5*l + 1.0) * fabs(result->val); |
|
1574
by jungman
oh, the painful torture... |
224 |
return GSL_ERROR_SELECT_2(stat_j0, stat_CF1); |
225 |
}
|
|
226 |
else { |
|
227 |
gsl_sf_result j1_result; |
|
|
2530
by bjg
unified error handling across library |
228 |
int stat_j1 = gsl_sf_bessel_j1_e(x, &j1_result); |
|
4850
by Patrick Alken
partial fix for bug #37209 + test case |
229 |
double pre = BESSEL_J_SMALL / jellp1; |
|
1574
by jungman
oh, the painful torture... |
230 |
result->val = j1_result.val * pre; |
231 |
result->err = j1_result.err * fabs(pre); |
|
|
4382
by Brian Gough
increase error estimate |
232 |
result->err += 4.0 * GSL_DBL_EPSILON * (0.5*l + 1.0) * fabs(result->val); |
|
1574
by jungman
oh, the painful torture... |
233 |
return GSL_ERROR_SELECT_2(stat_j1, stat_CF1); |
234 |
}
|
|
|
607
by jungman
cleaned up spherical bessels |
235 |
}
|
236 |
}
|
|
237 |
||
238 |
||
|
1574
by jungman
oh, the painful torture... |
239 |
int
|
|
2530
by bjg
unified error handling across library |
240 |
gsl_sf_bessel_jl_array(const int lmax, const double x, double * result_array) |
|
607
by jungman
cleaned up spherical bessels |
241 |
{
|
|
2530
by bjg
unified error handling across library |
242 |
/* CHECK_POINTER(result_array) */
|
243 |
||
244 |
if(lmax < 0 || x < 0.0) { |
|
|
1370
by jungman
work in progress |
245 |
int j; |
246 |
for(j=0; j<=lmax; j++) result_array[j] = 0.0; |
|
|
2530
by bjg
unified error handling across library |
247 |
GSL_ERROR ("error", GSL_EDOM); |
|
1370
by jungman
work in progress |
248 |
}
|
|
2171
by jungman
fixed jl(0,0) |
249 |
else if(x == 0.0) { |
250 |
int j; |
|
251 |
for(j=1; j<=lmax; j++) result_array[j] = 0.0; |
|
252 |
result_array[0] = 1.0; |
|
253 |
return GSL_SUCCESS; |
|
254 |
}
|
|
|
1370
by jungman
work in progress |
255 |
else { |
|
1536
by jungman
still broken... |
256 |
gsl_sf_result r_jellp1; |
257 |
gsl_sf_result r_jell; |
|
|
2530
by bjg
unified error handling across library |
258 |
int stat_0 = gsl_sf_bessel_jl_e(lmax+1, x, &r_jellp1); |
259 |
int stat_1 = gsl_sf_bessel_jl_e(lmax, x, &r_jell); |
|
|
1536
by jungman
still broken... |
260 |
double jellp1 = r_jellp1.val; |
261 |
double jell = r_jell.val; |
|
262 |
double jellm1; |
|
|
1370
by jungman
work in progress |
263 |
int ell; |
|
1536
by jungman
still broken... |
264 |
|
|
1370
by jungman
work in progress |
265 |
result_array[lmax] = jell; |
266 |
for(ell = lmax; ell >= 1; ell--) { |
|
267 |
jellm1 = -jellp1 + (2*ell + 1)/x * jell; |
|
268 |
jellp1 = jell; |
|
269 |
jell = jellm1; |
|
270 |
result_array[ell-1] = jellm1; |
|
271 |
}
|
|
|
1536
by jungman
still broken... |
272 |
|
|
1370
by jungman
work in progress |
273 |
return GSL_ERROR_SELECT_2(stat_0, stat_1); |
274 |
}
|
|
|
607
by jungman
cleaned up spherical bessels |
275 |
}
|
276 |
||
|
1536
by jungman
still broken... |
277 |
|
|
2530
by bjg
unified error handling across library |
278 |
int gsl_sf_bessel_jl_steed_array(const int lmax, const double x, double * jl_x) |
|
438
by jungman
added bessel_j.c |
279 |
{
|
|
2530
by bjg
unified error handling across library |
280 |
/* CHECK_POINTER(jl_x) */
|
281 |
||
282 |
if(lmax < 0 || x < 0.0) { |
|
|
1370
by jungman
work in progress |
283 |
int j; |
284 |
for(j=0; j<=lmax; j++) jl_x[j] = 0.0; |
|
|
2530
by bjg
unified error handling across library |
285 |
GSL_ERROR ("error", GSL_EDOM); |
|
438
by jungman
added bessel_j.c |
286 |
}
|
|
2171
by jungman
fixed jl(0,0) |
287 |
else if(x == 0.0) { |
288 |
int j; |
|
289 |
for(j=1; j<=lmax; j++) jl_x[j] = 0.0; |
|
290 |
jl_x[0] = 1.0; |
|
291 |
return GSL_SUCCESS; |
|
292 |
}
|
|
|
1536
by jungman
still broken... |
293 |
else if(x < 2.0*GSL_ROOT4_DBL_EPSILON) { |
|
438
by jungman
added bessel_j.c |
294 |
/* first two terms of Taylor series */
|
|
1370
by jungman
work in progress |
295 |
double inv_fact = 1.0; /* 1/(1 3 5 ... (2l+1)) */ |
296 |
double x_l = 1.0; /* x^l */ |
|
|
438
by jungman
added bessel_j.c |
297 |
int l; |
298 |
for(l=0; l<=lmax; l++) { |
|
299 |
jl_x[l] = x_l * inv_fact; |
|
|
1370
by jungman
work in progress |
300 |
jl_x[l] *= 1.0 - 0.5*x*x/(2.0*l+3.0); |
|
1520
by jungman
work in progress |
301 |
inv_fact /= 2.0*l+3.0; |
|
438
by jungman
added bessel_j.c |
302 |
x_l *= x; |
303 |
}
|
|
|
1536
by jungman
still broken... |
304 |
return GSL_SUCCESS; |
|
438
by jungman
added bessel_j.c |
305 |
}
|
306 |
else { |
|
|
607
by jungman
cleaned up spherical bessels |
307 |
/* Steed/Barnett algorithm [Comp. Phys. Comm. 21, 297 (1981)] */
|
|
1370
by jungman
work in progress |
308 |
double x_inv = 1.0/x; |
309 |
double W = 2.0*x_inv; |
|
310 |
double F = 1.0; |
|
311 |
double FP = (lmax+1.0) * x_inv; |
|
312 |
double B = 2.0*FP + x_inv; |
|
313 |
double end = B + 20000.0*W; |
|
314 |
double D = 1.0/B; |
|
|
438
by jungman
added bessel_j.c |
315 |
double del = -D; |
316 |
||
317 |
FP += del; |
|
318 |
||
319 |
/* continued fraction */
|
|
320 |
do { |
|
321 |
B += W; |
|
|
1370
by jungman
work in progress |
322 |
D = 1.0/(B-D); |
|
438
by jungman
added bessel_j.c |
323 |
del *= (B*D - 1.); |
324 |
FP += del; |
|
|
1370
by jungman
work in progress |
325 |
if(D < 0.0) F = -F; |
|
438
by jungman
added bessel_j.c |
326 |
if(B > end) { |
|
3468
by bjg
detabification |
327 |
GSL_ERROR ("error", GSL_EMAXITER); |
|
438
by jungman
added bessel_j.c |
328 |
}
|
329 |
}
|
|
|
1520
by jungman
work in progress |
330 |
while(fabs(del) >= fabs(FP) * GSL_DBL_EPSILON); |
|
438
by jungman
added bessel_j.c |
331 |
|
332 |
FP *= F; |
|
333 |
||
334 |
if(lmax > 0) { |
|
335 |
/* downward recursion */
|
|
336 |
double XP2 = FP; |
|
337 |
double PL = lmax * x_inv; |
|
338 |
int L = lmax; |
|
339 |
int LP; |
|
340 |
jl_x[lmax] = F; |
|
341 |
for(LP = 1; LP<=lmax; LP++) { |
|
|
3468
by bjg
detabification |
342 |
jl_x[L-1] = PL * jl_x[L] + XP2; |
343 |
FP = PL*jl_x[L-1] - jl_x[L]; |
|
344 |
XP2 = FP; |
|
345 |
PL -= x_inv; |
|
346 |
--L; |
|
|
438
by jungman
added bessel_j.c |
347 |
}
|
348 |
F = jl_x[0]; |
|
349 |
}
|
|
350 |
||
351 |
/* normalization */
|
|
|
3804
by bjg
use hypot internally to avoid overflow |
352 |
W = x_inv / hypot(FP, F); |
|
438
by jungman
added bessel_j.c |
353 |
jl_x[0] = W*F; |
354 |
if(lmax > 0) { |
|
355 |
int L; |
|
356 |
for(L=1; L<=lmax; L++) { |
|
|
3468
by bjg
detabification |
357 |
jl_x[L] *= W; |
|
438
by jungman
added bessel_j.c |
358 |
}
|
359 |
}
|
|
|
1536
by jungman
still broken... |
360 |
|
361 |
return GSL_SUCCESS; |
|
|
438
by jungman
added bessel_j.c |
362 |
}
|
363 |
}
|
|
364 |
||
365 |
||
|
2530
by bjg
unified error handling across library |
366 |
/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
|
367 |
||
368 |
#include "eval.h" |
|
369 |
||
370 |
double gsl_sf_bessel_j0(const double x) |
|
371 |
{
|
|
372 |
EVAL_RESULT(gsl_sf_bessel_j0_e(x, &result)); |
|
373 |
}
|
|
374 |
||
375 |
double gsl_sf_bessel_j1(const double x) |
|
376 |
{
|
|
377 |
EVAL_RESULT(gsl_sf_bessel_j1_e(x, &result)); |
|
378 |
}
|
|
379 |
||
380 |
double gsl_sf_bessel_j2(const double x) |
|
381 |
{
|
|
382 |
EVAL_RESULT(gsl_sf_bessel_j2_e(x, &result)); |
|
383 |
}
|
|
384 |
||
385 |
double gsl_sf_bessel_jl(const int l, const double x) |
|
386 |
{
|
|
387 |
EVAL_RESULT(gsl_sf_bessel_jl_e(l, x, &result)); |
|
388 |
}
|
|
389 |