nf_gammaFunctions.cc File Reference

#include "nf_specialFunctions.h"

Go to the source code of this file.

Macros
#define	MAXGAM   171.624376956302725
#define	MAXLGM   2.556348e305
#define	MAXSTIR   143.01608

Functions
static double	lgam (double x, int *sgngam, nfu_status *status)
double	nf_gammaFunction (double x, nfu_status *status)
double	nf_logGammaFunction (double x, nfu_status *status)
static double	stirf (double x, nfu_status *status)

Variables
static double	A []
static double	B []
static double	C []
static double	LOGPI = 1.14472988584940017414
static double	LS2PI = 0.91893853320467274178
static double	P []
static double	Q []
static double	SQTPI = 2.50662827463100050242E0
static double	STIR [5] = { 7.873113957930936284e-4, -2.2954996161337812638e-4, -2.6813261780578123283e-3, 3.472222216054586673e-3, 8.3333333333348225713e-2 }

Macro Definition Documentation

MAXGAM

#define MAXGAM   171.624376956302725

Definition at line 93 of file nf_gammaFunctions.cc.

MAXLGM

#define MAXLGM   2.556348e305

Definition at line 201 of file nf_gammaFunctions.cc.

MAXSTIR

#define MAXSTIR   143.01608

Definition at line 99 of file nf_gammaFunctions.cc.

Function Documentation

lgam()

static double lgam ( double  x, int *  sgngam, nfu_status *  status  )

static

Definition at line 219 of file nf_gammaFunctions.cc.

                                                                {
 
    double p, q, u, w, z;
    int i;
 
    *sgngam = 1;
 
    if( x < -34.0 ) {
        q = -x;
        w = lgam( q, sgngam, status );                  /* note this modifies *sgngam! */
        p = floor( q );
        if( p == q ) goto lgsing;
        i = (int) p;
        if( ( i & 1 ) == 0 ) {
            *sgngam = -1; }
        else {
            *sgngam = 1;
        }
        z = q - p;
        if( z > 0.5 ) {
            p += 1.0;
            z = p - q;
        }
        z = q * sin( M_PI * z );
        if( z == 0.0 ) goto lgsing;
        z = LOGPI - G4Log( z ) - w;
        return( z );
    }
 
    if( x < 13.0 ) {
        z = 1.0;
        p = 0.0;
        u = x;
        while( u >= 3.0 ) {
            p -= 1.0;
            u = x + p;
            z *= u;
        } // Loop checking, 11.06.2015, T. Koi
        while( u < 2.0 ) {
            if( u == 0.0 ) goto lgsing;
            z /= u;
            p += 1.0;
            u = x + p;
        } // Loop checking, 11.06.2015, T. Koi
        if( z < 0.0 ) {
            *sgngam = -1;
            z = -z; }
        else {
            *sgngam = 1;
        }
        if( u == 2.0 ) return( G4Log( z ) );
        p -= 2.0;
        x = x + p;
        p = x * nf_polevl( x, B, 5 ) / nf_p1evl( x, C, 6);
        return( G4Log( z ) + p );
    }
 
    if( x > MAXLGM ) goto lgsing;
    q = ( x - 0.5 ) * G4Log( x ) - x + LS2PI;
    if( x > 1.0e8 ) return( q );
 
    p = 1.0 / ( x * x );
    if( x >= 1000.0 ) {
        q += ( ( 7.9365079365079365079365e-4 * p - 2.7777777777777777777778e-3 ) * p + 0.0833333333333333333333 ) / x; }
    else {
        q += nf_polevl( p, A, 4 ) / x;
    }
    return( q );
 
lgsing:
    return( *sgngam * DBL_MAX );
}

References A, B, C, DBL_MAX, G4Log(), lgam(), LOGPI, LS2PI, M_PI, MAXLGM, nf_p1evl(), and nf_polevl().

Referenced by lgam(), and nf_logGammaFunction().

nf_gammaFunction()

double nf_gammaFunction	(	double	x,
		nfu_status *	status
	)

