EvtLi2Spence.cpp

Go to the documentation of this file.

/***********************************************************************
* Copyright 1998-2020 CERN for the benefit of the EvtGen authors       *
*                                                                      *
* This file is part of EvtGen.                                         *
*                                                                      *
* EvtGen is free software: you can redistribute it and/or modify       *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or    *
* (at your option) any later version.                                  *
*                                                                      *
* EvtGen 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 General Public License for more details.                         *
*                                                                      *
* You should have received a copy of the GNU General Public License    *
* along with EvtGen.  If not, see <https://www.gnu.org/licenses/>.     *
***********************************************************************/

#include <math.h>

const double a1 = -0.250000000000000;
const double a2 = -0.111111111111111;
const double a3 = -0.010000000000000;
const double a4 = -0.017006802721088;
const double a5 = -0.019444444444444;
const double a6 = -0.020661157024793;
const double a7 = -0.021417300648069;
const double a8 = -0.021948866377231;
const double a9 = -0.022349233811171;
const double a10 = -0.022663689135191;
const double zeta2 = 1.644934066848226;

double li2spence( double x )
{
    double result = 0;
    double xr = x;
    double xi = 0;
    double r2 = xr * xr + xi * xi;

    double y, z, z2, sp;

    if ( r2 > 1. && ( xr / r2 ) > 0.5 ) {
        y = ( x - 1.0 ) / x;
        z = -log( 1.0 - y );
        z2 = z * z;
        sp =
            z * ( 1.0 +
                  a1 * z *
                      ( 1.0 +
                        a2 * z *
                            ( 1.0 +
                              a3 * z2 *
                                  ( 1.0 +
                                    a4 * z2 *
                                        ( 1.0 +
                                          a5 * z2 *
                                              ( 1.0 +
                                                a6 * z2 *
                                                    ( 1.0 +
                                                      a7 * z2 *
                                                          ( 1.0 +
                                                            a8 * z2 *
                                                                ( 1.0 +
                                                                  a9 * z2 *
                                                                      ( 1.0 +
                                                                        a10 * z2 ) ) ) ) ) ) ) ) ) ) +
            zeta2 - log( x ) * log( 1.0 - x ) + 0.5 * log( x ) * log( x );
        result = sp;

    } else if ( r2 > 1 && ( xr / r2 ) <= 0.5 ) {
        y = 1.0 / x;
        z = -log( 1.0 - y );
        z2 = z * z;
        sp =
            -z *
                ( 1.0 +
                  a1 * z *
                      ( 1.0 +
                        a2 * z *
                            ( 1.0 +
                              a3 * z2 *
                                  ( 1.0 +
                                    a4 * z2 *
                                        ( 1.0 +
                                          a5 * z2 *
                                              ( 1.0 +
                                                a6 * z2 *
                                                    ( 1.0 +
                                                      a7 * z2 *
                                                          ( 1.0 +
                                                            a8 * z2 *
                                                                ( 1.0 +
                                                                  a9 * z2 *
                                                                      ( 1.0 +
                                                                        a10 * z2 ) ) ) ) ) ) ) ) ) ) +
            zeta2 - 0.5 * log( -x ) * log( -x );
        result = sp;

    } else if ( r2 == 1 && xi == 0 ) {
        sp = zeta2;
        result = sp;

    } else if ( r2 <= 1 && xr > 0.5 ) {
        y = 1.0 - x;
        z = -log( 1.0 - y );
        z2 = z * z;
        sp =
            -z *
                ( 1.0 +
                  a1 * z *
                      ( 1.0 +
                        a2 * z *
                            ( 1.0 +
                              a3 * z2 *
                                  ( 1.0 +
                                    a4 * z2 *
                                        ( 1.0 +
                                          a5 * z2 *
                                              ( 1.0 +
                                                a6 * z2 *
                                                    ( 1.0 +
                                                      a7 * z2 *
                                                          ( 1.0 +
                                                            a8 * z2 *
                                                                ( 1.0 +
                                                                  a9 * z2 *
                                                                      ( 1.0 +
                                                                        a10 * z2 ) ) ) ) ) ) ) ) ) ) +
            zeta2 - log( x ) * log( 1.0 - x );
        result = sp;

    } else if ( r2 <= 1 && xr <= 0.5 ) {
        y = x;
        z = -log( 1.0 - y );
        z2 = z * z;
        sp =
            z *
            ( 1.0 +
              a1 * z *
                  ( 1.0 +
                    a2 * z *
                        ( 1.0 +
                          a3 * z2 *
                              ( 1.0 +
                                a4 * z2 *
                                    ( 1.0 +
                                      a5 * z2 *
                                          ( 1.0 +
                                            a6 * z2 *
                                                ( 1.0 +
                                                  a7 * z2 *
                                                      ( 1.0 +
                                                        a8 * z2 *
                                                            ( 1.0 +
                                                              a9 * z2 *
                                                                  ( 1.0 +
                                                                    a10 * z2 ) ) ) ) ) ) ) ) ) );
        result = sp;
    }

    return result;
}