?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Update ooo320-m1
[ooovba.git] / scaddins / source / analysis / bessel.cxx
blobe354695d98ab7492c14a3e3878645bfce3baef46
1 /*************************************************************************
2 *
3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
4 *
5 * Copyright 2008 by Sun Microsystems, Inc.
6 *
7 * OpenOffice.org - a multi-platform office productivity suite
8 *
9 * $RCSfile: bessel.cxx,v $
10 * $Revision: 1.6 $
11 *
12 * This file is part of OpenOffice.org.
13 *
14 * OpenOffice.org is free software: you can redistribute it and/or modify
15 * it under the terms of the GNU Lesser General Public License version 3
16 * only, as published by the Free Software Foundation.
17 *
18 * OpenOffice.org is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU Lesser General Public License version 3 for more details
22 * (a copy is included in the LICENSE file that accompanied this code).
23 *
24 * You should have received a copy of the GNU Lesser General Public License
25 * version 3 along with OpenOffice.org. If not, see
26 * <http://www.openoffice.org/license.html>
27 * for a copy of the LGPLv3 License.
28 *
29 ************************************************************************/
30
31 #include "bessel.hxx"
32 #include "analysishelper.hxx"
33
34 #include <rtl/math.hxx>
35
36 using ::com::sun::star::lang::IllegalArgumentException;
37
38 namespace sca {
39 namespace analysis {
40
41 // ============================================================================
42
43 const double f_PI = 3.1415926535897932385;
44 const double f_2_PI = 2.0 * f_PI;
45 const double f_PI_DIV_2 = f_PI / 2.0;
46 const double f_PI_DIV_4 = f_PI / 4.0;
47 const double f_2_DIV_PI = 2.0 / f_PI;
48
49 const double THRESHOLD = 30.0; // Threshold for usage of approximation formula.
50 const double MAXEPSILON = 1e-10; // Maximum epsilon for end of iteration.
51 const sal_Int32 MAXITER = 100; // Maximum number of iterations.
52
53
54 // ============================================================================
55 // BESSEL J
56 // ============================================================================
57
58 /* The BESSEL function, first kind, unmodified:
59
60 inf (-1)^k (x/2)^(n+2k)
61 J_n(x) = SUM TERM(n,k) with TERM(n,k) := ---------------------
62 k=0 k! (n+k)!
63
64 Approximation for the BESSEL function, first kind, unmodified, for great x:
65
66 J_n(x) ~ sqrt( 2 / (PI x) ) cos( x - n PI/2 - PI/4 ) for x>=0.
67 */
68
69 // ----------------------------------------------------------------------------
70
71 double BesselJ( double x, sal_Int32 n ) throw( IllegalArgumentException )
72 {
73 if( n < 0 )
74 throw IllegalArgumentException();
75
76 double fResult = 0.0;
77 if( fabs( x ) <= THRESHOLD )
78 {
79 /* Start the iteration without TERM(n,0), which is set here.
80
81 TERM(n,0) = (x/2)^n / n!
82 */
83 double fTerm = pow( x / 2.0, (double)n ) / Fak( n );
84 sal_Int32 nK = 1; // Start the iteration with k=1.
85 fResult = fTerm; // Start result with TERM(n,0).
86
87 const double fSqrX = x * x / -4.0;
88
89 do
90 {
91 /* Calculation of TERM(n,k) from TERM(n,k-1):
92
93 (-1)^k (x/2)^(n+2k)
94 TERM(n,k) = ---------------------
95 k! (n+k)!
96
97 (-1)(-1)^(k-1) (x/2)^2 (x/2)^(n+2(k-1))
98 = -----------------------------------------
99 k (k-1)! (n+k) (n+k-1)!
100
101 -(x/2)^2 (-1)^(k-1) (x/2)^(n+2(k-1))
102 = ---------- * -----------------------------
103 k(n+k) (k-1)! (n+k-1)!
104
105 -(x^2/4)
106 = ---------- TERM(n,k-1)
107 k(n+k)
108 */
109 fTerm *= fSqrX; // defined above as -(x^2/4)
110 fTerm /= (nK * (nK + n));
111 fResult += fTerm;
112 }
113 while( (fabs( fTerm ) > MAXEPSILON) && (++nK < MAXITER) );
114 }
115 else
116 {
117 /* Approximation for the BESSEL function, first kind, unmodified:
118
119 J_n(x) ~ sqrt( 2 / (PI x) ) cos( x - n PI/2 - PI/4 ) for x>=0.
120
121 The BESSEL function J_n with n IN {0,2,4,...} is axially symmetric at
122 x=0, means J_n(x) = J_n(-x). Therefore the approximation for x<0 is:
123
124 J_n(x) = J_n(|x|) for x<0 and n IN {0,2,4,...}.
125
126 The BESSEL function J_n with n IN {1,3,5,...} is point-symmetric at
127 x=0, means J_n(x) = -J_n(-x). Therefore the approximation for x<0 is:
128
129 J_n(x) = -J_n(|x|) for x<0 and n IN {1,3,5,...}.
130 */
131 double fXAbs = fabs( x );
132 fResult = sqrt( f_2_DIV_PI / fXAbs ) * cos( fXAbs - n * f_PI_DIV_2 - f_PI_DIV_4 );
133 if( (n & 1) && (x < 0.0) )
134 fResult = -fResult;
135 }
136 return fResult;
137 }
138
139
140 // ============================================================================
141 // BESSEL I
142 // ============================================================================
143
144 /* The BESSEL function, first kind, modified:
145
146 inf (x/2)^(n+2k)
147 I_n(x) = SUM TERM(n,k) with TERM(n,k) := --------------
148 k=0 k! (n+k)!
149
150 Approximation for the BESSEL function, first kind, modified, for great x:
151
152 I_n(x) ~ e^x / sqrt( 2 PI x ) for x>=0.
153 */
154
155 // ----------------------------------------------------------------------------
156
157 double BesselI( double x, sal_Int32 n ) throw( IllegalArgumentException )
158 {
159 if( n < 0 )
160 throw IllegalArgumentException();
161
162 double fResult = 0.0;
163 if( fabs( x ) <= THRESHOLD )
164 {
165 /* Start the iteration without TERM(n,0), which is set here.
166
167 TERM(n,0) = (x/2)^n / n!
168 */
169 double fTerm = pow( x / 2.0, (double)n ) / Fak( n );
170 sal_Int32 nK = 1; // Start the iteration with k=1.
171 fResult = fTerm; // Start result with TERM(n,0).
172
173 const double fSqrX = x * x / 4.0;
174
175 do
176 {
177 /* Calculation of TERM(n,k) from TERM(n,k-1):
178
179 (x/2)^(n+2k)
180 TERM(n,k) = --------------
181 k! (n+k)!
182
183 (x/2)^2 (x/2)^(n+2(k-1))
184 = --------------------------
185 k (k-1)! (n+k) (n+k-1)!
186
187 (x/2)^2 (x/2)^(n+2(k-1))
188 = --------- * ------------------
189 k(n+k) (k-1)! (n+k-1)!
190
191 x^2/4
192 = -------- TERM(n,k-1)
193 k(n+k)
194 */
195 fTerm *= fSqrX; // defined above as x^2/4
196 fTerm /= (nK * (nK + n));
197 fResult += fTerm;
198 }
199 while( (fabs( fTerm ) > MAXEPSILON) && (++nK < MAXITER) );
200 }
201 else
202 {
203 /* Approximation for the BESSEL function, first kind, modified:
204
205 I_n(x) ~ e^x / sqrt( 2 PI x ) for x>=0.
206
207 The BESSEL function I_n with n IN {0,2,4,...} is axially symmetric at
208 x=0, means I_n(x) = I_n(-x). Therefore the approximation for x<0 is:
209
210 I_n(x) = I_n(|x|) for x<0 and n IN {0,2,4,...}.
211
212 The BESSEL function I_n with n IN {1,3,5,...} is point-symmetric at
213 x=0, means I_n(x) = -I_n(-x). Therefore the approximation for x<0 is:
214
215 I_n(x) = -I_n(|x|) for x<0 and n IN {1,3,5,...}.
216 */
217 double fXAbs = fabs( x );
218 fResult = exp( fXAbs ) / sqrt( f_2_PI * fXAbs );
219 if( (n & 1) && (x < 0.0) )
220 fResult = -fResult;
221 }
222 return fResult;
223 }
224
225
226 // ============================================================================
227
228 double Besselk0( double fNum ) throw( IllegalArgumentException )
229 {
230 double fRet;
231
232 if( fNum <= 2.0 )
233 {
234 double fNum2 = fNum * 0.5;
235 double y = fNum2 * fNum2;
236
237 fRet = -log( fNum2 ) * BesselI( fNum, 0 ) +
238 ( -0.57721566 + y * ( 0.42278420 + y * ( 0.23069756 + y * ( 0.3488590e-1 +
239 y * ( 0.262698e-2 + y * ( 0.10750e-3 + y * 0.74e-5 ) ) ) ) ) );
240 }
241 else
242 {
243 double y = 2.0 / fNum;
244
245 fRet = exp( -fNum ) / sqrt( fNum ) * ( 1.25331414 + y * ( -0.7832358e-1 +
246 y * ( 0.2189568e-1 + y * ( -0.1062446e-1 + y * ( 0.587872e-2 +
247 y * ( -0.251540e-2 + y * 0.53208e-3 ) ) ) ) ) );
248 }
249
250 return fRet;
251 }
252
253
254 double Besselk1( double fNum ) throw( IllegalArgumentException )
255 {
256 double fRet;
257
258 if( fNum <= 2.0 )
259 {
260 double fNum2 = fNum * 0.5;
261 double y = fNum2 * fNum2;
262
263 fRet = log( fNum2 ) * BesselI( fNum, 1 ) +
264 ( 1.0 + y * ( 0.15443144 + y * ( -0.67278579 + y * ( -0.18156897 + y * ( -0.1919402e-1 +
265 y * ( -0.110404e-2 + y * ( -0.4686e-4 ) ) ) ) ) ) )
266 / fNum;
267 }
268 else
269 {
270 double y = 2.0 / fNum;
271
272 fRet = exp( -fNum ) / sqrt( fNum ) * ( 1.25331414 + y * ( 0.23498619 +
273 y * ( -0.3655620e-1 + y * ( 0.1504268e-1 + y * ( -0.780353e-2 +
274 y * ( 0.325614e-2 + y * ( -0.68245e-3 ) ) ) ) ) ) );
275 }
276
277 return fRet;
278 }
279
280
281 double BesselK( double fNum, sal_Int32 nOrder ) throw( IllegalArgumentException )
282 {
283 switch( nOrder )
284 {
285 case 0: return Besselk0( fNum );
286 case 1: return Besselk1( fNum );
287 default:
288 {
289 double fBkp;
290
291 double fTox = 2.0 / fNum;
292 double fBkm = Besselk0( fNum );
293 double fBk = Besselk1( fNum );
294
295 for( sal_Int32 n = 1 ; n < nOrder ; n++ )
296 {
297 fBkp = fBkm + double( n ) * fTox * fBk;
298 fBkm = fBk;
299 fBk = fBkp;
300 }
301
302 return fBk;
303 }
304 }
305 }
306
307
308 double Bessely0( double fNum ) throw( IllegalArgumentException )
309 {
310 double fRet;
311
312 if( fNum < 8.0 )
313 {
314 double y = fNum * fNum;
315
316 double f1 = -2957821389.0 + y * ( 7062834065.0 + y * ( -512359803.6 +
317 y * ( 10879881.29 + y * ( -86327.92757 + y * 228.4622733 ) ) ) );
318
319 double f2 = 40076544269.0 + y * ( 745249964.8 + y * ( 7189466.438 +
320 y * ( 47447.26470 + y * ( 226.1030244 + y ) ) ) );
321
322 fRet = f1 / f2 + 0.636619772 * BesselJ( fNum, 0 ) * log( fNum );
323 }
324 else
325 {
326 double z = 8.0 / fNum;
327 double y = z * z;
328 double xx = fNum - 0.785398164;
329
330 double f1 = 1.0 + y * ( -0.1098628627e-2 + y * ( 0.2734510407e-4 +
331 y * ( -0.2073370639e-5 + y * 0.2093887211e-6 ) ) );
332
333 double f2 = -0.1562499995e-1 + y * ( 0.1430488765e-3 +
334 y * ( -0.6911147651e-5 + y * ( 0.7621095161e-6 +
335 y * ( -0.934945152e-7 ) ) ) );
336
337 fRet = sqrt( 0.636619772 / fNum ) * ( sin( xx ) * f1 + z * cos( xx ) * f2 );
338 }
339
340 return fRet;
341 }
342
343
344 double Bessely1( double fNum ) throw( IllegalArgumentException )
345 {
346 double fRet;
347
348 if( fNum < 8.0 )
349 {
350 double y = fNum * fNum;
351
352 double f1 = fNum * ( -0.4900604943e13 + y * ( 0.1275274390e13 +
353 y * ( -0.5153438139e11 + y * ( 0.7349264551e9 +
354 y * ( -0.4237922726e7 + y * 0.8511937935e4 ) ) ) ) );
355
356 double f2 = 0.2499580570e14 + y * ( 0.4244419664e12 +
357 y * ( 0.3733650367e10 + y * ( 0.2245904002e8 +
358 y * ( 0.1020426050e6 + y * ( 0.3549632885e3 + y ) ) ) ) );
359
360 fRet = f1 / f2 + 0.636619772 * ( BesselJ( fNum, 1 ) * log( fNum ) - 1.0 / fNum );
361 }
362 else
363 {
364 #if 0
365 // #i12430# don't know the intention of this piece of code...
366 double z = 8.0 / fNum;
367 double y = z * z;
368 double xx = fNum - 2.356194491;
369
370 double f1 = 1.0 + y * ( 0.183105e-2 + y * ( -0.3516396496e-4 +
371 y * ( 0.2457520174e-5 + y * ( -0.240337019e6 ) ) ) );
372
373 double f2 = 0.04687499995 + y * ( -0.2002690873e-3 +
374 y * ( 0.8449199096e-5 + y * ( -0.88228987e-6 +
375 y * 0.105787412e-6 ) ) );
376
377 fRet = sqrt( 0.636619772 / fNum ) * ( sin( xx ) * f1 + z * cos( xx ) * f2 );
378 #endif
379 // #i12430# ...but this seems to work much better.
380 fRet = sqrt( 0.636619772 / fNum ) * sin( fNum - 2.356194491 );
381 }
382
383 return fRet;
384 }
385
386
387 double BesselY( double fNum, sal_Int32 nOrder ) throw( IllegalArgumentException )
388 {
389 switch( nOrder )
390 {
391 case 0: return Bessely0( fNum );
392 case 1: return Bessely1( fNum );
393 default:
394 {
395 double fByp;
396
397 double fTox = 2.0 / fNum;
398 double fBym = Bessely0( fNum );
399 double fBy = Bessely1( fNum );
400
401 for( sal_Int32 n = 1 ; n < nOrder ; n++ )
402 {
403 fByp = double( n ) * fTox * fBy - fBym;
404 fBym = fBy;
405 fBy = fByp;
406 }
407
408 return fBy;
409 }
410 }
411 }
412
413
414 // ============================================================================
415
416 } // namespace analysis
417 } // namespace sca
418
vba interoperability for OOo