Definition at line 126 of file nf_gammaFunctions.cc.

                                                        {
 
    double p, q, z;
    int i, sgngam = 1;
 
    *status = nfu_badInput;
    if( !isfinite( x ) ) return( x );
    *status = nfu_Okay;
 
    q = fabs( x );
 
    if( q > 33.0 ) {
        if( x < 0.0 ) {
            p = floor( q );
            if( p == q ) goto goverf;
            i = (int) p;
            if( ( i & 1 ) == 0 ) sgngam = -1;
            z = q - p;
            if( z > 0.5 ) {
                p += 1.0;
                z = q - p;
            }
            z = q * sin( M_PI * z );
            if( z == 0.0 ) goto goverf;
            z = M_PI / ( fabs( z ) * stirf( q, status ) );
        }
        else {
            z = stirf( x, status );
        }
        return( sgngam * z );
    }
 
    z = 1.0;
    while( x >= 3.0 ) {
        x -= 1.0;
        z *= x;
    } // Loop checking, 11.06.2015, T. Koi
 
    while( x < 0.0 ) {
        if( x > -1.E-9 ) goto small;
        z /= x;
        x += 1.0;
    } // Loop checking, 11.06.2015, T. Koi
 
    while( x < 2.0 ) {
        if( x < 1.e-9 ) goto small;
        z /= x;
        x += 1.0;
    } // Loop checking, 11.06.2015, T. Koi
 
    if( x == 2.0 ) return( z );
 
    x -= 2.0;
    p = nf_polevl( x, P, 6 );
    q = nf_polevl( x, Q, 7 );
    return( z * p / q );
 
small:
    if( x == 0.0 ) goto goverf;
    return( z / ( ( 1.0 + 0.5772156649015329 * x ) * x ) );
 
goverf:
    return( sgngam * DBL_MAX );
}

References DBL_MAX, isfinite, M_PI, nf_polevl(), nfu_badInput, nfu_Okay, P, Q, and stirf().

Referenced by nf_incompleteGammaFunction(), and nf_incompleteGammaFunctionComplementary().

nf_logGammaFunction()

double nf_logGammaFunction	(	double	x,
		nfu_status *	status
	)

Definition at line 206 of file nf_gammaFunctions.cc.

                                                           {
/* Logarithm of gamma function */
 
    int sgngam;
 
    *status = nfu_badInput;
    if( !isfinite( x ) ) return( x );
    *status = nfu_Okay;
    return( lgam( x, &sgngam, status ) );
}

References isfinite, lgam(), nfu_badInput, and nfu_Okay.

stirf()

static double stirf ( double  x, nfu_status *  status  )

static

Definition at line 106 of file nf_gammaFunctions.cc.

                                                {
/* Gamma function computed by Stirling's formula. The polynomial STIR is valid for 33 <= x <= 172.  */
 
    double y, w, v;
 
    w = 1.0 / x;
    w = 1.0 + w * nf_polevl( w, STIR, 4 );
    y = G4Exp( x );
    if( x > MAXSTIR ) {                     /* Avoid overflow in pow() */
        v = G4Pow::GetInstance()->powA( x, 0.5 * x - 0.25 );
        y = v * (v / y); }
    else {
        y = G4Pow::GetInstance()->powA( x, x - 0.5 ) / y;
    }
    y = SQTPI * y * w;
    return( y );
}

References G4Exp(), G4Pow::GetInstance(), MAXSTIR, nf_polevl(), G4Pow::powA(), SQTPI, and STIR.

Referenced by nf_gammaFunction().

Variable Documentation

A

double A[]

static

Initial value:

= { 8.11614167470508450300E-4, -5.95061904284301438324E-4, 7.93650340457716943945E-4,
                     -2.77777777730099687205E-3, 8.33333333333331927722E-2 }

Definition at line 194 of file nf_gammaFunctions.cc.

Referenced by lgam().

B

double B[]

static

Initial value:

= { -1.37825152569120859100E3, -3.88016315134637840924E4, -3.31612992738871184744E5,
                      -1.16237097492762307383E6, -1.72173700820839662146E6, -8.53555664245765465627E5 }

Definition at line 196 of file nf_gammaFunctions.cc.

Referenced by lgam().

C

double C[]

static

Initial value:

= { -3.51815701436523470549E2, -1.70642106651881159223E4, -2.20528590553854454839E5,
                      -1.13933444367982507207E6, -2.53252307177582951285E6, -2.01889141433532773231E6 }

