Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Functions
nf_specialFunctions.h File Reference
#include <math.h>
#include <float.h>
#include <nf_utilities.h>

Go to the source code of this file.

Functions

double nf_polevl (double x, double coef[], int N)
 
double nf_p1evl (double x, double coef[], int N)
 
double nf_exponentialIntegral (int n, double x, nfu_status *status)
 
double nf_gammaFunction (double x, nfu_status *status)
 
double nf_logGammaFunction (double x, nfu_status *status)
 
double nf_incompleteGammaFunction (double a, double x, nfu_status *status)
 
double nf_incompleteGammaFunctionComplementary (double a, double x, nfu_status *status)
 
double nf_amc_log_factorial (int)
 
double nf_amc_factorial (int)
 
double nf_amc_wigner_3j (int, int, int, int, int, int)
 
double nf_amc_wigner_6j (int, int, int, int, int, int)
 
double nf_amc_wigner_9j (int, int, int, int, int, int, int, int, int)
 
double nf_amc_racah (int, int, int, int, int, int)
 
double nf_amc_clebsh_gordan (int, int, int, int, int)
 
double nf_amc_z_coefficient (int, int, int, int, int, int)
 
double nf_amc_zbar_coefficient (int, int, int, int, int, int)
 
double nf_amc_reduced_matrix_element (int, int, int, int, int, int, int)
 

Function Documentation

◆ nf_amc_clebsh_gordan()

double nf_amc_clebsh_gordan ( int  j1,
int  j2,
int  m1,
int  m2,
int  j3 
)

Definition at line 288 of file nf_angularMomentumCoupling.cc.

288 {
289/*
290* Clebsh-Gordan coefficient
291* = <j1,j2,m1,m2|j3,m1+m2>
292* = (-)^(j1-j2+m1+m2) * std::sqrt(2*j3+1) * / j1 j2 j3 \
293* \ m1 m2 -m1-m2 /
294*
295* Note: Last value m3 is preset to m1+m2. Any other value will evaluate to 0.0.
296*/
297
298 int m3, x1, x2, x3, y1, y2, y3;
299 double cg = 0.0;
300
301 if ( j1 < 0 || j2 < 0 || j3 < 0) return( 0.0 );
302 if ( j1 + j2 + j3 > 2 * MAX_FACTORIAL ) return( INFINITY );
303
304 m3 = m1 + m2;
305
306 if ( ( x1 = ( j1 + m1 ) / 2 + 1 ) <= 0 ) return( 0.0 );
307 if ( ( x2 = ( j2 + m2 ) / 2 + 1 ) <= 0 ) return( 0.0 );
308 if ( ( x3 = ( j3 - m3 ) / 2 + 1 ) <= 0 ) return( 0.0 );
309
310 if ( ( y1 = x1 - m1 ) <= 0 ) return( 0.0 );
311 if ( ( y2 = x2 - m2 ) <= 0 ) return( 0.0 );
312 if ( ( y3 = x3 + m3 ) <= 0 ) return( 0.0 );
313
314 if ( j3 == 0 ){
315 if ( j1 == j2 ) cg = ( 1.0 / std::sqrt( (double)j1 + 1.0 ) * ( ( y1 % 2 == 0 ) ? -1:1 ) );
316 }
317 else if ( (j1 == 0 || j2 == 0 ) ){
318 if ( ( j1 + j2 ) == j3 ) cg = 1.0;
319 }
320 else {
321 if( m3 == 0 && std::abs( m1 ) <= 1 ){
322 if( m1 == 0 ) cg = cg1( x1, x2, x3 );
323 else cg = cg2( x1 + y1 - y2, x3 - 1, x1 + x2 - 2, x1 - y2, j1, j2, j3, m2 );
324 }
325 else if ( m2 == 0 && std::abs( m1 ) <=1 ){
326 cg = cg2( x1 - y2 + y3, x2 - 1, x1 + x3 - 2, x3 - y1, j1, j3, j3, m1 );
327 }
328 else if ( m1 == 0 && std::abs( m3 ) <= 1 ){
329 cg = cg2( x1, x1 - 1, x2 + x3 - 2, x2 - y3, j2, j3, j3, -m3 );
330 }
331 else cg = cg3( x1, x2, x3, y1, y2, y3 );
332 }
333
334 return( cg );
335}

