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