13 #if defined __cplusplus
25 double x1, x2, y1,
y2;
29 double x2,
double *integral );
40 switch( interpolation ) {
42 *value = 0.5 * ( y1 + y2 ) * ( x2 - x1 );
45 if( ( y1 <= 0. ) || ( y2 <= 0. ) ) {
49 if( std::fabs( r - 1. ) < 1e-4 ) {
51 *value = y1 * ( x2 - x1 ) / ( 1. + r * ( -0.5 + r * ( 1. / 3. + r * ( -0.25 + .2 * r ) ) ) ); }
53 *value = ( y2 - y1 ) * ( x2 - x1 ) /
G4Log( r );
58 if( ( x1 <= 0. ) || ( x2 <= 0. ) ) {
62 if( std::fabs( r - 1. ) < 1e-4 ) {
64 r = r * ( -0.5 + r * ( 1. / 3. + r * ( -0.25 + .2 * r ) ) );
65 *value = x1 * ( y2 - y1 ) * r / ( 1. + r ) + y2 * ( x2 - x1 ); }
67 *value = ( y1 - y2 ) * ( x2 - x1 ) /
G4Log( r ) + x2 * y2 - x1 * y1;
72 if( ( x1 <= 0. ) || ( x2 <= 0. ) || ( y1 <= 0. ) || ( y2 <= 0. ) ) {
76 double a,
z, lx, ly,
s, f;
79 if( std::fabs( r - 1. ) < 1e-4 ) {
80 ly = ( y2 - y1 ) / y1;
81 ly = ly * ( 1. + ly * ( -0.5 + ly * ( 1. / 3. - 0.25 * ly ) ) ); }
86 if( std::fabs( r - 1. ) < 1e-4 ) {
87 lx = ( x2 - x1 ) / x1;
88 lx = lx * ( 1 + lx * ( -0.5 + lx * ( 1. / 3. - 0.25 * lx ) ) ); }
93 if( std::fabs( r - 1. ) < 1e-3 ) {
100 for( i = 0, s = 0.; i <
n; i++, a++, f-- ) s = ( 1. + s ) * a * z / f;
101 *value = y1 * ( x2 - x1 ) * ( 1. + s ); }
108 *value = y1 * ( x2 - x1 );
121 double sum = 0., dSum,
x, y, x1, x2, y1, y2, _sign = 1.;
134 if( n < 2 )
return( 0. );
137 for( i = 0, point = ptwXY->
points; i < n; i++, point++ ) {
138 if( point->
x >= xMin )
break;
140 if( i == n )
return( 0. );
161 for( ; i <
n; i++, point++ ) {
176 return( _sign * sum );
199 if( status !=
nfu_Okay )
return( status );
203 for( i = 0; i < ptwXY->
length; i++ ) ptwXY->
points[i].
y /= sum;
213 if( ptwXY->
length < 2 )
return( 0. );
222 double sum = 0.,
x, y, x1, x2, y1, y2, _sign = 1.;
230 if( n < 2 )
return( 0. );
239 for( i = 0, point = ptwXY->
points; i < n; ++i, ++point ) {
240 if( point->
x >= xMin )
break;
242 if( i == n )
return( 0. );
256 for( ; i <
n; ++i, ++point ) {
268 sum += ( x2 - x1 ) * y1 * 3 * ( x1 + x2 );
271 sum += ( x2 - x1 ) * ( y1 * ( 2 * x1 + x2 ) + y2 * ( x1 + 2 * x2 ) );
276 if( x2 == xMax )
break;
279 return( _sign * sum / 6 );
287 if( ptwXY->
length < 2 )
return( 0. );
296 double sum = 0.,
x, y, x1, x2, y1, y2, _sign = 1., sqrt_x1, sqrt_x2, inv_apb,
c;
304 if( n < 2 )
return( 0. );
313 for( i = 0, point = ptwXY->
points; i < n; ++i, ++point ) {
314 if( point->
x >= xMin )
break;
316 if( i == n )
return( 0. );
330 sqrt_x2 = std::sqrt( x2 );
331 for( ; i <
n; ++i, ++point ) {
342 sqrt_x2 = std::sqrt( x2 );
343 inv_apb = sqrt_x1 + sqrt_x2;
344 c = 2. * ( sqrt_x1 * sqrt_x2 + x1 + x2 );
347 sum += ( sqrt_x2 - sqrt_x1 ) * y1 * 2.5 * c;
350 sum += ( sqrt_x2 - sqrt_x1 ) * ( y1 * ( c + x1 * ( 1. + sqrt_x2 / inv_apb ) ) + y2 * ( c + x2 * ( 1. + sqrt_x1 / inv_apb ) ) );
355 if( x2 == xMax )
break;
358 return( 2. / 15. * _sign * sum );
366 double x1, y1, x2, y2, y2p, xg1, xg2, sum;
371 if( ( *status = groupBoundaries->
status ) !=
nfu_Okay )
return( NULL );
377 if( ptwX_norm == NULL ) {
381 *status = ptwX_norm->
status;
392 if( ( groupedData =
ptwX_new( ngs, status ) ) == NULL )
goto err;
393 xg1 = groupBoundaries->
points[0];
396 for( igs = 0, i = 1; igs < ngs; igs++ ) {
397 xg2 = groupBoundaries->
points[igs+1];
400 for( ; i < f->
length; i++, x1 = x2, y1 = y2 ) {
402 if( x2 > xg2 )
break;
405 sum += ( y1 + y2p ) * ( x2 - x1 );
410 sum /= ( xg2 - xg1 ); }
412 if( ptwX_norm->
points[igs] == 0. ) {
416 sum /= ptwX_norm->
points[igs];
419 groupedData->
points[igs] = 0.5 * sum;
425 return( groupedData );
429 if( groupedData != NULL )
ptwX_free( groupedData );
439 double x1, fy1, gy1, x2, fy2, gy2, fy2p, gy2p, xg1, xg2, sum;
445 if( ( *status = groupBoundaries->
status ) !=
nfu_Okay )
return( NULL );
452 if( ptwX_norm == NULL ) {
465 if( ( ff->length == 0 ) || ( gg->length == 0 ) ) {
475 if( ( groupedData =
ptwX_new( ngs, status ) ) == NULL )
goto err;
476 xg1 = groupBoundaries->
points[0];
479 gy1 = g->points[0].y;
480 for( igs = 0, i = 1; igs < ngs; igs++ ) {
481 xg2 = groupBoundaries->
points[igs+1];
484 for( ; i < f->
length; i++, x1 = x2, fy1 = fy2, gy1 = gy2 ) {
486 if( x2 > xg2 )
break;
489 gy2p = gy2 = g->points[i].y;
491 sum += ( ( fy1 + fy2p ) * ( gy1 + gy2p ) + fy1 * gy1 + fy2p * gy2p ) * ( x2 - x1 );
496 sum /= ( xg2 - xg1 ); }
498 if( ptwX_norm->
points[igs] == 0. ) {
502 sum /= ptwX_norm->
points[igs];
505 groupedData->
points[igs] = sum / 6.;
514 return( groupedData );
522 if( groupedData != NULL )
ptwX_free( groupedData );
532 double x1, fy1, gy1, hy1, x2, fy2, gy2, hy2, fy2p, gy2p, hy2p, xg1, xg2, sum;
533 ptwXYPoints *f = NULL, *ff, *fff = NULL, *
g = NULL, *gg = NULL, *h = NULL, *hh = NULL;
539 if( ( *status = groupBoundaries->
status ) !=
nfu_Okay )
return( NULL );
547 if( ptwX_norm == NULL ) {
561 if( ( ff->length == 0 ) || ( gg->length == 0 ) || ( hh->length == 0 ) )
return(
ptwX_createLine( ngs, ngs, 0, 0, status ) );
571 if( ( groupedData =
ptwX_new( ngs, status ) ) == NULL )
goto err;
572 xg1 = groupBoundaries->
points[0];
575 gy1 = g->points[0].y;
576 hy1 = h->points[0].y;
577 for( igs = 0, i = 1; igs < ngs; igs++ ) {
578 xg2 = groupBoundaries->
points[igs+1];
581 for( ; i < f->
length; i++, x1 = x2, fy1 = fy2, gy1 = gy2, hy1 = hy2 ) {
583 if( x2 > xg2 )
break;
586 gy2p = gy2 = g->points[i].y;
588 hy2p = hy2 = h->points[i].y;
590 sum += ( ( fy1 + fy2p ) * ( gy1 + gy2p ) * ( hy1 + hy2p ) + 2 * fy1 * gy1 * hy1 + 2 * fy2p * gy2p * hy2p ) * ( x2 - x1 );
595 sum /= ( xg2 - xg1 ); }
597 if( ptwX_norm->
points[igs] == 0. ) {
601 sum /= ptwX_norm->
points[igs];
604 groupedData->
points[igs] = sum / 12.;
616 return( groupedData );
626 if( groupedData != NULL )
ptwX_free( groupedData );
636 double integral = 0., sum;
639 if( ( runningIntegral =
ptwX_new( ptwXY->
length, status ) ) == NULL )
goto err;
642 for( i = 1; i < ptwXY->
length; i++ ) {
648 return( runningIntegral );
651 if( runningIntegral != NULL )
ptwX_free( runningIntegral );
658 double xMin,
double xMax,
int degree,
int recursionLimit,
double tolerance,
nfu_status *status ) {
660 int64_t i1, i2, n1 = ptwXY->
length;
662 double integral = 0., integral_,
sign = -1., xa, xb;
668 if( xMin == xMax )
return( 0. );
669 if( n1 < 2 )
return( 0. );
679 if( xMin >= ptwXY->
points[n1-1].
x )
return( 0. );
680 if( xMax <= ptwXY->points[0].
x )
return( 0. );
682 for( i1 = 0; i1 < ( n1 - 1 ); i1++ ) {
683 if( ptwXY->
points[i1+1].
x > xMin )
break;
685 for( i2 = n1 - 1; i2 > i1; i2-- ) {
686 if( ptwXY->
points[i2-1].
x < xMax )
break;
688 point = &(ptwXY->
points[i1]);
691 integrateWithFunctionInfo.
func = func;
692 integrateWithFunctionInfo.
argList = argList;
694 integrateWithFunctionInfo.
x2 = point->
x;
695 integrateWithFunctionInfo.
y2 = point->
y;
698 for( ; i1 < i2; i1++ ) {
699 integrateWithFunctionInfo.
x1 = integrateWithFunctionInfo.
x2;
700 integrateWithFunctionInfo.
y1 = integrateWithFunctionInfo.
y2;
702 integrateWithFunctionInfo.
x2 = point->
x;
703 integrateWithFunctionInfo.
y2 = point->
y;
705 if( xb > xMax ) xb = xMax;
707 xa, xb, recursionLimit, tolerance, &integral_, &evaluations );
708 if( *status !=
nfu_Okay )
return( 0. );
709 integral += integral_;
719 double x2,
double *integral ) {
737 integrateWithFunctionInfo->
x1, integrateWithFunctionInfo->
y1,
738 integrateWithFunctionInfo->
x2, integrateWithFunctionInfo->
y2 ) ) ==
nfu_Okay ) {
739 status = integrateWithFunctionInfo->
func( x, y, integrateWithFunctionInfo->
argList );
745 #if defined __cplusplus
static G4Pow * GetInstance()
G4double powA(G4double A, G4double y) const
ptwXYPoints * ptwXY_intersectionWith_ptwX(ptwXYPoints *ptwXY, ptwXPoints *ptwX, nfu_status *status)
ptwXY_interpolation interpolation
double ptwXY_integrateWithWeight_sqrt_x(ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status)
ptwXPoints * ptwX_createLine(int64_t size, int64_t length, double slope, double offset, nfu_status *status)
nfu_status ptwXY_interpolatePoint(ptwXY_interpolation interpolation, double x, double *y, double x1, double y1, double x2, double y2)
std::vector< ExP01TrackerHit * > a
ptwXPoints * ptwXY_groupThreeFunctions(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, ptwXYPoints *ptwXY3, ptwXPoints *groupBoundaries, ptwXY_group_normType normType, ptwXPoints *ptwX_norm, nfu_status *status)
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
double ptwXY_integrateDomainWithWeight_x(ptwXYPoints *ptwXY, nfu_status *status)
double ptwXY_integrateDomain(ptwXYPoints *ptwXY, nfu_status *status)
nfu_status(* ptwXY_createFromFunction_callback)(double x, double *y, void *argList)
struct ptwXY_integrateWithFunctionInfo_s ptwXY_integrateWithFunctionInfo
ptwXPoints * ptwX_new(int64_t size, nfu_status *status)
nfu_status ptwXY_f_integrate(ptwXY_interpolation interpolation, double x1, double y1, double x2, double y2, double *value)
ptwXYPoints * ptwXY_union(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int unionOptions)
const XML_Char int const XML_Char * value
static constexpr double degree
ptwXPoints * ptwX_free(ptwXPoints *ptwX)
enum nfu_status_e nfu_status
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
double ptwXY_integrateWithFunction(ptwXYPoints *ptwXY, ptwXY_createFromFunction_callback func, void *argList, double xMin, double xMax, int degree, int recursionLimit, double tolerance, nfu_status *status)
ptwXPoints * ptwXY_runningIntegral(ptwXYPoints *ptwXY, nfu_status *status)
double ptwXY_getXMin(ptwXYPoints *ptwXY)
static nfu_status ptwXY_integrateWithFunction3(double x, double *y, void *argList)
nfu_status ptwXY_normalize(ptwXYPoints *ptwXY1)
static nfu_status ptwXY_integrateWithFunction2(nf_Legendre_GaussianQuadrature_callback integrandFunction, void *argList, double x1, double x2, double *integral)
ptwXY_interpolation interpolation
double ptwXY_integrate(ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status)
double ptwXY_integrateDomainWithWeight_sqrt_x(ptwXYPoints *ptwXY, nfu_status *status)
enum ptwXY_interpolation_e ptwXY_interpolation
G4double G4Log(G4double x)
int64_t ptwX_length(ptwXPoints *ptwX)
ptwXPoints * ptwXY_groupOneFunction(ptwXYPoints *ptwXY, ptwXPoints *groupBoundaries, ptwXY_group_normType normType, ptwXPoints *ptwX_norm, nfu_status *status)
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
nfu_status ptwX_setPointAtIndex(ptwXPoints *ptwX, int64_t index, double x)
ptwXPoints * ptwXY_groupTwoFunctions(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, ptwXPoints *groupBoundaries, ptwXY_group_normType normType, ptwXPoints *ptwX_norm, nfu_status *status)
ptwXY_createFromFunction_callback func
double ptwXY_integrateWithWeight_x(ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status)
nfu_status ptwXY_tweakDomainsToMutualify(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, int epsilonFactor, double epsilon)
nfu_status nf_Legendre_GaussianQuadrature(int degree, double x1, double x2, nf_Legendre_GaussianQuadrature_callback func, void *argList, double *integral)
double ptwXY_getXMax(ptwXYPoints *ptwXY)
enum ptwXY_group_normType_e ptwXY_group_normType
nfu_status(* nf_Legendre_GaussianQuadrature_callback)(double x, double *y, void *argList)
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)