Referenced by nf_amc_reduced_matrix_element(), nf_amc_wigner_3j(), nf_amc_z_coefficient(), and nf_amc_zbar_coefficient().

◆ nf_amc_factorial()

double nf_amc_factorial ( int  n)

Definition at line 96 of file nf_angularMomentumCoupling.cc.

96 {
97/*
98* returns n! for pre-computed table. INFINITY is return if n is negative or too large.
99*/
100 return G4Exp( nf_amc_log_factorial( n ) );
101}
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:179
nf_amc_log_factorial
double nf_amc_log_factorial(int n)
Definition: nf_angularMomentumCoupling.cc:85

◆ nf_amc_log_factorial()

double nf_amc_log_factorial ( int  n)

Definition at line 85 of file nf_angularMomentumCoupling.cc.

85 {
86/*
87* returns ln( n! ).
88*/
89 if( n > MAX_FACTORIAL ) return( INFINITY );
90 if( n < 0 ) return( INFINITY );
91 return nf_amc_log_fact[n];
92}

Referenced by nf_amc_factorial().

◆ nf_amc_racah()

double nf_amc_racah ( int  j1,
int  j2,
int  l2,
int  l1,
int  j3,
int  l3 
)

Definition at line 253 of file nf_angularMomentumCoupling.cc.

253 {
254/*
255* Racah coefficient definition in Edmonds (AR Edmonds, "Angular Momentum in Quantum Mechanics", Princeton (1980) is
256* W(j1, j2, l2, l1 ; j3, l3) = (-1)^(j1+j2+l1+l2) * { j1 j2 j3 }
257* { l1 l2 l3 }
258* The call signature of W(...) appears jumbled, but hey, that's the convention.
259*
260* This convention is exactly that used by Blatt-Biedenharn (Rev. Mod. Phys. 24, 258 (1952)) too
261*/
262
263 double sig;
264
265 sig = ( ( ( j1 + j2 + l1 + l2 ) % 4 == 0 ) ? 1.0 : -1.0 );
266 return sig * nf_amc_wigner_6j( j1, j2, j3, l1, l2, l3 );
267}
nf_amc_wigner_6j
double nf_amc_wigner_6j(int j1, int j2, int j3, int j4, int j5, int j6)
Definition: nf_angularMomentumCoupling.cc:121

Referenced by nf_amc_wigner_9j(), nf_amc_z_coefficient(), and nf_amc_zbar_coefficient().

◆ nf_amc_reduced_matrix_element()

double nf_amc_reduced_matrix_element ( int  lt,
int  st,
int  jt,
int  l0,
int  j0,
int  l1,
int  j1 
)

Definition at line 473 of file nf_angularMomentumCoupling.cc.