Definition at line 198 of file nf_gammaFunctions.cc.

Referenced by lgam().

LOGPI

double LOGPI = 1.14472988584940017414

static

Definition at line 94 of file nf_gammaFunctions.cc.

Referenced by lgam().

LS2PI

double LS2PI = 0.91893853320467274178

static

Definition at line 200 of file nf_gammaFunctions.cc.

Referenced by lgam().

P

double P[]

static

Initial value:

= { 1.60119522476751861407E-4, 1.19135147006586384913E-3, 1.04213797561761569935E-2, 4.76367800457137231464E-2,
                      2.07448227648435975150E-1, 4.94214826801497100753E-1, 9.99999999999999996796E-1 }

Definition at line 89 of file nf_gammaFunctions.cc.

Referenced by G4VarNtp::addParticle(), G4LFission::ApplyYourself(), G4LEHadronProtonElastic::ApplyYourself(), G4LEnp::ApplyYourself(), G4LEpp::ApplyYourself(), G4DNAMolecularReactionData::ArrehniusParam(), G4QuasiElRatios::CalcQF2IN_Ratio(), G4ChipsAntiBaryonInelasticXS::CalculateCrossSection(), G4ChipsHyperonInelasticXS::CalculateCrossSection(), G4ChipsKaonMinusInelasticXS::CalculateCrossSection(), G4ChipsKaonPlusInelasticXS::CalculateCrossSection(), G4ChipsNeutronInelasticXS::CalculateCrossSection(), G4ChipsPionMinusInelasticXS::CalculateCrossSection(), G4ChipsPionPlusInelasticXS::CalculateCrossSection(), G4ChipsProtonInelasticXS::CalculateCrossSection(), G4GNASHTransitions::CalculateProbability(), G4PreCompoundTransitions::CalculateProbability(), G4QuasiElRatios::ChExer(), G4DiffractiveExcitation::ChooseP(), G4QGSDiffractiveExcitation::ChooseP(), G4GeomTools::ClosestPointOnSegment(), G4GeomTools::ClosestPointOnTriangle(), G4NistManager::ConstructNewIdealGasMaterial(), G4NistManager::ConstructNewMaterial(), G4ChipsAntiBaryonInelasticXS::CrossSectionFormula(), G4ChipsHyperonInelasticXS::CrossSectionFormula(), G4ChipsKaonMinusInelasticXS::CrossSectionFormula(), G4ChipsKaonPlusInelasticXS::CrossSectionFormula(), G4ChipsNeutronInelasticXS::CrossSectionFormula(), G4ChipsPionMinusInelasticXS::CrossSectionFormula(), G4ChipsPionPlusInelasticXS::CrossSectionFormula(), G4ChipsProtonInelasticXS::CrossSectionFormula(), G4ChipsAntiBaryonInelasticXS::CrossSectionLin(), G4ChipsHyperonInelasticXS::CrossSectionLin(), G4ChipsKaonMinusInelasticXS::CrossSectionLin(), G4ChipsKaonPlusInelasticXS::CrossSectionLin(), G4ChipsNeutronInelasticXS::CrossSectionLin(), G4ChipsPionMinusInelasticXS::CrossSectionLin(), G4ChipsPionPlusInelasticXS::CrossSectionLin(), G4ChipsProtonInelasticXS::CrossSectionLin(), G4ChipsAntiBaryonInelasticXS::CrossSectionLog(), G4ChipsHyperonInelasticXS::CrossSectionLog(), G4ChipsKaonMinusInelasticXS::CrossSectionLog(), G4ChipsKaonPlusInelasticXS::CrossSectionLog(), G4ChipsNeutronInelasticXS::CrossSectionLog(), G4ChipsPionMinusInelasticXS::CrossSectionLog(), G4ChipsPionPlusInelasticXS::CrossSectionLog(), G4ChipsProtonInelasticXS::CrossSectionLog(), G4MicroElecLOPhononModel::CrossSectionPerVolume(), G4FermiPhaseSpaceDecay::Decay(), G4GeomTools::DistancePointSegment(), G4PreCompoundAlpha::FactorialFactor(), G4PreCompoundDeuteron::FactorialFactor(), G4PreCompoundHe3::FactorialFactor(), G4PreCompoundTriton::FactorialFactor(), G4INCL::Particle::getBeta(), G4DELPHIMagField::GetFieldValue(), G4HadronNucleonXsc::HadronNucleonXscPDG(), G4HETCFragment::IntegrateEmissionProbability(), G4HETCAlpha::K(), G4HETCDeuteron::K(), G4HETCHe3::K(), G4HETCNeutron::K(), G4HETCProton::K(), G4HETCTriton::K(), G4INCL::NKbElasticChannel::KaonMomentum(), G4INCL::NKbToLpiChannel::KaonMomentum(), G4INCL::NKbToNKbChannel::KaonMomentum(), G4INCL::NKbToSpiChannel::KaonMomentum(), G4INCL::NpiToLKChannel::KaonMomentum(), G4INCL::NpiToSKChannel::KaonMomentum(), G4GDMLReadMaterials::MaterialRead(), nf_gammaFunction(), nf_Legendre_evauluateAtMu(), nf_Legendre_to_ptwXY2(), G4KL3DecayChannel::PhaseSpace(), G4GeomTools::PointInTriangle(), G4DNAMolecularReactionData::PolynomialParam(), G4PreCompoundIon::ProbabilityDistributionFunction(), G4PreCompoundNucleon::ProbabilityDistributionFunction(), G4FermiPhaseSpaceDecay::PtwoBody(), G4GDMLWriteMaterials::PWrite(), G4StatMFChannel::RotateMomentum(), G4QuasiElRatios::Scatter(), G4DNAMolecularReactionData::SetArrehniusParameterization(), G4ReactionTableMessenger::SetNewValue(), G4DNAMolecularReactionData::SetPolynomialParameterization(), and G4GDMLWriteStructure::TraverseVolumeTree().

