11 #if defined __cplusplus 
   19 static nfu_status 
ptwXY_clip2( ptwXYPoints *ptwXY1, 
double y, 
double x1, 
double y1, 
double x2, 
double y2 );
 
   21 static nfu_status 
ptwXY_thin2( ptwXYPoints *thinned, 
char *thin, 
double accuracy, int64_t i1, int64_t i2 );
 
   25 nfu_status 
ptwXY_clip( ptwXYPoints *ptwXY1, 
double yMin, 
double yMax ) {
 
   37     if( ptwXY1->interpolation == ptwXY_interpolationOther ) 
return( nfu_otherInterpolation );
 
   46         y2 = ptwXY1->points[0].y;
 
   48             ptwXY1->points[0].y = yMin; }
 
   49         else if( y2 > yMax ) {
 
   50             ptwXY1->points[0].y = yMax;
 
   53         if( ( clipped = 
ptwXY_new( ptwXY1->interpolation, &(ptwXY1->interpolationOtherInfo), 
 
   54                 ptwXY1->biSectionMax, ptwXY1->accuracy, n, 10, &status, ptwXY1->userFlag ) ) == NULL )
 
   55             return( ptwXY1->status = status );
 
   56         for( i = 0; i < 
n; i++ ) {
 
   57             x2 = ptwXY1->points[i].x;
 
   58             y2 = ptwXY1->points[i].y;
 
   62                     if( points->y > yMin ) {
 
   63                         if( ( status = 
ptwXY_clip2( clipped, yMin, points->x, points->y, x2, y2 ) ) != nfu_Okay ) 
goto Err;
 
   68                 for( i++; i < 
n; i++ ) 
if( !( ptwXY1->points[i].y < yMin ) ) 
break;
 
   70                     x2 = ptwXY1->points[i].x;
 
   71                     y2 = ptwXY1->points[i].y;
 
   72                     if( ( status = 
ptwXY_clip2( clipped, yMin, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) 
goto Err;
 
   74                         if( ( status = 
ptwXY_clip2( clipped, yMax, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) 
goto Err;
 
   76                 else if( j != n - 1 ) {
 
   77                     if( ( status = 
ptwXY_setValueAtX( clipped, ptwXY1->points[n - 1].x, yMin ) ) != nfu_Okay ) 
goto Err;
 
   80             else if( y2 > yMax ) {
 
   83                     if( points->y < yMax ) {
 
   84                         if( ( status = 
ptwXY_clip2( clipped, yMax, points->x, points->y, x2, y2 ) ) != nfu_Okay ) 
goto Err;
 
   89                 for( i++; i < 
n; i++ ) 
if( !( ptwXY1->points[i].y > yMax ) ) 
break;
 
   91                     x2 = ptwXY1->points[i].x;
 
   92                     y2 = ptwXY1->points[i].y;
 
   93                     if( ( status = 
ptwXY_clip2( clipped, yMax, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) 
goto Err;
 
   95                         if( ( status = 
ptwXY_clip2( clipped, yMin, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) 
goto Err;
 
   97                 else if( j != n - 1 ) {
 
   98                     if( ( status = 
ptwXY_setValueAtX( clipped, ptwXY1->points[n - 1].x, yMax ) ) != nfu_Okay ) 
goto Err;
 
  106         ptwXY1->length = clipped->length;   
 
  108         n = ptwXY1->allocatedSize;
 
  109         ptwXY1->allocatedSize = clipped->allocatedSize;
 
  110         clipped->allocatedSize = 
n;
 
  111         points = clipped->points;
 
  112         clipped->points = ptwXY1->points;
 
  113         ptwXY1->points = points;
 
  117     return( ptwXY1->status );
 
  121     return( ptwXY1->status = status );
 
  126 static nfu_status 
ptwXY_clip2( ptwXYPoints *clipped, 
double y, 
double x1, 
double y1, 
double x2, 
double y2 ) {
 
  129     nfu_status status = nfu_Okay;
 
  131     x = ( y - y1 ) * ( x2 - x1 ) / ( y2 - y1 ) + x1;
 
  144 nfu_status 
ptwXY_thicken( ptwXYPoints *ptwXY1, 
int sectionSubdivideMax, 
double dxMax, 
double fxMax ) {
 
  146     double x1, x2 = 0., y1, y2 = 0., fx = 1.1, 
x, dx, dxp, lfx, y;    
 
  147     int64_t i, notFirstPass = 0;
 
  148     int nfx, nDone, doLinear;
 
  151     if( ptwXY1->interpolation == ptwXY_interpolationOther ) 
return( nfu_otherInterpolation );
 
  152     if( ( sectionSubdivideMax < 1 ) || ( dxMax < 0. ) || ( fxMax < 1. ) ) 
return( nfu_badInput );
 
  153     if( sectionSubdivideMax > ptwXY_sectionSubdivideMax ) sectionSubdivideMax = ptwXY_sectionSubdivideMax;
 
  155     for( i = ptwXY1->length - 1; i >= 0; i-- ) {
 
  156         x1 = ptwXY1->points[i].x;
 
  157         y1 = ptwXY1->points[i].y;
 
  168                         nfx = sectionSubdivideMax; }
 
  170                         nfx = ( (int) ( lfx / 
G4Log( fxMax ) ) ) + 1;
 
  171                         if( nfx > sectionSubdivideMax ) nfx = sectionSubdivideMax;
 
  173                     if( nfx > 0 ) fx = 
G4Exp( lfx / nfx );
 
  175                     if( dx < ( fx - 1 ) * x1 ) doLinear = 1; }
 
  188                     if( dx <= ( fx - 1 ) * 
x ) {
 
  193                     dxp = ( fx - 1. ) * 
x;
 
  196                 if( ( x2 - 
x ) < 0.05 * std::fabs( dxp ) ) 
break;
 
  197                 if( ( status = 
ptwXY_interpolatePoint( ptwXY1->interpolation, 
x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) 
return( status );
 
  211 ptwXYPoints *
ptwXY_thin( ptwXYPoints *ptwXY1, 
double accuracy, nfu_status *status ) {
 
  213     int64_t i, j, length = ptwXY1->length;
 
  214     ptwXYPoints *thinned = NULL;
 
  218     if( length < 3 ) 
return( 
ptwXY_clone( ptwXY1, status ) );   
 
  220     *status = nfu_otherInterpolation;
 
  221     if( ptwXY1->interpolation == ptwXY_interpolationOther ) 
return( NULL );
 
  222     if( accuracy < ptwXY1->accuracy ) accuracy = ptwXY1->accuracy;
 
  223     if( ( thinned = 
ptwXY_new( ptwXY1->interpolation, &(ptwXY1->interpolationOtherInfo), 
 
  224         ptwXY1->biSectionMax, accuracy, length, ptwXY1->overflowLength, status, ptwXY1->userFlag ) ) == NULL ) 
return( NULL );
 
  226     thinned->points[0] = ptwXY1->points[0];                     
 
  227     y1 = ptwXY1->points[0].y;
 
  228     y2 = ptwXY1->points[1].y;
 
  229     for( i = 2, j = 1; i < length; i++ ) {
 
  230         y3 = ptwXY1->points[i].y;
 
  231         if( ( y1 != y2 ) || ( y2 != y3 ) ) {
 
  232             thinned->points[j++] = ptwXY1->points[i - 1];
 
  237     thinned->points[j++] = ptwXY1->points[length - 1];
 
  239     if( ptwXY1->interpolation != ptwXY_interpolationFlat ) {    
 
  240         length = thinned->length = j;
 
  241         if( ( thin = (
char *) 
nfu_calloc( 1, (
size_t) length ) ) == NULL ) 
goto Err;
 
  242         if( ( *status = 
ptwXY_thin2( thinned, thin, accuracy, 0, length - 1 ) ) != nfu_Okay ) 
goto Err;
 
  243         for( j = 1; j < length; j++ ) 
if( thin[j] != 0 ) 
break;
 
  244         for( i = j + 1; i < length; i++ ) {
 
  246                 thinned->points[j] = thinned->points[i];
 
  258     if( thin != NULL ) 
nfu_free( thin );
 
  264 static nfu_status 
ptwXY_thin2( ptwXYPoints *thinned, 
char *thin, 
double accuracy, int64_t i1, int64_t i2 ) {
 
  267     double y, 
s, dY, dYMax = 0., dYR, dYRMax = 0;
 
  268     double x1 = thinned->points[i1].x, y1 = thinned->points[i1].y, x2 = thinned->points[i2].x, y2 = thinned->points[i2].y;
 
  271     if( i1 + 1 >= i2 ) 
return( nfu_Okay );
 
  272     for( i = i1 + 1; i < i2; i++ ) {
 
  273         if( ( status = 
ptwXY_interpolatePoint( thinned->interpolation, thinned->points[i].x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) 
return( status );
 
  274         s = 0.5 * ( std::fabs( y ) + std::fabs( thinned->points[i].y ) );
 
  275         dY = std::fabs( y - thinned->points[i].y );
 
  277         if( s != 0 ) dYR = dY / 
s;
 
  278         if( ( dYR > dYRMax ) || ( ( dYR >= 0.9999 * dYRMax ) && ( dY > dYMax ) ) ) {    
 
  280             if( dY > dYMax ) dYMax = dY;
 
  281             if( dYR > dYRMax ) dYRMax = dYR;
 
  284     if( dYRMax < accuracy ) {
 
  285         for( i = i1 + 1; i < i2; i++ ) thin[i] = 1; }
 
  287         if( ( status = 
ptwXY_thin2( thinned, thin, accuracy, i1, iMax ) ) != nfu_Okay ) 
return( status );
 
  288         status = 
ptwXY_thin2( thinned, thin, accuracy, iMax, i2 );
 
  298     double dx = x2 - x1, dndx;
 
  301         dx = ( x2 - x1 ) / sectionSubdivideMax; }
 
  305         if( ( dndx  - ndx ) > 1e-6 ) ndx++;
 
  306         if( ndx > sectionSubdivideMax ) ndx = sectionSubdivideMax;
 
  307         if( ndx > 0 ) dx /= ndx;
 
  323     if( ptwXY->status != nfu_Okay ) 
return( ptwXY->status );
 
  325     for( i1 = 0; i1 < ptwXY->length; i1++ ) {
 
  326         if( ptwXY->points[i1].y != 0 ) 
break;
 
  329     for( i2 = ptwXY->length - 1; i2 >= 0; i2-- ) {
 
  330         if( ptwXY->points[i2].y != 0 ) 
break;
 
  333     if( i2 < ptwXY->length ) i2++;
 
  336             for( i = i1; i < i2; i++ ) ptwXY->points[i - i1] = ptwXY->points[i];
 
  338         ptwXY->length = i2 - i1; }
 
  340         ptwXY->points[1] = ptwXY->points[ptwXY->length - 1];
 
  349 ptwXYPoints *
ptwXY_union( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, 
int unionOptions ) {
 
  351     int64_t overflowSize, i, i1 = 0, i2 = 0, n1 = ptwXY1->length, n2 = ptwXY2->length, length;
 
  352     int fillWithFirst = unionOptions & ptwXY_union_fill, trim = unionOptions & ptwXY_union_trim;
 
  354     double x1 = 0., x2 = 0., y1 = 0., y2 = 0., y, biSectionMax, accuracy;
 
  356     if( ( *status = ptwXY1->status ) != nfu_Okay ) 
return( NULL );
 
  357     if( ( *status = ptwXY2->status ) != nfu_Okay ) 
return( NULL );
 
  358     *status = nfu_otherInterpolation;
 
  359     if( ptwXY1->interpolation == ptwXY_interpolationOther ) 
return( NULL );
 
  366     if( ( n1 == 1 ) || ( n2 == 1 ) ) {
 
  367         *status = nfu_tooFewPoints;
 
  373                 if( ptwXY1->points[0].x < ptwXY2->points[0].x ) {
 
  375                         if( ptwXY1->points[i1].x >= ptwXY2->points[0].x ) 
break;
 
  376                         if( fillWithFirst ) {
 
  377                             if( i1 < ( ptwXY1->length - 1 ) ) {
 
  378                                 x1 = ptwXY1->points[i1].x;
 
  379                                 y1 = ptwXY1->points[i1].y;
 
  380                                 x2 = ptwXY1->points[i1+1].x;
 
  381                                 y2 = ptwXY1->points[i1+1].y; 
 
  388                         if( ptwXY2->points[i2].x >= ptwXY1->points[0].x ) 
break;
 
  392                 if( ptwXY1->points[n1-1].x > ptwXY2->points[n2-1].x ) {
 
  394                         if( ptwXY1->points[n1-1].x <= ptwXY2->points[n2-1].x ) 
break;
 
  399                         if( ptwXY2->points[n2-1].x <= ptwXY1->points[n1-1].x ) 
break;
 
  410     overflowSize = ptwXY1->overflowAllocatedSize;
 
  411     if( overflowSize < ptwXY2->overflowAllocatedSize ) overflowSize = ptwXY2->overflowAllocatedSize;
 
  412     length = ( n1 - i1 ) + ( n2 - i2 );
 
  413     if( length == 0 ) length = ptwXY_minimumSize;
 
  414     biSectionMax = ptwXY1->biSectionMax;
 
  415     if( biSectionMax < ptwXY2->biSectionMax ) biSectionMax = ptwXY2->biSectionMax;
 
  416     accuracy = ptwXY1->accuracy;
 
  417     if( accuracy < ptwXY2->accuracy ) accuracy = ptwXY2->accuracy;
 
  418     n = 
ptwXY_new( ptwXY1->interpolation, NULL, biSectionMax, accuracy, length, overflowSize, status, ptwXY1->userFlag );
 
  419     if( n == NULL ) 
return( NULL );
 
  421     for( i = 0; ( i1 < n1 ) && ( i2 < n2 ); i++ ) {
 
  423         if( ptwXY1->points[i1].x <= ptwXY2->points[i2].x ) {
 
  424             n->points[i].x = ptwXY1->points[i1].x;
 
  425             if( fillWithFirst ) {
 
  426                 y = ptwXY1->points[i1].y;
 
  427                 if( i1 < ( ptwXY1->length - 1 ) ) {
 
  428                     x1 = ptwXY1->points[i1].x;
 
  429                     y1 = ptwXY1->points[i1].y;
 
  430                     x2 = ptwXY1->points[i1+1].x;
 
  431                     y2 = ptwXY1->points[i1+1].y; }
 
  437             if( ptwXY1->points[i1].x == ptwXY2->points[i2].x ) i2++;
 
  440             n->points[i].x = ptwXY2->points[i2].x;
 
  441             if( fillWithFirst && ( ( y1 != 0. ) || ( y2 != 0. ) ) ) {
 
  442                 if( ( *status = 
ptwXY_interpolatePoint( ptwXY1->interpolation, ptwXY2->points[i2].x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) {
 
  453     for( ; i1 < n1; i1++, i++ ) {
 
  454         n->points[i].x = ptwXY1->points[i1].x;
 
  455         if( fillWithFirst ) y = ptwXY1->points[i1].y;
 
  458     for( ; i2 < n2; i2++, i++ ) {
 
  459         n->points[i].x = ptwXY2->points[i2].x;
 
  460         if( fillWithFirst && trim && ( n->points[i].x <= x2 ) ) {
 
  461             if( ( *status = 
ptwXY_interpolatePoint( ptwXY1->interpolation, n->points[i].x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) {
 
  470     if( unionOptions & ptwXY_union_mergeClosePoints ) {
 
  483     int64_t i1, length = ptwXY->length;
 
  487     if( ptwXY->status != nfu_Okay ) 
return( ptwXY->status );
 
  488     if( xScale == 0 ) 
return( nfu_XNotAscending );
 
  492     for( i1 = 0, p1 = ptwXY->points; i1 < length; i1++, p1++ ) {
 
  493         p1->x = xScale * p1->x + xOffset;
 
  494         p1->y = yScale * p1->y + yOffset;
 
  498         int64_t length_2 = length / 2;
 
  501         for( i1 = 0, p1 = ptwXY->points, p2 = &(ptwXY->points[length-1]); i1 < length_2; i1++ ) {
 
  508     return( ptwXY->status );
 
  511 #if defined __cplusplus 
nfu_status ptwXY_mergeClosePoints(ptwXYPoints *ptwXY, double epsilon)
nfu_status ptwXY_thicken(ptwXYPoints *ptwXY1, int sectionSubdivideMax, double dxMax, double fxMax)
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
static double ptwXY_thicken_linear_dx(int sectionSubdivideMax, double dxMax, double x1, double x2)
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
nfu_status ptwXY_clip(ptwXYPoints *ptwXY1, double yMin, double yMax)
static nfu_status ptwXY_clip2(ptwXYPoints *ptwXY1, double y, double x1, double y1, double x2, double y2)
ptwXYPoints * ptwXY_thin(ptwXYPoints *ptwXY1, double accuracy, nfu_status *status)
double ptwXY_getYMax(ptwXYPoints *ptwXY)
double ptwXY_getYMin(ptwXYPoints *ptwXY)
G4double G4Log(G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision. 
ptwXYPoints * ptwXY_clone(ptwXYPoints *ptwXY, nfu_status *status)
ptwXYPoint * ptwXY_getPointAtIndex_Unsafely(ptwXYPoints *ptwXY, int64_t index)
nfu_status ptwXY_clear(ptwXYPoints *ptwXY)
static nfu_status ptwXY_thin2(ptwXYPoints *thinned, char *thin, double accuracy, int64_t i1, int64_t i2)
nfu_status ptwXY_scaleOffsetXAndY(ptwXYPoints *ptwXY, double xScale, double xOffset, double yScale, double yOffset)
const G4double x[NPOINTSGL]
ptwXYPoints * ptwXY_union(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int unionOptions)
nfu_status ptwXY_trim(ptwXYPoints *ptwXY)
ptwXYPoints * ptwXY_new(ptwXY_interpolation interpolation, ptwXY_interpolationOtherInfo const *interpolationOtherInfo, double biSectionMax, double accuracy, int64_t primarySize, int64_t secondarySize, nfu_status *status, int userFlag)
nfu_status ptwXY_interpolatePoint(ptwXY_interpolation interpolation, double x, double *y, double x1, double y1, double x2, double y2)
void * nfu_calloc(size_t size, size_t n)