473 {
474/*
475* Reduced Matrix Element for Tensor Operator
476* = < l1j1 || T(YL,sigma_S)J || l0j0 >
477*
478* M.B.Johnson, L.W.Owen, G.R.Satchler
479* Phys. Rev. 142, 748 (1966)
480* Note: definition differs from JOS by the factor sqrt(2j1+1)
481*/
482 int llt;
483 double x1, x2, x3, reduced_mat, clebsh_gordan;
484
485 if ( parity( lt ) != parity( l0 ) * parity( l1 ) ) return( 0.0 );
486 if ( std::abs( l0 - l1 ) > lt || ( l0 + l1 ) < lt ) return( 0.0 );
487 if ( std::abs( ( j0 - j1 ) / 2 ) > jt || ( ( j0 + j1 ) / 2 ) < jt ) return( 0.0 );
488
489 llt = 2 * lt;
490 jt *= 2;
491 st *= 2;
492
493 if( ( clebsh_gordan = nf_amc_clebsh_gordan( j1, j0, 1, -1, jt ) ) == INFINITY ) return( INFINITY );
494
495 reduced_mat = 1.0 / std::sqrt( 4 * M_PI ) * clebsh_gordan / std::sqrt( jt + 1.0 ) /* BRB jt + 1 <= 0? */
496 * std::sqrt( ( j0 + 1.0 ) * ( j1 + 1.0 ) * ( llt + 1.0 ) )
497 * parity( ( j1 - j0 ) / 2 ) * parity( ( -l0 + l1 + lt ) / 2 ) * parity( ( j0 - 1 ) / 2 );
498
499 if( st == 2 ){
500 x1 = ( l0 - j0 / 2.0 ) * ( j0 + 1.0 );
501 x2 = ( l1 - j1 / 2.0 ) * ( j1 + 1.0 );
502 if ( jt == llt ){
503 x3 = ( lt == 0 ) ? 0 : ( x1 - x2 ) / std::sqrt( lt * ( lt + 1.0 ) );
504 }
505 else if ( jt == ( llt - st ) ){
506 x3 = ( lt == 0 ) ? 0 : -( lt + x1 + x2 ) / std::sqrt( lt * ( 2.0 * lt + 1.0 ) );
507 }
508 else if ( jt == ( llt + st ) ){
509 x3 = ( lt + 1 - x1 - x2 ) / std::sqrt( ( 2.0 * lt + 1.0 ) * ( lt + 1.0 ) );
510 }
511 else{
512 x3 = 1.0;
513 }
514 }
515 else x3 = 1.0;
516 reduced_mat *= x3;
517
518 return( reduced_mat );
519}
M_PI
#define M_PI
Definition: SbMath.h:33
nf_amc_clebsh_gordan
double nf_amc_clebsh_gordan(int j1, int j2, int m1, int m2, int j3)
Definition: nf_angularMomentumCoupling.cc:288

◆ nf_amc_wigner_3j()

double nf_amc_wigner_3j ( int  j1,
int  j2,
int  j3,
int  j4,
int  j5,
int  j6 
)

Definition at line 105 of file nf_angularMomentumCoupling.cc.

105 {
106/*
107* Wigner's 3J symbol (similar to Clebsh-Gordan)
108* = / j1 j2 j3 \
109* \ j4 j5 j6 /
110*/
111 double cg;
112
113 if( ( j4 + j5 + j6 ) != 0 ) return( 0.0 );
114 if( ( cg = nf_amc_clebsh_gordan( j1, j2, j4, j5, j3 ) ) == 0.0 ) return ( 0.0 );
115 if( cg == INFINITY ) return( cg );
116 return( ( ( ( j1 - j2 - j6 ) % 4 == 0 ) ? 1.0 : -1.0 ) * cg / std::sqrt( j3 + 1.0 ) ); /* BRB j3 + 1 <= 0? */
117}

◆ nf_amc_wigner_6j()

double nf_amc_wigner_6j ( int  j1,
int  j2,
int  j3,
int  j4,
int  j5,
int  j6 
)

Definition at line 121 of file nf_angularMomentumCoupling.cc.

121 {
122/*
123* Wigner's 6J symbol (similar to Racah)
124* = { j1 j2 j3 }
125* { j4 j5 j6 }
126*/
127 int i, x[6];
128
129 x[0] = j1; x[1] = j2; x[2] = j3; x[3] = j4; x[4] = j5; x[5] = j6;
130 for( i = 0; i < 6; i++ ) if ( x[i] == 0 ) return( w6j0( i, x ) );
131
132 return( w6j1( x ) );
133}

Referenced by nf_amc_racah().

◆ nf_amc_wigner_9j()

double nf_amc_wigner_9j ( int  j1,
int  j2,
int  j3,
int  j4,
int  j5,
int  j6,
int  j7,
int  j8,
int  j9 
)

