Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lend/src/nf_incompleteGammaFunctions.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 /*                            igam.c
  2  *
  3  *    Incomplete gamma integral
  4  *
  5  *
  6  * DESCRIPTION:
  7  *
  8  * The function is defined by
  9  *
 10  *                      x
 11  *                       -
 12  *                      | |  -t  a-1
 13  *  igam(a,x)  =        |   e   t   dt.
 14  *                    | |
 15  *                     -
 16  *                      0
 17  *
 18  *
 19  * In this implementation both arguments must be positive.
 20  * The integral is evaluated by either a power series or
 21  * continued fraction expansion, depending on the relative
 22  * values of a and x.
 23  *
 24  * ACCURACY:
 25  *
 26  *                      Relative error:
 27  * arithmetic   domain     # trials      peak         rms
 28  *    IEEE      0,30       200000       3.6e-14     2.9e-15
 29  *    IEEE      0,100      300000       9.9e-14     1.5e-14
 30  */
 31 /*                            igamc()
 32  *
 33  *    Complemented incomplete gamma integral
 34  *
 35  *
 36  * DESCRIPTION:
 37  *
 38  * The function is defined by
 39  *
 40  *
 41  *  igamc(a,x)   =   1 - igam(a,x)
 42  *
 43  *                       inf.
 44  *                         -
 45  *                        | |  -t  a-1
 46  *               =        |   e   t   dt.
 47  *                      | |
 48  *                       -
 49  *                        x
 50  *
 51  *
 52  * In this implementation both arguments must be positive.
 53  * The integral is evaluated by either a power series or
 54  * continued fraction expansion, depending on the relative
 55  * values of a and x.
 56  *
 57  * ACCURACY:
 58  *
 59  * Tested at random a, x.
 60  *                a         x                      Relative error:
 61  * arithmetic   domain   domain     # trials      peak         rms
 62  *    IEEE     0.5,100   0,100      200000       1.9e-14     1.7e-15
 63  *    IEEE     0.01,0.5  0,100      200000       1.4e-13     1.6e-15
 64  */
 65 
 66 /*
 67 Cephes Math Library Release 2.8:  June, 2000
 68 Copyright 1985, 1987, 2000 by Stephen L. Moshier
 69 */
 70 
 71 #include "nf_specialFunctions.h"
 72 
 73 #if defined __cplusplus
 74 #include <cmath>
 75 #include "G4Exp.hh"
 76 #include "G4Log.hh"
 77 namespace GIDI {
 78 using namespace GIDI;
 79 using namespace std;
 80 #endif
 81 
 82 static double big = 4.503599627370496e15;
 83 static double biginv =  2.22044604925031308085e-16;
 84 
 85 /*
 86 ************************************************************
 87 */
 88 double nf_incompleteGammaFunctionComplementary( double a, double x, nfu_status *status ) {
 89 
 90     double ans, ax, c, yc, r, t, y, z;
 91     double pk, pkm1, pkm2, qk, qkm1, qkm2;
 92 
 93     *status = nfu_badInput;
 94     if( !isfinite( x ) ) return( x );
 95     *status = nfu_Okay;
 96 
 97     if( ( x <= 0 ) || ( a <= 0 ) ) return( 1.0 );
 98     if( ( x < 1.0 ) || ( x < a ) ) return( nf_gammaFunction( a, status ) - nf_incompleteGammaFunction( a, x, status ) );
 99 
100     ax = G4Exp( a * G4Log( x ) - x );
101     if( ax == 0. ) return( 0.0 );
102 
103     if( x < 10000. ) {
104         y = 1.0 - a;                        /* continued fraction */
105         z = x + y + 1.0;
106         c = 0.0;
107         pkm2 = 1.0;
108         qkm2 = x;
109         pkm1 = x + 1.0;
110         qkm1 = z * x;
111         ans = pkm1 / qkm1;
112 
113         do {
114             c += 1.0;
115             y += 1.0;
116             z += 2.0;
117             yc = y * c;
118             pk = pkm1 * z  -  pkm2 * yc;
119             qk = qkm1 * z  -  qkm2 * yc;
120             if( qk != 0 ) {
121                 r = pk / qk;
122                 t = fabs( ( ans - r ) / r );
123                 ans = r; }
124             else {
125                 t = 1.0;
126             }
127             pkm2 = pkm1;
128             pkm1 = pk;
129             qkm2 = qkm1;
130             qkm1 = qk;
131             if( fabs( pk ) > big ) {
132                 pkm2 *= biginv;
133                 pkm1 *= biginv;
134                 qkm2 *= biginv;
135                 qkm1 *= biginv;
136             }
137         } while( t > DBL_EPSILON ); } // Loop checking, 11.06.2015, T. Koi
138     else {                                  /* Asymptotic expansion. */
139         y = 1. / x;
140         r = a;
141         c = 1.;
142         ans = 1.;
143         do {
144             a -= 1.;
145             c *= a * y;
146             ans += c;
147         } while( fabs( c ) > 100 * ans * DBL_EPSILON ); // Loop checking, 11.06.2015, T. Koi
148     }
149 
150     return( ans * ax );
151 }
152 /*
153 ************************************************************
154 */
155 double nf_incompleteGammaFunction( double a, double x, nfu_status *status ) {
156 /* left tail of incomplete gamma function:
157 *
158 *          inf.      k
159 *   a  -x   -       x
160 *  x  e     >   ----------
161 *           -     -
162 *          k=0   | (a+k+1)
163 */
164     double ans, ax, c, r;
165 
166     *status = nfu_badInput;
167     if( !isfinite( x ) ) return( x );
168     *status = nfu_Okay;
169 
170     if( ( x <= 0 ) || ( a <= 0 ) ) return( 0.0 );
171     if( ( x > 1.0 ) && ( x > a ) ) return( nf_gammaFunction( a, status ) - nf_incompleteGammaFunctionComplementary( a, x, status ) );
172 
173     ax = G4Exp( a * G4Log( x ) - x );       /* Compute  x**a * exp(-x) */
174     if( ax == 0. ) return( 0.0 );
175 
176     r = a;                                                      /* power series */
177     c = 1.0;
178     ans = 1.0;
179     do {
180         r += 1.0;
181         c *= x / r;
182         ans += c;
183     } while( c > ans * DBL_EPSILON ); // Loop checking, 11.06.2015, T. Koi
184 
185     return( ans * ax / a );
186 }
187 
188 #if defined __cplusplus
189 }
190 #endif
191