| Geant4 Cross Reference |
1 /* 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_
16 nfu_status status;
17 nf_Legendre_GaussianQuadrature_callback in
18 void *argList;
19 nf_GnG_adaptiveQuadrature_callback quadrat
20 double estimate;
21 int evaluations, maxDepth, maxDepthReached
22 } nf_GnG_adaptiveQuadrature_info;
23
24 static double initialPoints[] = { 0.2311, 0.48
25 static int numberOfInitialPoints = sizeof( ini
26
27 static double nf_GnG_adaptiveQuadrature2( nf_G
28 /*
29 ==============================================
30 */
31 nfu_status nf_GnG_adaptiveQuadrature( nf_GnG_a
32 void *argList, double x1, double x2, int m
33 /*
34 * See W. Gander and W. Gautschi, "Adaptive q
35 */
36 int i1;
37 double estimate = 0., y1, integral_, coars
38 nfu_status status = nfu_Okay;
39 nf_GnG_adaptiveQuadrature_info adaptiveQua
40
41 *integral = 0.;
42 *evaluations = 0;
43 if( x1 == x2 ) return( nfu_Okay );
44
45 if( tolerance < 10 * DBL_EPSILON ) toleran
46 if( maxDepth > nf_GnG_adaptiveQuadrature_M
47
48 for( i1 = 0; i1 < numberOfInitialPoints; i
49 if( ( status = integrandFunction( x1 +
50 estimate += y1;
51 }
52 if( ( status = quadratureFunction( integra
53 estimate = 0.5 * ( estimate * ( x2 - x1 )
54 if( estimate == 0. ) estimate = x2 - x1;
55 adaptiveQuadrature_info.estimate = toleran
56
57 if( ( status = quadratureFunction( integra
58 integral_ = nf_GnG_adaptiveQuadrature2( &a
59
60 for( i1 = 0; i1 < 2; i1++ ) { /* Est
61 if( integral_ == 0. ) break;
62 y1 = integral_ / estimate;
63 if( ( y1 > 0.1 ) && ( y1 < 10. ) ) bre
64
65 estimate = integral_;
66 adaptiveQuadrature_info.estimate = tol
67 *evaluations += adaptiveQuadrature_inf
68 adaptiveQuadrature_info.evaluations =
69 integral_ = nf_GnG_adaptiveQuadrature2
70 }
71
72 *evaluations += adaptiveQuadrature_info.ev
73 if( adaptiveQuadrature_info.status == nfu_
74 return( adaptiveQuadrature_info.status );
75 }
76 /*
77 ==============================================
78 */
79 static double nf_GnG_adaptiveQuadrature2( nf_G
80
81 double xm, intregral1, intregral2, fine, e
82
83 if( adaptiveQuadrature_info->status != nfu
84 if( x1 == x2 ) return( 0. );
85
86 adaptiveQuadrature_info->evaluations++;
87 depth++;
88 if( depth > adaptiveQuadrature_info->maxDe
89
90 xm = 0.5 * ( x1 + x2 );
91 if( ( adaptiveQuadrature_info->status = ad
92 adaptiveQuadrature_info->argList, x1,
93 if( ( adaptiveQuadrature_info->status = ad
94 adaptiveQuadrature_info->argList, xm,
95 fine = intregral1 + intregral2;
96 extrapolate = ( 16. * fine - coarse ) / 15
97 if( extrapolate != 0 ) {
98 if( adaptiveQuadrature_info->estimate
99 }
100 if( depth > adaptiveQuadrature_info->maxDe
101 return( nf_GnG_adaptiveQuadrature2( adapti
102 nf_GnG_adaptiveQuadrature2( adapti
103 }
104
105 #if defined __cplusplus
106 }
107 #endif
108