Geant4 Cross Reference |
1 /* 1 /* 2 # <<BEGIN-copyright>> 2 # <<BEGIN-copyright>> 3 # <<END-copyright>> 3 # <<END-copyright>> 4 */ 4 */ 5 5 6 #include <float.h> 6 #include <float.h> 7 7 8 #include "nf_integration.h" 8 #include "nf_integration.h" 9 9 10 #if defined __cplusplus 10 #if defined __cplusplus 11 namespace GIDI { 11 namespace GIDI { 12 using namespace GIDI; 12 using namespace GIDI; 13 #endif 13 #endif 14 14 15 typedef struct nf_GnG_adaptiveQuadrature_info_ 15 typedef struct nf_GnG_adaptiveQuadrature_info_s { 16 nfu_status status; 16 nfu_status status; 17 nf_Legendre_GaussianQuadrature_callback in 17 nf_Legendre_GaussianQuadrature_callback integrandFunction; 18 void *argList; 18 void *argList; 19 nf_GnG_adaptiveQuadrature_callback quadrat 19 nf_GnG_adaptiveQuadrature_callback quadratureFunction; 20 double estimate; 20 double estimate; 21 int evaluations, maxDepth, maxDepthReached 21 int evaluations, maxDepth, maxDepthReached; 22 } nf_GnG_adaptiveQuadrature_info; 22 } nf_GnG_adaptiveQuadrature_info; 23 23 24 static double initialPoints[] = { 0.2311, 0.48 24 static double initialPoints[] = { 0.2311, 0.4860, 0.6068, 0.8913, 0.9501 }; 25 static int numberOfInitialPoints = sizeof( ini 25 static int numberOfInitialPoints = sizeof( initialPoints ) / sizeof( initialPoints[0] ); 26 26 27 static double nf_GnG_adaptiveQuadrature2( nf_G 27 static double nf_GnG_adaptiveQuadrature2( nf_GnG_adaptiveQuadrature_info *adaptiveQuadrature_info, double currentIntrgral, double x1, double x2, int depth ); 28 /* 28 /* 29 ============================================== 29 ============================================================ 30 */ 30 */ 31 nfu_status nf_GnG_adaptiveQuadrature( nf_GnG_a 31 nfu_status nf_GnG_adaptiveQuadrature( nf_GnG_adaptiveQuadrature_callback quadratureFunction, nf_Legendre_GaussianQuadrature_callback integrandFunction, 32 void *argList, double x1, double x2, int m 32 void *argList, double x1, double x2, int maxDepth, double tolerance, double *integral, long *evaluations ) { 33 /* 33 /* 34 * See W. Gander and W. Gautschi, "Adaptive q 34 * See W. Gander and W. Gautschi, "Adaptive quadrature--revisited", BIT 40 (2000), 84-101. 35 */ 35 */ 36 int i1; 36 int i1; 37 double estimate = 0., y1, integral_, coars 37 double estimate = 0., y1, integral_, coarse; 38 nfu_status status = nfu_Okay; 38 nfu_status status = nfu_Okay; 39 nf_GnG_adaptiveQuadrature_info adaptiveQua 39 nf_GnG_adaptiveQuadrature_info adaptiveQuadrature_info = { nfu_Okay, integrandFunction, argList, quadratureFunction, 0., 0, maxDepth, 0 }; 40 40 41 *integral = 0.; 41 *integral = 0.; 42 *evaluations = 0; 42 *evaluations = 0; 43 if( x1 == x2 ) return( nfu_Okay ); 43 if( x1 == x2 ) return( nfu_Okay ); 44 44 45 if( tolerance < 10 * DBL_EPSILON ) toleran 45 if( tolerance < 10 * DBL_EPSILON ) tolerance = 10 * DBL_EPSILON; 46 if( maxDepth > nf_GnG_adaptiveQuadrature_M 46 if( maxDepth > nf_GnG_adaptiveQuadrature_MaxMaxDepth ) maxDepth = nf_GnG_adaptiveQuadrature_MaxMaxDepth; 47 47 48 for( i1 = 0; i1 < numberOfInitialPoints; i 48 for( i1 = 0; i1 < numberOfInitialPoints; i1++ ) { 49 if( ( status = integrandFunction( x1 + 49 if( ( status = integrandFunction( x1 + ( x2 - x1 ) * initialPoints[i1], &y1, argList ) ) != nfu_Okay ) return( status ); 50 estimate += y1; 50 estimate += y1; 51 } 51 } 52 if( ( status = quadratureFunction( integra 52 if( ( status = quadratureFunction( integrandFunction, argList, x1, x2, &integral_ ) ) != nfu_Okay ) return( status ); 53 estimate = 0.5 * ( estimate * ( x2 - x1 ) 53 estimate = 0.5 * ( estimate * ( x2 - x1 ) / numberOfInitialPoints + integral_ ); 54 if( estimate == 0. ) estimate = x2 - x1; 54 if( estimate == 0. ) estimate = x2 - x1; 55 adaptiveQuadrature_info.estimate = toleran 55 adaptiveQuadrature_info.estimate = tolerance * estimate / DBL_EPSILON; 56 56 57 if( ( status = quadratureFunction( integra 57 if( ( status = quadratureFunction( integrandFunction, argList, x1, x2, &coarse ) ) != nfu_Okay ) return( status ); 58 integral_ = nf_GnG_adaptiveQuadrature2( &a 58 integral_ = nf_GnG_adaptiveQuadrature2( &adaptiveQuadrature_info, coarse, x1, x2, 0 ); 59 59 60 for( i1 = 0; i1 < 2; i1++ ) { /* Est 60 for( i1 = 0; i1 < 2; i1++ ) { /* Estimate may be off by more than a factor of 10. Iterate at most 2 times. */ 61 if( integral_ == 0. ) break; 61 if( integral_ == 0. ) break; 62 y1 = integral_ / estimate; 62 y1 = integral_ / estimate; 63 if( ( y1 > 0.1 ) && ( y1 < 10. ) ) bre 63 if( ( y1 > 0.1 ) && ( y1 < 10. ) ) break; 64 64 65 estimate = integral_; 65 estimate = integral_; 66 adaptiveQuadrature_info.estimate = tol 66 adaptiveQuadrature_info.estimate = tolerance * integral_ / DBL_EPSILON; 67 *evaluations += adaptiveQuadrature_inf 67 *evaluations += adaptiveQuadrature_info.evaluations; 68 adaptiveQuadrature_info.evaluations = 68 adaptiveQuadrature_info.evaluations = 0; 69 integral_ = nf_GnG_adaptiveQuadrature2 69 integral_ = nf_GnG_adaptiveQuadrature2( &adaptiveQuadrature_info, integral_, x1, x2, 0 ); 70 } 70 } 71 71 72 *evaluations += adaptiveQuadrature_info.ev 72 *evaluations += adaptiveQuadrature_info.evaluations; 73 if( adaptiveQuadrature_info.status == nfu_ 73 if( adaptiveQuadrature_info.status == nfu_Okay ) *integral = integral_; 74 return( adaptiveQuadrature_info.status ); 74 return( adaptiveQuadrature_info.status ); 75 } 75 } 76 /* 76 /* 77 ============================================== 77 ============================================================ 78 */ 78 */ 79 static double nf_GnG_adaptiveQuadrature2( nf_G 79 static double nf_GnG_adaptiveQuadrature2( nf_GnG_adaptiveQuadrature_info *adaptiveQuadrature_info, double coarse, double x1, double x2, int depth ) { 80 80 81 double xm, intregral1, intregral2, fine, e 81 double xm, intregral1, intregral2, fine, extrapolate; 82 82 83 if( adaptiveQuadrature_info->status != nfu 83 if( adaptiveQuadrature_info->status != nfu_Okay ) return( 0. ); 84 if( x1 == x2 ) return( 0. ); 84 if( x1 == x2 ) return( 0. ); 85 85 86 adaptiveQuadrature_info->evaluations++; 86 adaptiveQuadrature_info->evaluations++; 87 depth++; 87 depth++; 88 if( depth > adaptiveQuadrature_info->maxDe 88 if( depth > adaptiveQuadrature_info->maxDepthReached ) adaptiveQuadrature_info->maxDepthReached = depth; 89 89 90 xm = 0.5 * ( x1 + x2 ); 90 xm = 0.5 * ( x1 + x2 ); 91 if( ( adaptiveQuadrature_info->status = ad 91 if( ( adaptiveQuadrature_info->status = adaptiveQuadrature_info->quadratureFunction( adaptiveQuadrature_info->integrandFunction, 92 adaptiveQuadrature_info->argList, x1, 92 adaptiveQuadrature_info->argList, x1, xm, &intregral1 ) ) != nfu_Okay ) return( 0. ); 93 if( ( adaptiveQuadrature_info->status = ad 93 if( ( adaptiveQuadrature_info->status = adaptiveQuadrature_info->quadratureFunction( adaptiveQuadrature_info->integrandFunction, 94 adaptiveQuadrature_info->argList, xm, 94 adaptiveQuadrature_info->argList, xm, x2, &intregral2 ) ) != nfu_Okay ) return( 0. ); 95 fine = intregral1 + intregral2; 95 fine = intregral1 + intregral2; 96 extrapolate = ( 16. * fine - coarse ) / 15 96 extrapolate = ( 16. * fine - coarse ) / 15.; 97 if( extrapolate != 0 ) { 97 if( extrapolate != 0 ) { 98 if( adaptiveQuadrature_info->estimate 98 if( adaptiveQuadrature_info->estimate + ( extrapolate - fine ) == adaptiveQuadrature_info->estimate ) return( fine ); 99 } 99 } 100 if( depth > adaptiveQuadrature_info->maxDe 100 if( depth > adaptiveQuadrature_info->maxDepth ) return( fine ); 101 return( nf_GnG_adaptiveQuadrature2( adapti 101 return( nf_GnG_adaptiveQuadrature2( adaptiveQuadrature_info, intregral1, x1, xm, depth ) + 102 nf_GnG_adaptiveQuadrature2( adapti 102 nf_GnG_adaptiveQuadrature2( adaptiveQuadrature_info, intregral2, xm, x2, depth ) ); 103 } 103 } 104 104 105 #if defined __cplusplus 105 #if defined __cplusplus 106 } 106 } 107 #endif 107 #endif 108 108