Line data Source code
1 : /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2 : /*
3 : * This file is part of the LibreOffice project.
4 : *
5 : * This Source Code Form is subject to the terms of the Mozilla Public
6 : * License, v. 2.0. If a copy of the MPL was not distributed with this
7 : * file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 : *
9 : * This file incorporates work covered by the following license notice:
10 : *
11 : * Licensed to the Apache Software Foundation (ASF) under one or more
12 : * contributor license agreements. See the NOTICE file distributed
13 : * with this work for additional information regarding copyright
14 : * ownership. The ASF licenses this file to you under the Apache
15 : * License, Version 2.0 (the "License"); you may not use this file
16 : * except in compliance with the License. You may obtain a copy of
17 : * the License at http://www.apache.org/licenses/LICENSE-2.0 .
18 : */
19 :
20 : #include "bessel.hxx"
21 : #include "analysishelper.hxx"
22 :
23 : #include <rtl/math.hxx>
24 :
25 : using ::com::sun::star::lang::IllegalArgumentException;
26 : using ::com::sun::star::sheet::NoConvergenceException;
27 :
28 : namespace sca {
29 : namespace analysis {
30 :
31 : const double f_PI = 3.1415926535897932385;
32 : const double f_PI_DIV_2 = f_PI / 2.0;
33 : const double f_PI_DIV_4 = f_PI / 4.0;
34 : const double f_2_DIV_PI = 2.0 / f_PI;
35 :
36 :
37 : // BESSEL J
38 :
39 :
40 : /* The BESSEL function, first kind, unmodified:
41 : The algorithm follows
42 : http://www.reference-global.com/isbn/978-3-11-020354-7
43 : Numerical Mathematics 1 / Numerische Mathematik 1,
44 : An algorithm-based introduction / Eine algorithmisch orientierte Einfuehrung
45 : Deuflhard, Peter; Hohmann, Andreas
46 : Berlin, New York (Walter de Gruyter) 2008
47 : 4. ueberarb. u. erw. Aufl. 2008
48 : eBook ISBN: 978-3-11-020355-4
49 : Chapter 6.3.2 , algorithm 6.24
50 : The source is in German.
51 : The BesselJ-function is a special case of the adjoint summation with
52 : a_k = 2*(k-1)/x for k=1,...
53 : b_k = -1, for all k, directly substituted
54 : m_0=1, m_k=2 for k even, and m_k=0 for k odd, calculated on the fly
55 : alpha_k=1 for k=N and alpha_k=0 otherwise
56 : */
57 :
58 0 : double BesselJ( double x, sal_Int32 N ) throw (IllegalArgumentException, NoConvergenceException)
59 :
60 : {
61 0 : if( N < 0 )
62 0 : throw IllegalArgumentException();
63 0 : if (x==0.0)
64 0 : return (N==0) ? 1.0 : 0.0;
65 :
66 : /* The algorithm works only for x>0, therefore remember sign. BesselJ
67 : with integer order N is an even function for even N (means J(-x)=J(x))
68 : and an odd function for odd N (means J(-x)=-J(x)).*/
69 0 : double fSign = (N % 2 == 1 && x < 0) ? -1.0 : 1.0;
70 0 : double fX = fabs(x);
71 :
72 0 : const double fMaxIteration = 9000000.0; //experimental, for to return in < 3 seconds
73 0 : double fEstimateIteration = fX * 1.5 + N;
74 0 : bool bAsymptoticPossible = pow(fX,0.4) > N;
75 0 : if (fEstimateIteration > fMaxIteration)
76 : {
77 0 : if (bAsymptoticPossible)
78 0 : return fSign * sqrt(f_2_DIV_PI/fX)* cos(fX-N*f_PI_DIV_2-f_PI_DIV_4);
79 : else
80 0 : throw NoConvergenceException();
81 : }
82 :
83 0 : double epsilon = 1.0e-15; // relative error
84 0 : bool bHasfound = false;
85 0 : double k= 0.0;
86 : // e_{-1} = 0; e_0 = alpha_0 / b_2
87 : double u ; // u_0 = e_0/f_0 = alpha_0/m_0 = alpha_0
88 :
89 : // first used with k=1
90 : double m_bar; // m_bar_k = m_k * f_bar_{k-1}
91 : double g_bar; // g_bar_k = m_bar_k - a_{k+1} + g_{k-1}
92 : double g_bar_delta_u; // g_bar_delta_u_k = f_bar_{k-1} * alpha_k
93 : // - g_{k-1} * delta_u_{k-1} - m_bar_k * u_{k-1}
94 : // f_{-1} = 0.0; f_0 = m_0 / b_2 = 1/(-1) = -1
95 0 : double g = 0.0; // g_0= f_{-1} / f_0 = 0/(-1) = 0
96 0 : double delta_u = 0.0; // dummy initialize, first used with * 0
97 0 : double f_bar = -1.0; // f_bar_k = 1/f_k, but only used for k=0
98 :
99 0 : if (N==0)
100 : {
101 : //k=0; alpha_0 = 1.0
102 0 : u = 1.0; // u_0 = alpha_0
103 : // k = 1.0; at least one step is necessary
104 : // m_bar_k = m_k * f_bar_{k-1} ==> m_bar_1 = 0.0
105 0 : g_bar_delta_u = 0.0; // alpha_k = 0.0, m_bar = 0.0; g= 0.0
106 0 : g_bar = - 2.0/fX; // k = 1.0, g = 0.0
107 0 : delta_u = g_bar_delta_u / g_bar;
108 0 : u = u + delta_u ; // u_k = u_{k-1} + delta_u_k
109 0 : g = -1.0 / g_bar; // g_k=b_{k+2}/g_bar_k
110 0 : f_bar = f_bar * g; // f_bar_k = f_bar_{k-1}* g_k
111 0 : k = 2.0;
112 : // From now on all alpha_k = 0.0 and k > N+1
113 : }
114 : else
115 : { // N >= 1 and alpha_k = 0.0 for k<N
116 0 : u=0.0; // u_0 = alpha_0
117 0 : for (k =1.0; k<= N-1; k = k + 1.0)
118 : {
119 0 : m_bar=2.0 * fmod(k-1.0, 2.0) * f_bar;
120 0 : g_bar_delta_u = - g * delta_u - m_bar * u; // alpha_k = 0.0
121 0 : g_bar = m_bar - 2.0*k/fX + g;
122 0 : delta_u = g_bar_delta_u / g_bar;
123 0 : u = u + delta_u;
124 0 : g = -1.0/g_bar;
125 0 : f_bar=f_bar * g;
126 : }
127 : // Step alpha_N = 1.0
128 0 : m_bar=2.0 * fmod(k-1.0, 2.0) * f_bar;
129 0 : g_bar_delta_u = f_bar - g * delta_u - m_bar * u; // alpha_k = 1.0
130 0 : g_bar = m_bar - 2.0*k/fX + g;
131 0 : delta_u = g_bar_delta_u / g_bar;
132 0 : u = u + delta_u;
133 0 : g = -1.0/g_bar;
134 0 : f_bar = f_bar * g;
135 0 : k = k + 1.0;
136 : }
137 : // Loop until desired accuracy, always alpha_k = 0.0
138 0 : do
139 : {
140 0 : m_bar = 2.0 * fmod(k-1.0, 2.0) * f_bar;
141 0 : g_bar_delta_u = - g * delta_u - m_bar * u;
142 0 : g_bar = m_bar - 2.0*k/fX + g;
143 0 : delta_u = g_bar_delta_u / g_bar;
144 0 : u = u + delta_u;
145 0 : g = -1.0/g_bar;
146 0 : f_bar = f_bar * g;
147 0 : bHasfound = (fabs(delta_u)<=fabs(u)*epsilon);
148 0 : k = k + 1.0;
149 : }
150 0 : while (!bHasfound && k <= fMaxIteration);
151 0 : if (bHasfound)
152 0 : return u * fSign;
153 : else
154 0 : throw NoConvergenceException(); // unlikely to happen
155 : }
156 :
157 :
158 : // BESSEL I
159 :
160 :
161 : /* The BESSEL function, first kind, modified:
162 :
163 : inf (x/2)^(n+2k)
164 : I_n(x) = SUM TERM(n,k) with TERM(n,k) := --------------
165 : k=0 k! (n+k)!
166 :
167 : No asymptotic approximation used, see issue 43040.
168 : */
169 :
170 0 : double BesselI( double x, sal_Int32 n ) throw( IllegalArgumentException, NoConvergenceException )
171 : {
172 0 : const sal_Int32 nMaxIteration = 2000;
173 0 : const double fXHalf = x / 2.0;
174 0 : if( n < 0 )
175 0 : throw IllegalArgumentException();
176 :
177 0 : double fResult = 0.0;
178 :
179 : /* Start the iteration without TERM(n,0), which is set here.
180 :
181 : TERM(n,0) = (x/2)^n / n!
182 : */
183 0 : sal_Int32 nK = 0;
184 0 : double fTerm = 1.0;
185 : // avoid overflow in Fak(n)
186 0 : for( nK = 1; nK <= n; ++nK )
187 : {
188 0 : fTerm = fTerm / static_cast< double >( nK ) * fXHalf;
189 : }
190 0 : fResult = fTerm; // Start result with TERM(n,0).
191 0 : if( fTerm != 0.0 )
192 : {
193 0 : nK = 1;
194 0 : const double fEpsilon = 1.0E-15;
195 0 : do
196 : {
197 : /* Calculation of TERM(n,k) from TERM(n,k-1):
198 :
199 : (x/2)^(n+2k)
200 : TERM(n,k) = --------------
201 : k! (n+k)!
202 :
203 : (x/2)^2 (x/2)^(n+2(k-1))
204 : = --------------------------
205 : k (k-1)! (n+k) (n+k-1)!
206 :
207 : (x/2)^2 (x/2)^(n+2(k-1))
208 : = --------- * ------------------
209 : k(n+k) (k-1)! (n+k-1)!
210 :
211 : x^2/4
212 : = -------- TERM(n,k-1)
213 : k(n+k)
214 : */
215 0 : fTerm = fTerm * fXHalf / static_cast<double>(nK) * fXHalf / static_cast<double>(nK+n);
216 0 : fResult += fTerm;
217 0 : nK++;
218 : }
219 0 : while( (fabs( fTerm ) > fabs(fResult) * fEpsilon) && (nK < nMaxIteration) );
220 :
221 : }
222 0 : return fResult;
223 : }
224 :
225 :
226 0 : double Besselk0( double fNum ) throw( IllegalArgumentException, NoConvergenceException )
227 : {
228 : double fRet;
229 :
230 0 : if( fNum <= 2.0 )
231 : {
232 0 : double fNum2 = fNum * 0.5;
233 0 : double y = fNum2 * fNum2;
234 :
235 0 : fRet = -log( fNum2 ) * BesselI( fNum, 0 ) +
236 0 : ( -0.57721566 + y * ( 0.42278420 + y * ( 0.23069756 + y * ( 0.3488590e-1 +
237 0 : y * ( 0.262698e-2 + y * ( 0.10750e-3 + y * 0.74e-5 ) ) ) ) ) );
238 : }
239 : else
240 : {
241 0 : double y = 2.0 / fNum;
242 :
243 0 : fRet = exp( -fNum ) / sqrt( fNum ) * ( 1.25331414 + y * ( -0.7832358e-1 +
244 0 : y * ( 0.2189568e-1 + y * ( -0.1062446e-1 + y * ( 0.587872e-2 +
245 0 : y * ( -0.251540e-2 + y * 0.53208e-3 ) ) ) ) ) );
246 : }
247 :
248 0 : return fRet;
249 : }
250 :
251 :
252 0 : double Besselk1( double fNum ) throw( IllegalArgumentException, NoConvergenceException )
253 : {
254 : double fRet;
255 :
256 0 : if( fNum <= 2.0 )
257 : {
258 0 : double fNum2 = fNum * 0.5;
259 0 : double y = fNum2 * fNum2;
260 :
261 0 : fRet = log( fNum2 ) * BesselI( fNum, 1 ) +
262 0 : ( 1.0 + y * ( 0.15443144 + y * ( -0.67278579 + y * ( -0.18156897 + y * ( -0.1919402e-1 +
263 0 : y * ( -0.110404e-2 + y * ( -0.4686e-4 ) ) ) ) ) ) )
264 0 : / fNum;
265 : }
266 : else
267 : {
268 0 : double y = 2.0 / fNum;
269 :
270 0 : fRet = exp( -fNum ) / sqrt( fNum ) * ( 1.25331414 + y * ( 0.23498619 +
271 0 : y * ( -0.3655620e-1 + y * ( 0.1504268e-1 + y * ( -0.780353e-2 +
272 0 : y * ( 0.325614e-2 + y * ( -0.68245e-3 ) ) ) ) ) ) );
273 : }
274 :
275 0 : return fRet;
276 : }
277 :
278 :
279 0 : double BesselK( double fNum, sal_Int32 nOrder ) throw( IllegalArgumentException, NoConvergenceException )
280 : {
281 0 : switch( nOrder )
282 : {
283 0 : case 0: return Besselk0( fNum );
284 0 : case 1: return Besselk1( fNum );
285 : default:
286 : {
287 0 : double fTox = 2.0 / fNum;
288 0 : double fBkm = Besselk0( fNum );
289 0 : double fBk = Besselk1( fNum );
290 :
291 0 : for( sal_Int32 n = 1 ; n < nOrder ; n++ )
292 : {
293 0 : const double fBkp = fBkm + double( n ) * fTox * fBk;
294 0 : fBkm = fBk;
295 0 : fBk = fBkp;
296 : }
297 :
298 0 : return fBk;
299 : }
300 : }
301 : }
302 :
303 :
304 : // BESSEL Y
305 :
306 :
307 : /* The BESSEL function, second kind, unmodified:
308 : The algorithm for order 0 and for order 1 follows
309 : http://www.reference-global.com/isbn/978-3-11-020354-7
310 : Numerical Mathematics 1 / Numerische Mathematik 1,
311 : An algorithm-based introduction / Eine algorithmisch orientierte Einfuehrung
312 : Deuflhard, Peter; Hohmann, Andreas
313 : Berlin, New York (Walter de Gruyter) 2008
314 : 4. ueberarb. u. erw. Aufl. 2008
315 : eBook ISBN: 978-3-11-020355-4
316 : Chapter 6.3.2 , algorithm 6.24
317 : The source is in German.
318 : See #i31656# for a commented version of the implementation, attachment #desc6
319 : http://www.openoffice.org/nonav/issues/showattachment.cgi/63609/Comments%20to%20the%20implementation%20of%20the%20Bessel%20functions.odt
320 : */
321 :
322 0 : double Bessely0( double fX ) throw( IllegalArgumentException, NoConvergenceException )
323 : {
324 0 : if (fX <= 0)
325 0 : throw IllegalArgumentException();
326 0 : const double fMaxIteration = 9000000.0; // should not be reached
327 0 : if (fX > 5.0e+6) // iteration is not considerable better then approximation
328 0 : return sqrt(1/f_PI/fX)
329 0 : *(rtl::math::sin(fX)-rtl::math::cos(fX));
330 0 : const double epsilon = 1.0e-15;
331 0 : const double EulerGamma = 0.57721566490153286060;
332 0 : double alpha = log(fX/2.0)+EulerGamma;
333 0 : double u = alpha;
334 :
335 0 : double k = 1.0;
336 0 : double m_bar = 0.0;
337 0 : double g_bar_delta_u = 0.0;
338 0 : double g_bar = -2.0 / fX;
339 0 : double delta_u = g_bar_delta_u / g_bar;
340 0 : double g = -1.0/g_bar;
341 0 : double f_bar = -1 * g;
342 :
343 0 : double sign_alpha = 1.0;
344 : double km1mod2;
345 0 : bool bHasFound = false;
346 0 : k = k + 1;
347 0 : do
348 : {
349 0 : km1mod2 = fmod(k-1.0,2.0);
350 0 : m_bar=(2.0*km1mod2) * f_bar;
351 0 : if (km1mod2 == 0.0)
352 0 : alpha = 0.0;
353 : else
354 : {
355 0 : alpha = sign_alpha * (4.0/k);
356 0 : sign_alpha = -sign_alpha;
357 : }
358 0 : g_bar_delta_u = f_bar * alpha - g * delta_u - m_bar * u;
359 0 : g_bar = m_bar - (2.0*k)/fX + g;
360 0 : delta_u = g_bar_delta_u / g_bar;
361 0 : u = u+delta_u;
362 0 : g = -1.0 / g_bar;
363 0 : f_bar = f_bar*g;
364 0 : bHasFound = (fabs(delta_u)<=fabs(u)*epsilon);
365 0 : k=k+1;
366 : }
367 0 : while (!bHasFound && k<fMaxIteration);
368 0 : if (bHasFound)
369 0 : return u*f_2_DIV_PI;
370 : else
371 0 : throw NoConvergenceException(); // not likely to happen
372 : }
373 :
374 : // See #i31656# for a commented version of this implementation, attachment #desc6
375 : // http://www.openoffice.org/nonav/issues/showattachment.cgi/63609/Comments%20to%20the%20implementation%20of%20the%20Bessel%20functions.odt
376 0 : double Bessely1( double fX ) throw( IllegalArgumentException, NoConvergenceException )
377 : {
378 0 : if (fX <= 0)
379 0 : throw IllegalArgumentException();
380 0 : const double fMaxIteration = 9000000.0; // should not be reached
381 0 : if (fX > 5.0e+6) // iteration is not considerable better then approximation
382 0 : return - sqrt(1/f_PI/fX)
383 0 : *(rtl::math::sin(fX)+rtl::math::cos(fX));
384 0 : const double epsilon = 1.0e-15;
385 0 : const double EulerGamma = 0.57721566490153286060;
386 0 : double alpha = 1.0/fX;
387 0 : double f_bar = -1.0;
388 0 : double u = alpha;
389 0 : double k = 1.0;
390 0 : double m_bar = 0.0;
391 0 : alpha = 1.0 - EulerGamma - log(fX/2.0);
392 0 : double g_bar_delta_u = -alpha;
393 0 : double g_bar = -2.0 / fX;
394 0 : double delta_u = g_bar_delta_u / g_bar;
395 0 : u = u + delta_u;
396 0 : double g = -1.0/g_bar;
397 0 : f_bar = f_bar * g;
398 0 : double sign_alpha = -1.0;
399 : double km1mod2; //will be (k-1) mod 2
400 : double q; // will be (k-1) div 2
401 0 : bool bHasFound = false;
402 0 : k = k + 1.0;
403 0 : do
404 : {
405 0 : km1mod2 = fmod(k-1.0,2.0);
406 0 : m_bar=(2.0*km1mod2) * f_bar;
407 0 : q = (k-1.0)/2.0;
408 0 : if (km1mod2 == 0.0) // k is odd
409 : {
410 0 : alpha = sign_alpha * (1.0/q + 1.0/(q+1.0));
411 0 : sign_alpha = -sign_alpha;
412 : }
413 : else
414 0 : alpha = 0.0;
415 0 : g_bar_delta_u = f_bar * alpha - g * delta_u - m_bar * u;
416 0 : g_bar = m_bar - (2.0*k)/fX + g;
417 0 : delta_u = g_bar_delta_u / g_bar;
418 0 : u = u+delta_u;
419 0 : g = -1.0 / g_bar;
420 0 : f_bar = f_bar*g;
421 0 : bHasFound = (fabs(delta_u)<=fabs(u)*epsilon);
422 0 : k=k+1;
423 : }
424 0 : while (!bHasFound && k<fMaxIteration);
425 0 : if (bHasFound)
426 0 : return -u*2.0/f_PI;
427 : else
428 0 : throw NoConvergenceException();
429 : }
430 :
431 0 : double BesselY( double fNum, sal_Int32 nOrder ) throw( IllegalArgumentException, NoConvergenceException )
432 : {
433 0 : switch( nOrder )
434 : {
435 0 : case 0: return Bessely0( fNum );
436 0 : case 1: return Bessely1( fNum );
437 : default:
438 : {
439 0 : double fTox = 2.0 / fNum;
440 0 : double fBym = Bessely0( fNum );
441 0 : double fBy = Bessely1( fNum );
442 :
443 0 : for( sal_Int32 n = 1 ; n < nOrder ; n++ )
444 : {
445 0 : const double fByp = double( n ) * fTox * fBy - fBym;
446 0 : fBym = fBy;
447 0 : fBy = fByp;
448 : }
449 :
450 0 : return fBy;
451 : }
452 : }
453 : }
454 :
455 : } // namespace analysis
456 : } // namespace sca
457 :
458 : /* vim:set shiftwidth=4 softtabstop=4 expandtab: */
|