Definition at line 226 of file nf_angularMomentumCoupling.cc.

226 {
227/*
228* Wigner's 9J symbol
229* / j1 j2 j3 \
230* = | j4 j5 j6 |
231* \ j7 j8 j9 /
232*
233*/
234 int i, i0, i1;
235 double rac;
236
237 i0 = max3( std::abs( j1 - j9 ), std::abs( j2 - j6 ), std::abs( j4 - j8 ) );
238 i1 = min3( ( j1 + j9 ), ( j2 + j6 ), ( j4 + j8 ) );
239
240 rac = 0.0;
241 for ( i = i0; i <= i1; i += 2 ){
242 rac += nf_amc_racah( j1, j4, j9, j8, j7, i )
243 * nf_amc_racah( j2, j5, i, j4, j8, j6 )
244 * nf_amc_racah( j9, i, j3, j2, j1, j6 ) * ( i + 1 );
245 if( rac == INFINITY ) return( INFINITY );
246 }
247
248 return( ( ( (int)( ( j1 + j3 + j5 + j8 ) / 2 + j2 + j4 + j9 ) % 4 == 0 ) ? 1.0 : -1.0 ) * rac );
249}
nf_amc_racah
double nf_amc_racah(int j1, int j2, int l2, int l1, int j3, int l3)
Definition: nf_angularMomentumCoupling.cc:253

◆ nf_amc_z_coefficient()

double nf_amc_z_coefficient ( int  l1,
int  j1,
int  l2,
int  j2,
int  s,
int  ll 
)

Definition at line 437 of file nf_angularMomentumCoupling.cc.

437 {
438/*
439* Biedenharn's Z-coefficient coefficient
440* = Z(l1 j1 l2 j2 | S L )
441*/
442 double z, clebsh_gordan = nf_amc_clebsh_gordan( l1, l2, 0, 0, ll ), racah = nf_amc_racah( l1, j1, l2, j2, s, ll );
443
444 if( ( clebsh_gordan == INFINITY ) || ( racah == INFINITY ) ) return( INFINITY );
445 z = ( ( ( -l1 + l2 + ll ) % 8 == 0 ) ? 1.0 : -1.0 )
446 * std::sqrt( l1 + 1.0 ) * std::sqrt( l2 + 1.0 ) * std::sqrt( j1 + 1.0 ) * std::sqrt( j2 + 1.0 ) * clebsh_gordan * racah;
447
448 return( z );
449}

◆ nf_amc_zbar_coefficient()

double nf_amc_zbar_coefficient ( int  l1,
int  j1,
int  l2,
int  j2,
int  s,
int  ll 
)

Definition at line 453 of file nf_angularMomentumCoupling.cc.

453 {
454/*
455* Lane & Thomas's Zbar-coefficient coefficient
456* = Zbar(l1 j1 l2 j2 | S L )
457* = (-i)^( -l1 + l2 + ll ) * Z(l1 j1 l2 j2 | S L )
458*
459* Lane & Thomas Rev. Mod. Phys. 30, 257-353 (1958).
460* Note, Lane & Thomas define this because they did not like the different phase convention in Blatt & Biedenharn's Z coefficient. They changed it to get better time-reversal behavior.
461* Froehner uses Lane & Thomas convention as does T. Kawano.
462*/
463 double zbar, clebsh_gordan = nf_amc_clebsh_gordan( l1, l2, 0, 0, ll ), racah = nf_amc_racah( l1, j1, l2, j2, s, ll );
464
465 if( ( clebsh_gordan == INFINITY ) || ( racah == INFINITY ) ) return( INFINITY );
466 zbar = std::sqrt( l1 + 1.0 ) * std::sqrt( l2 + 1.0 ) * std::sqrt( j1 + 1.0 ) * std::sqrt( j2 + 1.0 ) * clebsh_gordan * racah;
467
468 return( zbar );
469}

◆ nf_exponentialIntegral()

