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)