Q

double Q[]

static

Initial value:

= { -2.31581873324120129819E-5, 5.39605580493303397842E-4, -4.45641913851797240494E-3, 1.18139785222060435552E-2,
                       3.58236398605498653373E-2, -2.34591795718243348568E-1, 7.14304917030273074085E-2, 1.00000000000000000320E0 }

Definition at line 91 of file nf_gammaFunctions.cc.

Referenced by G4LFission::ApplyYourself(), G4LEHadronProtonElastic::ApplyYourself(), G4LEnp::ApplyYourself(), G4LEpp::ApplyYourself(), G4ParticleHPCaptureFS::ApplyYourself(), G4ErrorSurfaceTrajState::BuildErrorMatrix(), G4Generator2BN::Calculatedsdkdt(), G4ECDecay::DecayIt(), G4EquilibriumEvaporator::deExcite(), G4NRESP71M03::DKINMA(), G4PenelopeRayleighModel::DumpFormFactorTable(), G4PenelopeRayleighModelMI::DumpFormFactorTable(), G4ionEffectiveCharge::EffectiveCharge(), G4ErrorFreeTrajState::G4ErrorFreeTrajState(), G4IntraNucleiCascader::generateCascade(), G4DELPHIMagField::GetFieldValue(), G4GaussXTRadiator::GetStackFactor(), G4TransparentRegXTRadiator::GetStackFactor(), G4XTRGammaRadModel::GetStackFactor(), G4XTRRegularRadModel::GetStackFactor(), G4XTRTransparentRegRadModel::GetStackFactor(), G4BetheHeitler5DModel::MaxDiffCrossSection(), MCGIDI_outputChannel_getFinalQ(), MCGIDI_outputChannel_parseFromTOM(), nf_gammaFunction(), G4PenelopeIonisationModel::SampleFinalStateElectron(), G4PenelopeIonisationModel::SampleFinalStatePositron(), G4AntiNuclElastic::SampleInvariantT(), and G4CollisionOutput::tuneSelectedPair().

SQTPI

double SQTPI = 2.50662827463100050242E0

static

Definition at line 95 of file nf_gammaFunctions.cc.

Referenced by stirf().

STIR

double STIR[5] = { 7.873113957930936284e-4, -2.2954996161337812638e-4, -2.6813261780578123283e-3, 3.472222216054586673e-3, 8.3333333333348225713e-2 }

static

Definition at line 98 of file nf_gammaFunctions.cc.

Referenced by stirf().