Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
nf_GnG_adaptiveQuadrature.cc File Reference
#include <float.h>
#include "nf_integration.h"
Include dependency graph for nf_GnG_adaptiveQuadrature.cc:

Go to the source code of this file.

Classes

struct  nf_GnG_adaptiveQuadrature_info_s
 

Typedefs

typedef struct
nf_GnG_adaptiveQuadrature_info_s 
nf_GnG_adaptiveQuadrature_info
 

Functions

static double nf_GnG_adaptiveQuadrature2 (nf_GnG_adaptiveQuadrature_info *adaptiveQuadrature_info, double currentIntrgral, double x1, double x2, int depth)
 
nfu_status nf_GnG_adaptiveQuadrature (nf_GnG_adaptiveQuadrature_callback quadratureFunction, nf_Legendre_GaussianQuadrature_callback integrandFunction, void *argList, double x1, double x2, int maxDepth, double tolerance, double *integral, long *evaluations)
 

Variables

static double initialPoints [] = { 0.2311, 0.4860, 0.6068, 0.8913, 0.9501 }
 
static int numberOfInitialPoints = sizeof( initialPoints ) / sizeof( initialPoints[0] )
 

Typedef Documentation

Function Documentation

nfu_status nf_GnG_adaptiveQuadrature ( nf_GnG_adaptiveQuadrature_callback  quadratureFunction,
nf_Legendre_GaussianQuadrature_callback  integrandFunction,
void argList,
double  x1,
double  x2,
int  maxDepth,
double  tolerance,
double *  integral,
long *  evaluations 
)

Definition at line 31 of file nf_GnG_adaptiveQuadrature.cc.

32  {
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;
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 }
static double nf_GnG_adaptiveQuadrature2(nf_GnG_adaptiveQuadrature_info *adaptiveQuadrature_info, double currentIntrgral, double x1, double x2, int depth)
static int numberOfInitialPoints
#define nf_GnG_adaptiveQuadrature_MaxMaxDepth
static double initialPoints[]
enum nfu_status_e nfu_status
#define DBL_EPSILON
Definition: templates.hh:87

Here is the call graph for this function:

Here is the caller graph for this function:

static double nf_GnG_adaptiveQuadrature2 ( nf_GnG_adaptiveQuadrature_info adaptiveQuadrature_info,
double  currentIntrgral,
double  x1,
double  x2,
int  depth 
)
static

Definition at line 79 of file nf_GnG_adaptiveQuadrature.cc.

79  {
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 }
static double nf_GnG_adaptiveQuadrature2(nf_GnG_adaptiveQuadrature_info *adaptiveQuadrature_info, double currentIntrgral, double x1, double x2, int depth)
nf_Legendre_GaussianQuadrature_callback integrandFunction
nf_GnG_adaptiveQuadrature_callback quadratureFunction

Here is the caller graph for this function:

Variable Documentation

double initialPoints[] = { 0.2311, 0.4860, 0.6068, 0.8913, 0.9501 }
static

Definition at line 24 of file nf_GnG_adaptiveQuadrature.cc.

int numberOfInitialPoints = sizeof( initialPoints ) / sizeof( initialPoints[0] )
static

Definition at line 25 of file nf_GnG_adaptiveQuadrature.cc.