double nf_exponentialIntegral ( int  n,
double  x,
nfu_status *  status 
)

Definition at line 28 of file nf_exponentialIntegral.cc.

28 {
29
30 int i, ii, nm1;
31 double a, b, c, d, del, fact, h, psi;
32 double ans = 0.0;
33
34 *status = nfu_badInput;
35 if( !isfinite( x ) ) return( x );
36 *status = nfu_Okay;
37
38 nm1 = n - 1;
39 if( ( n < 0 ) || ( x < 0.0 ) || ( ( x == 0.0 ) && ( ( n == 0 ) || ( n == 1 ) ) ) ) {
40 *status = nfu_badInput; }
41 else {
42 if( n == 0 ) {
43 ans = G4Exp( -x ) / x; } /* Special case */
44 else if( x == 0.0 ) {
45 ans = 1.0 / nm1; } /* Another special case */
46 else if( x > 1.0 ) { /* Lentz's algorithm */
47 b = x + n;
48 c = 1.0 / FPMIN;
49 d = 1.0 / b;
50 h = d;
51 for( i = 1; i <= MAXIT; i++ ) {
52 a = -i * ( nm1 + i );
53 b += 2.0;
54 d = 1.0 / ( a * d + b ); /* Denominators cannot be zero */
55 c = b + a / c;
56 del = c * d;
57 h *= del;
58 if( fabs( del - 1.0 ) < EPS ) return( h * G4Exp( -x ) );
59 }
60 *status = nfu_failedToConverge; }
61 else {
62 ans = ( nm1 != 0 ) ? 1.0 / nm1 : -G4Log(x) - EULER; /* Set first term */
63 fact = 1.0;
64 for( i = 1; i <= MAXIT; i++ ) {
65 fact *= -x / i;
66 if( i != nm1 ) {
67 del = -fact / ( i - nm1 ); }
68 else {
69 psi = -EULER; /* Compute psi(n) */
70 for( ii = 1; ii <= nm1; ii++ ) psi += 1.0 / ii;
71 del = fact * ( -G4Log( x ) + psi );
72 }
73 ans += del;
74 if( fabs( del ) < fabs( ans ) * EPS ) return( ans );
75 }
76 *status = nfu_failedToConverge;
77 }
78 }
79 return( ans );
80}
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:226
FPMIN
#define FPMIN
Definition: nf_exponentialIntegral.cc:22
EPS
#define EPS
Definition: nf_exponentialIntegral.cc:23
EULER
#define EULER
Definition: nf_exponentialIntegral.cc:20
MAXIT
#define MAXIT
Definition: nf_exponentialIntegral.cc:21
isfinite
#define isfinite
Definition: nf_stringToDoubles.cc:15
nfu_Okay
@ nfu_Okay
Definition: nf_utilities.h:25
nfu_failedToConverge
@ nfu_failedToConverge
Definition: nf_utilities.h:30
nfu_badInput
@ nfu_badInput
Definition: nf_utilities.h:29

◆ nf_gammaFunction()

double nf_gammaFunction ( double  x,
nfu_status *  status 
)

Definition at line 126 of file nf_gammaFunctions.cc.

