LCOV - code coverage report
Current view: top level - scaddins/source/analysis - bessel.cxx (source / functions) Hit Total Coverage
Test: commit e02a6cb2c3e2b23b203b422e4e0680877f232636 Lines: 0 213 0.0 %
Date: 2014-04-14 Functions: 0 8 0.0 %
Legend: Lines: hit not hit

          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             :             double      fBkp;
     288             : 
     289           0 :             double      fTox = 2.0 / fNum;
     290           0 :             double      fBkm = Besselk0( fNum );
     291           0 :             double      fBk = Besselk1( fNum );
     292             : 
     293           0 :             for( sal_Int32 n = 1 ; n < nOrder ; n++ )
     294             :             {
     295           0 :                 fBkp = fBkm + double( n ) * fTox * fBk;
     296           0 :                 fBkm = fBk;
     297           0 :                 fBk = fBkp;
     298             :             }
     299             : 
     300           0 :             return fBk;
     301             :         }
     302             :     }
     303             : }
     304             : 
     305             : 
     306             : // BESSEL Y
     307             : 
     308             : 
     309             : /*  The BESSEL function, second kind, unmodified:
     310             :     The algorithm for order 0 and for order 1 follows
     311             :     http://www.reference-global.com/isbn/978-3-11-020354-7
     312             :     Numerical Mathematics 1 / Numerische Mathematik 1,
     313             :     An algorithm-based introduction / Eine algorithmisch orientierte Einfuehrung
     314             :     Deuflhard, Peter; Hohmann, Andreas
     315             :     Berlin, New York (Walter de Gruyter) 2008
     316             :     4. ueberarb. u. erw. Aufl. 2008
     317             :     eBook ISBN: 978-3-11-020355-4
     318             :     Chapter 6.3.2 , algorithm 6.24
     319             :     The source is in German.
     320             :     See #i31656# for a commented version of the implementation, attachment #desc6
     321             :     http://www.openoffice.org/nonav/issues/showattachment.cgi/63609/Comments%20to%20the%20implementation%20of%20the%20Bessel%20functions.odt
     322             : */
     323             : 
     324           0 : double Bessely0( double fX ) throw( IllegalArgumentException, NoConvergenceException )
     325             : {
     326           0 :     if (fX <= 0)
     327           0 :         throw IllegalArgumentException();
     328           0 :     const double fMaxIteration = 9000000.0; // should not be reached
     329           0 :     if (fX > 5.0e+6) // iteration is not considerable better then approximation
     330           0 :         return sqrt(1/f_PI/fX)
     331           0 :                 *(rtl::math::sin(fX)-rtl::math::cos(fX));
     332           0 :     const double epsilon = 1.0e-15;
     333           0 :     const double EulerGamma = 0.57721566490153286060;
     334           0 :     double alpha = log(fX/2.0)+EulerGamma;
     335           0 :     double u = alpha;
     336             : 
     337           0 :     double k = 1.0;
     338           0 :     double m_bar = 0.0;
     339           0 :     double g_bar_delta_u = 0.0;
     340           0 :     double g_bar = -2.0 / fX;
     341           0 :     double delta_u = g_bar_delta_u / g_bar;
     342           0 :     double g = -1.0/g_bar;
     343           0 :     double f_bar = -1 * g;
     344             : 
     345           0 :     double sign_alpha = 1.0;
     346             :     double km1mod2;
     347           0 :     bool bHasFound = false;
     348           0 :     k = k + 1;
     349           0 :     do
     350             :     {
     351           0 :         km1mod2 = fmod(k-1.0,2.0);
     352           0 :         m_bar=(2.0*km1mod2) * f_bar;
     353           0 :         if (km1mod2 == 0.0)
     354           0 :             alpha = 0.0;
     355             :         else
     356             :         {
     357           0 :            alpha = sign_alpha * (4.0/k);
     358           0 :            sign_alpha = -sign_alpha;
     359             :         }
     360           0 :         g_bar_delta_u = f_bar * alpha - g * delta_u - m_bar * u;
     361           0 :         g_bar = m_bar - (2.0*k)/fX + g;
     362           0 :         delta_u = g_bar_delta_u / g_bar;
     363           0 :         u = u+delta_u;
     364           0 :         g = -1.0 / g_bar;
     365           0 :         f_bar = f_bar*g;
     366           0 :         bHasFound = (fabs(delta_u)<=fabs(u)*epsilon);
     367           0 :         k=k+1;
     368             :     }
     369           0 :     while (!bHasFound && k<fMaxIteration);
     370           0 :     if (bHasFound)
     371           0 :         return u*f_2_DIV_PI;
     372             :     else
     373           0 :         throw NoConvergenceException(); // not likely to happen
     374             : }
     375             : 
     376             : // See #i31656# for a commented version of this implementation, attachment #desc6
     377             : // http://www.openoffice.org/nonav/issues/showattachment.cgi/63609/Comments%20to%20the%20implementation%20of%20the%20Bessel%20functions.odt
     378           0 : double Bessely1( double fX ) throw( IllegalArgumentException, NoConvergenceException )
     379             : {
     380           0 :     if (fX <= 0)
     381           0 :         throw IllegalArgumentException();
     382           0 :     const double fMaxIteration = 9000000.0; // should not be reached
     383           0 :     if (fX > 5.0e+6) // iteration is not considerable better then approximation
     384           0 :         return - sqrt(1/f_PI/fX)
     385           0 :                 *(rtl::math::sin(fX)+rtl::math::cos(fX));
     386           0 :     const double epsilon = 1.0e-15;
     387           0 :     const double EulerGamma = 0.57721566490153286060;
     388           0 :     double alpha = 1.0/fX;
     389           0 :     double f_bar = -1.0;
     390           0 :     double u = alpha;
     391           0 :     double k = 1.0;
     392           0 :     double m_bar = 0.0;
     393           0 :     alpha = 1.0 - EulerGamma - log(fX/2.0);
     394           0 :     double g_bar_delta_u = -alpha;
     395           0 :     double g_bar = -2.0 / fX;
     396           0 :     double delta_u = g_bar_delta_u / g_bar;
     397           0 :     u = u + delta_u;
     398           0 :     double g = -1.0/g_bar;
     399           0 :     f_bar = f_bar * g;
     400           0 :     double sign_alpha = -1.0;
     401             :     double km1mod2; //will be (k-1) mod 2
     402             :     double q; // will be (k-1) div 2
     403           0 :     bool bHasFound = false;
     404           0 :     k = k + 1.0;
     405           0 :     do
     406             :     {
     407           0 :         km1mod2 = fmod(k-1.0,2.0);
     408           0 :         m_bar=(2.0*km1mod2) * f_bar;
     409           0 :         q = (k-1.0)/2.0;
     410           0 :         if (km1mod2 == 0.0) // k is odd
     411             :         {
     412           0 :            alpha = sign_alpha * (1.0/q + 1.0/(q+1.0));
     413           0 :            sign_alpha = -sign_alpha;
     414             :         }
     415             :         else
     416           0 :             alpha = 0.0;
     417           0 :         g_bar_delta_u = f_bar * alpha - g * delta_u - m_bar * u;
     418           0 :         g_bar = m_bar - (2.0*k)/fX + g;
     419           0 :         delta_u = g_bar_delta_u / g_bar;
     420           0 :         u = u+delta_u;
     421           0 :         g = -1.0 / g_bar;
     422           0 :         f_bar = f_bar*g;
     423           0 :         bHasFound = (fabs(delta_u)<=fabs(u)*epsilon);
     424           0 :         k=k+1;
     425             :     }
     426           0 :     while (!bHasFound && k<fMaxIteration);
     427           0 :     if (bHasFound)
     428           0 :         return -u*2.0/f_PI;
     429             :     else
     430           0 :         throw NoConvergenceException();
     431             : }
     432             : 
     433           0 : double BesselY( double fNum, sal_Int32 nOrder ) throw( IllegalArgumentException, NoConvergenceException )
     434             : {
     435           0 :     switch( nOrder )
     436             :     {
     437           0 :         case 0:     return Bessely0( fNum );
     438           0 :         case 1:     return Bessely1( fNum );
     439             :         default:
     440             :         {
     441             :             double      fByp;
     442             : 
     443           0 :             double      fTox = 2.0 / fNum;
     444           0 :             double      fBym = Bessely0( fNum );
     445           0 :             double      fBy = Bessely1( fNum );
     446             : 
     447           0 :             for( sal_Int32 n = 1 ; n < nOrder ; n++ )
     448             :             {
     449           0 :                 fByp = double( n ) * fTox * fBy - fBym;
     450           0 :                 fBym = fBy;
     451           0 :                 fBy = fByp;
     452             :             }
     453             : 
     454           0 :             return fBy;
     455             :         }
     456             :     }
     457             : }
     458             : 
     459             : } // namespace analysis
     460             : } // namespace sca
     461             : 
     462             : /* vim:set shiftwidth=4 softtabstop=4 expandtab: */

Generated by: LCOV version 1.10