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