126 {
127
128 double p, q, z;
129 int i, sgngam = 1;
130
131 *status = nfu_badInput;
132 if( !isfinite( x ) ) return( x );
133 *status = nfu_Okay;
134
135 q = fabs( x );
136
137 if( q > 33.0 ) {
138 if( x < 0.0 ) {
139 p = floor( q );
140 if( p == q ) goto goverf;
141 i = (int) p;
142 if( ( i & 1 ) == 0 ) sgngam = -1;
143 z = q - p;
144 if( z > 0.5 ) {
145 p += 1.0;
146 z = q - p;
147 }
148 z = q * sin( M_PI * z );
149 if( z == 0.0 ) goto goverf;
150 z = M_PI / ( fabs( z ) * stirf( q, status ) );
151 }
152 else {
153 z = stirf( x, status );
154 }
155 return( sgngam * z );
156 }
157
158 z = 1.0;
159 while( x >= 3.0 ) {
160 x -= 1.0;
161 z *= x;
162 } // Loop checking, 11.06.2015, T. Koi
163
164 while( x < 0.0 ) {
165 if( x > -1.E-9 ) goto small;
166 z /= x;
167 x += 1.0;
168 } // Loop checking, 11.06.2015, T. Koi
169
170 while( x < 2.0 ) {
171 if( x < 1.e-9 ) goto small;
172 z /= x;
173 x += 1.0;
174 } // Loop checking, 11.06.2015, T. Koi
175
176 if( x == 2.0 ) return( z );
177
178 x -= 2.0;
179 p = nf_polevl( x, P, 6 );
180 q = nf_polevl( x, Q, 7 );
181 return( z * p / q );
182
183small:
184 if( x == 0.0 ) goto goverf;
185 return( z / ( ( 1.0 + 0.5772156649015329 * x ) * x ) );
186
187goverf:
188 return( sgngam * DBL_MAX );
189}
nf_polevl
double nf_polevl(double x, double coef[], int N)
Definition: nf_polevl.cc:46
DBL_MAX
#define DBL_MAX
Definition: templates.hh:62

Referenced by nf_incompleteGammaFunction(), and nf_incompleteGammaFunctionComplementary().

◆ nf_incompleteGammaFunction()

double nf_incompleteGammaFunction ( double  a,
double  x,
nfu_status *  status 
)

Definition at line 155 of file nf_incompleteGammaFunctions.cc.

155 {
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}
nf_incompleteGammaFunctionComplementary
double nf_incompleteGammaFunctionComplementary(double a, double x, nfu_status *status)
Definition: nf_incompleteGammaFunctions.cc:88
nf_gammaFunction
double nf_gammaFunction(double x, nfu_status *status)
Definition: nf_gammaFunctions.cc:126
DBL_EPSILON
#define DBL_EPSILON
Definition: templates.hh:66

Referenced by nf_incompleteGammaFunctionComplementary().

◆ nf_incompleteGammaFunctionComplementary()

double nf_incompleteGammaFunctionComplementary ( double  a,
double  x,
nfu_status *  status 
)

Definition at line 88 of file nf_incompleteGammaFunctions.cc.

88 {
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}
nf_incompleteGammaFunction
double nf_incompleteGammaFunction(double a, double x, nfu_status *status)
Definition: nf_incompleteGammaFunctions.cc:155

Referenced by nf_incompleteGammaFunction().

◆ nf_logGammaFunction()

double nf_logGammaFunction ( double  x,
nfu_status *  status 
)

Definition at line 206 of file nf_gammaFunctions.cc.

206 {
207/* Logarithm of gamma function */
208
209 int sgngam;
210
211 *status = nfu_badInput;
212 if( !isfinite( x ) ) return( x );
213 *status = nfu_Okay;
214 return( lgam( x, &sgngam, status ) );
215}

◆ nf_p1evl()

double nf_p1evl ( double  x,
double  coef[],
int  N 
)

Definition at line 67 of file nf_polevl.cc.

67 {
68
69 double ans;
70 double *p;
71 int i;
72
73 p = coef;
74 ans = x + *p++;
75 i = N-1;
76
77 do {
78 ans = ans * x + *p++; }
79 while( --i ); // Loop checking, 11.06.2015, T. Koi
80
81 return( ans );
82}

◆ nf_polevl()

double nf_polevl ( double  x,
double  coef[],
int  N 
)

Definition at line 46 of file nf_polevl.cc.

46 {
47
48 double ans;
49 int i;
50 double *p;
51
52 p = coef;
53 ans = *p++;
54 i = N;
55
56 do {
57 ans = ans * x + *p++; }
58 while( --i ); // Loop checking, 11.06.2015, T. Koi
59
60 return( ans );
61}

Referenced by nf_gammaFunction().