xref: /AOO41X/main/scaddins/source/analysis/bessel.cxx (revision cdf0e10c4e3984b49a9502b011690b615761d4a3)

/*************************************************************************
 *
 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
 *
 * Copyright 2000, 2010 Oracle and/or its affiliates.
 *
 * OpenOffice.org - a multi-platform office productivity suite
 *
 * This file is part of OpenOffice.org.
 *
 * OpenOffice.org is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License version 3
 * only, as published by the Free Software Foundation.
 *
 * OpenOffice.org is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License version 3 for more details
 * (a copy is included in the LICENSE file that accompanied this code).
 *
 * You should have received a copy of the GNU Lesser General Public License
 * version 3 along with OpenOffice.org.  If not, see
 * <http://www.openoffice.org/license.html>
 * for a copy of the LGPLv3 License.
 *
 ************************************************************************/

#include "bessel.hxx"
#include "analysishelper.hxx"

#include <rtl/math.hxx>

using ::com::sun::star::lang::IllegalArgumentException;
using ::com::sun::star::sheet::NoConvergenceException;

namespace sca {
namespace analysis {

// ============================================================================

const double f_PI       = 3.1415926535897932385;
const double f_2_PI     = 2.0 * f_PI;
const double f_PI_DIV_2 = f_PI / 2.0;
const double f_PI_DIV_4 = f_PI / 4.0;
const double f_2_DIV_PI = 2.0 / f_PI;

const double THRESHOLD  = 30.0;     // Threshold for usage of approximation formula.
const double MAXEPSILON = 1e-10;    // Maximum epsilon for end of iteration.
const sal_Int32 MAXITER = 100;      // Maximum number of iterations.

// ============================================================================
// BESSEL J
// ============================================================================

/*  The BESSEL function, first kind, unmodified:
    The algorithm follows
    http://www.reference-global.com/isbn/978-3-11-020354-7
    Numerical Mathematics 1 / Numerische Mathematik 1,
    An algorithm-based introduction / Eine algorithmisch orientierte Einfuehrung
    Deuflhard, Peter; Hohmann, Andreas
    Berlin, New York (Walter de Gruyter) 2008
    4. ueberarb. u. erw. Aufl. 2008
    eBook ISBN: 978-3-11-020355-4
    Chapter 6.3.2 , algorithm 6.24
    The source is in German.
    The BesselJ-function is a special case of the adjoint summation with
    a_k = 2*(k-1)/x for k=1,...
    b_k = -1, for all k, directly substituted
    m_0=1, m_k=2 for k even, and m_k=0 for k odd, calculated on the fly
    alpha_k=1 for k=N and alpha_k=0 otherwise
*/

// ----------------------------------------------------------------------------

double BesselJ( double x, sal_Int32 N ) throw (IllegalArgumentException, NoConvergenceException)

{
    if( N < 0 )
        throw IllegalArgumentException();
    if (x==0.0)
        return (N==0) ? 1.0 : 0.0;

    /*  The algorithm works only for x>0, therefore remember sign. BesselJ
        with integer order N is an even function for even N (means J(-x)=J(x))
        and an odd function for odd N (means J(-x)=-J(x)).*/
    double fSign = (N % 2 == 1 && x < 0) ? -1.0 : 1.0;
    double fX = fabs(x);

    const double fMaxIteration = 9000000.0; //experimental, for to return in < 3 seconds
    double fEstimateIteration = fX * 1.5 + N;
    bool bAsymptoticPossible = pow(fX,0.4) > N;
    if (fEstimateIteration > fMaxIteration)
    {
        if (bAsymptoticPossible)
            return fSign * sqrt(f_2_DIV_PI/fX)* cos(fX-N*f_PI_DIV_2-f_PI_DIV_4);
        else
            throw NoConvergenceException();
    }

    double epsilon = 1.0e-15; // relative error
    bool bHasfound = false;
    double k= 0.0;
    // e_{-1} = 0; e_0 = alpha_0 / b_2
    double  u ; // u_0 = e_0/f_0 = alpha_0/m_0 = alpha_0

    // first used with k=1
    double m_bar;         // m_bar_k = m_k * f_bar_{k-1}
    double g_bar;         // g_bar_k = m_bar_k - a_{k+1} + g_{k-1}
    double g_bar_delta_u; // g_bar_delta_u_k = f_bar_{k-1} * alpha_k
                          // - g_{k-1} * delta_u_{k-1} - m_bar_k * u_{k-1}
    // f_{-1} = 0.0; f_0 = m_0 / b_2 = 1/(-1) = -1
    double g = 0.0;       // g_0= f_{-1} / f_0 = 0/(-1) = 0
    double delta_u = 0.0; // dummy initialize, first used with * 0
    double f_bar = -1.0;  // f_bar_k = 1/f_k, but only used for k=0

    if (N==0)
    {
        //k=0; alpha_0 = 1.0
        u = 1.0; // u_0 = alpha_0
        // k = 1.0; at least one step is necessary
        // m_bar_k = m_k * f_bar_{k-1} ==> m_bar_1 = 0.0
        g_bar_delta_u = 0.0;    // alpha_k = 0.0, m_bar = 0.0; g= 0.0
        g_bar = - 2.0/fX;       // k = 1.0, g = 0.0
        delta_u = g_bar_delta_u / g_bar;
        u = u + delta_u ;       // u_k = u_{k-1} + delta_u_k
        g = -1.0 / g_bar;       // g_k=b_{k+2}/g_bar_k
        f_bar = f_bar * g;      // f_bar_k = f_bar_{k-1}* g_k
        k = 2.0;
        // From now on all alpha_k = 0.0 and k > N+1
    }
    else
    {   // N >= 1 and alpha_k = 0.0 for k<N
        u=0.0; // u_0 = alpha_0
        for (k =1.0; k<= N-1; k = k + 1.0)
        {
            m_bar=2.0 * fmod(k-1.0, 2.0) * f_bar;
            g_bar_delta_u = - g * delta_u - m_bar * u; // alpha_k = 0.0
            g_bar = m_bar - 2.0*k/fX + g;
            delta_u = g_bar_delta_u / g_bar;
            u = u + delta_u;
            g = -1.0/g_bar;
            f_bar=f_bar * g;
        }
        // Step alpha_N = 1.0
        m_bar=2.0 * fmod(k-1.0, 2.0) * f_bar;
        g_bar_delta_u = f_bar - g * delta_u - m_bar * u; // alpha_k = 1.0
        g_bar = m_bar - 2.0*k/fX + g;
        delta_u = g_bar_delta_u / g_bar;
        u = u + delta_u;
        g = -1.0/g_bar;
        f_bar = f_bar * g;
        k = k + 1.0;
    }
    // Loop until desired accuracy, always alpha_k = 0.0
    do
    {
        m_bar = 2.0 * fmod(k-1.0, 2.0) * f_bar;
        g_bar_delta_u = - g * delta_u - m_bar * u;
        g_bar = m_bar - 2.0*k/fX + g;
        delta_u = g_bar_delta_u / g_bar;
        u = u + delta_u;
        g = -1.0/g_bar;
        f_bar = f_bar * g;
        bHasfound = (fabs(delta_u)<=fabs(u)*epsilon);
        k = k + 1.0;
    }
    while (!bHasfound && k <= fMaxIteration);
    if (bHasfound)
        return u * fSign;
    else
        throw NoConvergenceException(); // unlikely to happen
}

// ============================================================================
// BESSEL I
// ============================================================================

/*  The BESSEL function, first kind, modified:

                     inf                                  (x/2)^(n+2k)
        I_n(x)  =  SUM   TERM(n,k)   with   TERM(n,k) := --------------
                     k=0                                   k! (n+k)!

    No asymptotic approximation used, see issue 43040.
 */

// ----------------------------------------------------------------------------

double BesselI( double x, sal_Int32 n ) throw( IllegalArgumentException, NoConvergenceException )
{
    const double fEpsilon = 1.0E-15;
    const sal_Int32 nMaxIteration = 2000;
    const double fXHalf = x / 2.0;
    if( n < 0 )
        throw IllegalArgumentException();

    double fResult = 0.0;

    /*  Start the iteration without TERM(n,0), which is set here.

            TERM(n,0) = (x/2)^n / n!
     */
    sal_Int32 nK = 0;
    double fTerm = 1.0;
    // avoid overflow in Fak(n)
    for( nK = 1; nK <= n; ++nK )
    {
        fTerm = fTerm / static_cast< double >( nK ) * fXHalf;
    }
    fResult = fTerm;    // Start result with TERM(n,0).
    if( fTerm != 0.0 )
    {
        nK = 1;
        do
        {
            /*  Calculation of TERM(n,k) from TERM(n,k-1):

                                   (x/2)^(n+2k)
                    TERM(n,k)  =  --------------
                                    k! (n+k)!

                                   (x/2)^2 (x/2)^(n+2(k-1))
                               =  --------------------------
                                   k (k-1)! (n+k) (n+k-1)!

                                   (x/2)^2     (x/2)^(n+2(k-1))
                               =  --------- * ------------------
                                   k(n+k)      (k-1)! (n+k-1)!

                                   x^2/4
                               =  -------- TERM(n,k-1)
                                   k(n+k)
            */
        fTerm = fTerm * fXHalf / static_cast<double>(nK) * fXHalf / static_cast<double>(nK+n);
        fResult += fTerm;
        nK++;
        }
        while( (fabs( fTerm ) > fabs(fResult) * fEpsilon) && (nK < nMaxIteration) );

    }
    return fResult;
}


// ============================================================================

double Besselk0( double fNum ) throw( IllegalArgumentException, NoConvergenceException )
{
	double	fRet;

	if( fNum <= 2.0 )
	{
		double	fNum2 = fNum * 0.5;
		double	y = fNum2 * fNum2;

        fRet = -log( fNum2 ) * BesselI( fNum, 0 ) +
				( -0.57721566 + y * ( 0.42278420 + y * ( 0.23069756 + y * ( 0.3488590e-1 +
					y * ( 0.262698e-2 + y * ( 0.10750e-3 + y * 0.74e-5 ) ) ) ) ) );
	}
	else
	{
		double	y = 2.0 / fNum;

		fRet = exp( -fNum ) / sqrt( fNum ) * ( 1.25331414 + y * ( -0.7832358e-1 +
				y * ( 0.2189568e-1 + y * ( -0.1062446e-1 + y * ( 0.587872e-2 +
				y * ( -0.251540e-2 + y * 0.53208e-3 ) ) ) ) ) );
	}

	return fRet;
}


double Besselk1( double fNum ) throw( IllegalArgumentException, NoConvergenceException )
{
	double	fRet;

	if( fNum <= 2.0 )
	{
		double	fNum2 = fNum * 0.5;
		double	y = fNum2 * fNum2;

        fRet = log( fNum2 ) * BesselI( fNum, 1 ) +
				( 1.0 + y * ( 0.15443144 + y * ( -0.67278579 + y * ( -0.18156897 + y * ( -0.1919402e-1 +
					y * ( -0.110404e-2 + y * ( -0.4686e-4 ) ) ) ) ) ) )
				/ fNum;
	}
	else
	{
		double	y = 2.0 / fNum;

		fRet = exp( -fNum ) / sqrt( fNum ) * ( 1.25331414 + y * ( 0.23498619 +
				y * ( -0.3655620e-1 + y * ( 0.1504268e-1 + y * ( -0.780353e-2 +
				y * ( 0.325614e-2 + y * ( -0.68245e-3 ) ) ) ) ) ) );
	}

	return fRet;
}


double BesselK( double fNum, sal_Int32 nOrder ) throw( IllegalArgumentException, NoConvergenceException )
{
	switch( nOrder )
	{
		case 0:		return Besselk0( fNum );
		case 1:		return Besselk1( fNum );
		default:
		{
			double		fBkp;

			double		fTox = 2.0 / fNum;
			double		fBkm = Besselk0( fNum );
			double		fBk = Besselk1( fNum );

			for( sal_Int32 n = 1 ; n < nOrder ; n++ )
			{
				fBkp = fBkm + double( n ) * fTox * fBk;
				fBkm = fBk;
				fBk = fBkp;
			}

			return fBk;
		}
	}
}

// ============================================================================
// BESSEL Y
// ============================================================================

/*  The BESSEL function, second kind, unmodified:
    The algorithm for order 0 and for order 1 follows
    http://www.reference-global.com/isbn/978-3-11-020354-7
    Numerical Mathematics 1 / Numerische Mathematik 1,
    An algorithm-based introduction / Eine algorithmisch orientierte Einfuehrung
    Deuflhard, Peter; Hohmann, Andreas
    Berlin, New York (Walter de Gruyter) 2008
    4. ueberarb. u. erw. Aufl. 2008
    eBook ISBN: 978-3-11-020355-4
    Chapter 6.3.2 , algorithm 6.24
    The source is in German.
    See #i31656# for a commented version of the implementation, attachment #desc6
    http://www.openoffice.org/nonav/issues/showattachment.cgi/63609/Comments%20to%20the%20implementation%20of%20the%20Bessel%20functions.odt
*/

double Bessely0( double fX ) throw( IllegalArgumentException, NoConvergenceException )
{
    if (fX <= 0)
        throw IllegalArgumentException();
    const double fMaxIteration = 9000000.0; // should not be reached
    if (fX > 5.0e+6) // iteration is not considerable better then approximation
        return sqrt(1/f_PI/fX)
                *(rtl::math::sin(fX)-rtl::math::cos(fX));
    const double epsilon = 1.0e-15;
    const double EulerGamma = 0.57721566490153286060;
    double alpha = log(fX/2.0)+EulerGamma;
    double u = alpha;

    double k = 1.0;
    double m_bar = 0.0;
    double g_bar_delta_u = 0.0;
    double g_bar = -2.0 / fX;
    double delta_u = g_bar_delta_u / g_bar;
    double g = -1.0/g_bar;
    double f_bar = -1 * g;

    double sign_alpha = 1.0;
    double km1mod2;
    bool bHasFound = false;
    k = k + 1;
    do
    {
        km1mod2 = fmod(k-1.0,2.0);
        m_bar=(2.0*km1mod2) * f_bar;
        if (km1mod2 == 0.0)
            alpha = 0.0;
        else
        {
           alpha = sign_alpha * (4.0/k);
           sign_alpha = -sign_alpha;
        }
        g_bar_delta_u = f_bar * alpha - g * delta_u - m_bar * u;
        g_bar = m_bar - (2.0*k)/fX + g;
        delta_u = g_bar_delta_u / g_bar;
        u = u+delta_u;
        g = -1.0 / g_bar;
        f_bar = f_bar*g;
        bHasFound = (fabs(delta_u)<=fabs(u)*epsilon);
        k=k+1;
    }
    while (!bHasFound && k<fMaxIteration);
    if (bHasFound)
        return u*f_2_DIV_PI;
    else
        throw NoConvergenceException(); // not likely to happen
}

// See #i31656# for a commented version of this implementation, attachment #desc6
// http://www.openoffice.org/nonav/issues/showattachment.cgi/63609/Comments%20to%20the%20implementation%20of%20the%20Bessel%20functions.odt
double Bessely1( double fX ) throw( IllegalArgumentException, NoConvergenceException )
{
    if (fX <= 0)
        throw IllegalArgumentException();
    const double fMaxIteration = 9000000.0; // should not be reached
    if (fX > 5.0e+6) // iteration is not considerable better then approximation
        return - sqrt(1/f_PI/fX)
                *(rtl::math::sin(fX)+rtl::math::cos(fX));
    const double epsilon = 1.0e-15;
    const double EulerGamma = 0.57721566490153286060;
    double alpha = 1.0/fX;
    double f_bar = -1.0;
    double g = 0.0;
    double u = alpha;
    double k = 1.0;
    double m_bar = 0.0;
    alpha = 1.0 - EulerGamma - log(fX/2.0);
    double g_bar_delta_u = -alpha;
    double g_bar = -2.0 / fX;
    double delta_u = g_bar_delta_u / g_bar;
    u = u + delta_u;
    g = -1.0/g_bar;
    f_bar = f_bar * g;
    double sign_alpha = -1.0;
    double km1mod2; //will be (k-1) mod 2
    double q; // will be (k-1) div 2
    bool bHasFound = false;
    k = k + 1.0;
    do
    {
        km1mod2 = fmod(k-1.0,2.0);
        m_bar=(2.0*km1mod2) * f_bar;
        q = (k-1.0)/2.0;
        if (km1mod2 == 0.0) // k is odd
        {
           alpha = sign_alpha * (1.0/q + 1.0/(q+1.0));
           sign_alpha = -sign_alpha;
        }
        else
            alpha = 0.0;
        g_bar_delta_u = f_bar * alpha - g * delta_u - m_bar * u;
        g_bar = m_bar - (2.0*k)/fX + g;
        delta_u = g_bar_delta_u / g_bar;
        u = u+delta_u;
        g = -1.0 / g_bar;
        f_bar = f_bar*g;
        bHasFound = (fabs(delta_u)<=fabs(u)*epsilon);
        k=k+1;
    }
    while (!bHasFound && k<fMaxIteration);
    if (bHasFound)
        return -u*2.0/f_PI;
    else
        throw NoConvergenceException();
}

double BesselY( double fNum, sal_Int32 nOrder ) throw( IllegalArgumentException, NoConvergenceException )
{
    switch( nOrder )
    {
        case 0:     return Bessely0( fNum );
        case 1:     return Bessely1( fNum );
        default:
        {
            double      fByp;

            double      fTox = 2.0 / fNum;
            double      fBym = Bessely0( fNum );
            double      fBy = Bessely1( fNum );

            for( sal_Int32 n = 1 ; n < nOrder ; n++ )
            {
                fByp = double( n ) * fTox * fBy - fBym;
                fBym = fBy;
                fBy = fByp;
            }

            return fBy;
        }
    }
}

// ============================================================================

} // namespace analysis
} // namespace sca