12 #if defined __cplusplus
23 static void ptwXY_initialOverflowPoint( ptwXYOverflowPoint *overflowPoint, ptwXYOverflowPoint *prior, ptwXYOverflowPoint *next );
24 static nfu_status
ptwXY_mergeFrom( ptwXYPoints *ptwXY,
int incY,
int length,
double *xs,
double *ys );
29 ptwXYPoints *
ptwXY_new( ptwXY_interpolation interpolation, ptwXY_interpolationOtherInfo
const *interpolationOtherInfo,
double biSectionMax,
30 double accuracy, int64_t primarySize, int64_t secondarySize, nfu_status *status,
int userFlag ) {
32 ptwXYPoints *ptwXY = (ptwXYPoints *)
nfu_calloc(
sizeof( ptwXYPoints ), 1 );
34 *status = nfu_mallocError;
35 if( ptwXY == NULL )
return( NULL );
36 ptwXY_setup( ptwXY, interpolation, interpolationOtherInfo, biSectionMax, accuracy, primarySize,
37 secondarySize, userFlag );
38 if( ( *status = ptwXY->status ) != nfu_Okay ) {
39 ptwXY = (ptwXYPoints *)
nfu_free( ptwXY );
46 nfu_status
ptwXY_setup( ptwXYPoints *ptwXY, ptwXY_interpolation interpolation, ptwXY_interpolationOtherInfo
const *interpolationOtherInfo,
47 double biSectionMax,
double accuracy, int64_t primarySize, int64_t secondarySize,
int userFlag ) {
49 ptwXY->status = nfu_Okay;
50 ptwXY->typeX = ptwXY_sigma_none;
51 ptwXY->typeY = ptwXY_sigma_none;
52 ptwXY->interpolation = interpolation;
53 ptwXY->interpolationOtherInfo.interpolationString = NULL;
54 ptwXY->interpolationOtherInfo.getValueFunc = NULL;
55 ptwXY->interpolationOtherInfo.argList = NULL;
56 switch( interpolation ) {
57 case ptwXY_interpolationLinLin :
59 case ptwXY_interpolationLinLog :
61 case ptwXY_interpolationLogLin :
63 case ptwXY_interpolationLogLog :
65 case ptwXY_interpolationFlat :
67 case ptwXY_interpolationOther :
68 if( interpolationOtherInfo == NULL ) {
69 ptwXY->status = nfu_otherInterpolation; }
71 if( interpolationOtherInfo->interpolationString == NULL ) {
72 ptwXY->status = nfu_otherInterpolation; }
74 if( ( ptwXY->interpolationOtherInfo.interpolationString = strdup( interpolationOtherInfo->interpolationString ) ) == NULL ) {
75 ptwXY->status = nfu_mallocError;
78 ptwXY->interpolationOtherInfo.getValueFunc = interpolationOtherInfo->getValueFunc;
79 ptwXY->interpolationOtherInfo.argList = interpolationOtherInfo->argList;
84 ptwXY->biSectionMax = ptwXY_maxBiSectionMax;
86 ptwXY->accuracy = ptwXY_minAccuracy;
90 ptwXY->allocatedSize = 0;
91 ptwXY->overflowLength = 0;
92 ptwXY->overflowAllocatedSize = 0;
93 ptwXY->mallocFailedSize = 0;
98 ptwXY->overflowPoints = NULL;
103 return( ptwXY->status );
108 ptwXYPoints *
ptwXY_create( ptwXY_interpolation interpolation, ptwXY_interpolationOtherInfo
const *interpolationOtherInfo,
109 double biSectionMax,
double accuracy, int64_t primarySize, int64_t secondarySize, int64_t length,
double const *xy,
110 nfu_status *status,
int userFlag ) {
114 if( primarySize < length ) primarySize = length;
115 if( ( ptwXY =
ptwXY_new( interpolation, interpolationOtherInfo, biSectionMax, accuracy, primarySize,
116 secondarySize, status, userFlag ) ) != NULL ) {
126 ptwXYPoints *
ptwXY_createFrom_Xs_Ys( ptwXY_interpolation interpolation, ptwXY_interpolationOtherInfo
const *interpolationOtherInfo,
127 double biSectionMax,
double accuracy, int64_t primarySize, int64_t secondarySize, int64_t length,
double const *Xs,
128 double const *Ys, nfu_status *status,
int userFlag ) {
133 if( primarySize < length ) primarySize = length;
134 if( ( ptwXY =
ptwXY_new( interpolation, interpolationOtherInfo, biSectionMax, accuracy, primarySize,
135 secondarySize, status, userFlag ) ) != NULL ) {
136 for( i = 0; i < length; i++ ) {
137 ptwXY->points[i].x = Xs[i];
138 ptwXY->points[i].y = Ys[i];
140 ptwXY->length = length;
148 nfu_status
ptwXY_copy( ptwXYPoints *dest, ptwXYPoints *src ) {
151 ptwXYPoint *pointFrom, *pointTo;
152 ptwXYOverflowPoint *o, *overflowHeader = &(src->overflowHeader);
154 if( dest->status != nfu_Okay )
return( dest->status );
155 if( src->status != nfu_Okay )
return( src->status );
158 if( dest->interpolation == ptwXY_interpolationOther ) {
159 if( dest->interpolationOtherInfo.interpolationString != NULL ) {
160 dest->interpolationOtherInfo.interpolationString = (
char const *)
nfu_free( (
void *) dest->interpolationOtherInfo.interpolationString );
163 dest->interpolation = ptwXY_interpolationLinLin;
165 if( dest->status != nfu_Okay )
return( dest->status );
166 dest->interpolation = src->interpolation;
167 if( dest->interpolation == ptwXY_interpolationOther ) {
168 if( src->interpolationOtherInfo.interpolationString != NULL ) {
169 if( ( dest->interpolationOtherInfo.interpolationString = strdup( src->interpolationOtherInfo.interpolationString ) ) == NULL )
170 return( dest->status = nfu_mallocError );
173 dest->interpolationOtherInfo.interpolationString = src->interpolationOtherInfo.interpolationString;
175 dest->interpolationOtherInfo.getValueFunc = src->interpolationOtherInfo.getValueFunc;
176 dest->interpolationOtherInfo.argList = src->interpolationOtherInfo.argList;
177 dest->userFlag = src->userFlag;
178 dest->biSectionMax = src->biSectionMax;
179 dest->accuracy = src->accuracy;
180 dest->minFractional_dx = src->minFractional_dx;
181 pointFrom = src->points;
182 o = src->overflowHeader.next;
183 pointTo = dest->points;
185 while( o != overflowHeader ) {
186 if( i < nonOverflowLength ) {
187 if( pointFrom->x < o->point.x ) {
188 *pointTo = *pointFrom;
201 for( ; i < nonOverflowLength; i++, pointFrom++, pointTo++ ) *pointTo = *pointFrom;
202 dest->length = src->length;
203 return( dest->status );
208 ptwXYPoints *
ptwXY_clone( ptwXYPoints *ptwXY, nfu_status *status ) {
210 return(
ptwXY_slice( ptwXY, 0, ptwXY->length, ptwXY->overflowAllocatedSize, status ) );
219 if( interpolationTo == ptwXY_interpolationOther ) {
220 *status = nfu_otherInterpolation;
223 if( ( n1 =
ptwXY_clone( ptwXY, status ) ) != NULL ) {
224 if( n1->interpolation == ptwXY_interpolationOther )
nfu_free( (
void *) n1->interpolationOtherInfo.interpolationString );
225 n1->interpolation = interpolationTo;
226 switch( interpolationTo ) {
227 case ptwXY_interpolationLinLin :
229 case ptwXY_interpolationLinLog :
231 case ptwXY_interpolationLogLin :
233 case ptwXY_interpolationLogLog :
235 case ptwXY_interpolationFlat :
237 case ptwXY_interpolationOther :
240 n1->interpolationOtherInfo.getValueFunc = NULL;
241 n1->interpolationOtherInfo.argList = NULL;
248 ptwXYPoints *
ptwXY_slice( ptwXYPoints *ptwXY, int64_t index1, int64_t index2, int64_t secondarySize, nfu_status *status ) {
253 *status = nfu_badSelf;
254 if( ptwXY->status != nfu_Okay )
return( NULL );
256 *status = nfu_badIndex;
257 if( index2 < index1 )
return( NULL );
258 if( index1 < 0 ) index1 = 0;
259 if( index2 > ptwXY->length ) index2 = ptwXY->length;
261 length = index2 - index1;
263 if( ( n =
ptwXY_new( ptwXY->interpolation, &(ptwXY->interpolationOtherInfo), ptwXY->biSectionMax,
264 ptwXY->accuracy, length, secondarySize, status, ptwXY->userFlag ) ) == NULL )
return( NULL );
266 *status = n->status = ptwXY->status;
267 for( i = index1; i < index2; i++ ) n->points[i - index1] = ptwXY->points[i];
274 ptwXYPoints *
ptwXY_xSlice( ptwXYPoints *ptwXY,
double xMin,
double xMax, int64_t secondarySize,
int fill, nfu_status *status ) {
278 ptwXYPoints *
n = NULL;
280 if( ( *status = ptwXY->status ) != nfu_Okay )
return( NULL );
283 n =
ptwXY_new( ptwXY->interpolation, &(ptwXY->interpolationOtherInfo), ptwXY->biSectionMax,
284 ptwXY->accuracy, 0, secondarySize, status, ptwXY->userFlag ); }
286 if( ( n =
ptwXY_clone( ptwXY, status ) ) == NULL )
return( NULL );
287 if( ( n->points[0].x < xMin ) || ( n->points[n->length - 1].x > xMax ) ) {
288 if( fill && ( n->points[n->length - 1].x > xMax ) ) {
292 if( fill && ( n->points[0].x < xMin ) ) {
297 for( i1 = 0; i1 < n->length; i1++ )
if( n->points[i1].x >= xMin )
break;
298 for( i2 = n->length - 1; i2 > 0; i2-- )
if( n->points[i2].x <= xMax )
break;
301 for( i = i1; i < i2; i++ ) n->points[i- i1] = n->points[i];
315 ptwXYPoints *
ptwXY_xMinSlice( ptwXYPoints *ptwXY,
double xMin, int64_t secondarySize,
int fill, nfu_status *status ) {
317 double xMax = 1.1 * xMin + 1;
319 if( xMin < 0 ) xMax = 0.9 * xMin + 1;
321 return(
ptwXY_xSlice( ptwXY, xMin, xMax, secondarySize, fill, status ) );
326 ptwXYPoints *
ptwXY_xMaxSlice( ptwXYPoints *ptwXY,
double xMax, int64_t secondarySize,
int fill, nfu_status *status ) {
328 double xMin = 0.9 * xMax - 1;
330 if( xMax < 0 ) xMin = 1.1 * xMax - 1;
332 return(
ptwXY_xSlice( ptwXY, xMin, xMax, secondarySize, fill, status ) );
339 return( ptwXY->interpolation );
346 return( ptwXY->interpolationOtherInfo.interpolationString );
353 return( ptwXY->status );
360 return( ptwXY->userFlag );
367 ptwXY->userFlag = userFlag;
374 return( ptwXY->accuracy );
381 if( accuracy < ptwXY_minAccuracy ) accuracy = ptwXY_minAccuracy;
382 if( accuracy < ptwXY->accuracy ) accuracy = ptwXY->accuracy;
383 if( accuracy > 1 ) accuracy = 1.;
384 ptwXY->accuracy = accuracy;
385 return( ptwXY->accuracy );
392 return( ptwXY->biSectionMax );
399 if( biSectionMax < 0 ) {
401 else if( biSectionMax > ptwXY_maxBiSectionMax ) {
402 biSectionMax = ptwXY_maxBiSectionMax;
404 ptwXY->biSectionMax = biSectionMax;
405 return( ptwXY->biSectionMax );
414 if( ptwXY->status != nfu_Okay )
return( ptwXY->status );
416 if( size < ptwXY_minimumSize ) size = ptwXY_minimumSize;
417 if( size < ptwXY->length ) size = ptwXY->length;
418 if( size != ptwXY->allocatedSize ) {
419 if( size > ptwXY->allocatedSize ) {
420 ptwXY->points = (ptwXYPoint *)
nfu_realloc( (
size_t) size *
sizeof( ptwXYPoint ), ptwXY->points ); }
421 else if( ( ptwXY->allocatedSize > 2 * size ) || forceSmallerResize ) {
422 ptwXY->points = (ptwXYPoint *)
nfu_realloc( (
size_t) size *
sizeof( ptwXYPoint ), ptwXY->points ); }
424 size = ptwXY->allocatedSize;
426 if( ptwXY->points == NULL ) {
428 ptwXY->mallocFailedSize = size;
430 ptwXY->status = nfu_mallocError;
432 ptwXY->allocatedSize = size;
434 return( ptwXY->status );
443 nfu_status status = nfu_Okay;
445 if( ptwXY->status != nfu_Okay )
return( ptwXY->status );
447 if( size < ptwXY_minimumOverflowSize ) size = ptwXY_minimumOverflowSize;
448 if( size < ptwXY->overflowLength ) status =
ptwXY_coalescePoints( ptwXY, ptwXY->length + ptwXY->overflowAllocatedSize, NULL, 0 );
449 if( status == nfu_Okay ) {
450 if( size != ptwXY->overflowAllocatedSize ) {
451 ptwXY->overflowPoints = (ptwXYOverflowPoint *)
nfu_realloc( (
size_t) size *
sizeof( ptwXYOverflowPoint ), ptwXY->overflowPoints );
452 if( ptwXY->overflowPoints == NULL ) {
454 ptwXY->overflowLength = 0;
455 ptwXY->mallocFailedSize = size;
457 ptwXY->status = nfu_mallocError;
460 ptwXY->overflowAllocatedSize = size; }
462 ptwXY->status = status;
464 return( ptwXY->status );
469 nfu_status
ptwXY_coalescePoints( ptwXYPoints *ptwXY, int64_t size, ptwXYPoint *newPoint,
int forceSmallerResize ) {
472 int64_t length = ptwXY->length + ( ( newPoint != NULL ) ? 1 : 0 );
473 ptwXYOverflowPoint *last = ptwXY->overflowHeader.prior;
474 ptwXYPoint *pointsFrom, *pointsTo;
476 if( ptwXY->status != nfu_Okay )
return( ptwXY->status );
477 if( ptwXY->overflowLength == 0 )
return( nfu_Okay );
479 if( size < length ) size = length;
480 if( size > ptwXY->allocatedSize ) {
484 pointsTo = &(ptwXY->points[length - 1]);
485 while( last != &(ptwXY->overflowHeader) ) {
487 if( newPoint != NULL ) {
488 if( ( pointsFrom >= ptwXY->points ) && ( pointsFrom->x > last->point.x ) ) {
489 if( newPoint->x > pointsFrom->x ) addNewPoint = 1; }
491 if( newPoint->x > last->point.x ) addNewPoint = 1;
493 if( addNewPoint == 1 ) {
494 *pointsTo = *newPoint;
498 if( addNewPoint == 0 ) {
499 if( ( pointsFrom >= ptwXY->points ) && ( pointsFrom->x > last->point.x ) ) {
500 *pointsTo = *pointsFrom;
503 *pointsTo = last->point;
509 while( ( newPoint != NULL ) && ( pointsFrom >= ptwXY->points ) ) {
510 if( newPoint->x > pointsFrom->x ) {
511 *pointsTo = *newPoint;
514 *pointsTo = *pointsFrom;
519 if( newPoint != NULL ) *pointsTo = *newPoint;
520 ptwXY->overflowHeader.prior = &(ptwXY->overflowHeader);
521 ptwXY->overflowHeader.next = &(ptwXY->overflowHeader);
522 ptwXY->length = length;
523 ptwXY->overflowLength = 0;
538 if( ptwXY->status != nfu_Okay )
return( ptwXY->status );
541 ptwXY->overflowLength = 0;
542 ptwXY->overflowHeader.prior = &(ptwXY->overflowHeader);
543 ptwXY->overflowHeader.next = &(ptwXY->overflowHeader);
554 if( ptwXY->interpolation == ptwXY_interpolationOther ) {
555 if( ptwXY->interpolationOtherInfo.interpolationString != NULL )
556 ptwXY->interpolationOtherInfo.interpolationString = (
char const *)
nfu_free( (
void *) ptwXY->interpolationOtherInfo.interpolationString );
558 ptwXY->interpolation = ptwXY_interpolationLinLin;
559 ptwXY->interpolationOtherInfo.getValueFunc = NULL;
560 ptwXY->interpolationOtherInfo.argList = NULL;
562 ptwXY->allocatedSize = 0;
563 ptwXY->points = (ptwXYPoint *)
nfu_free( ptwXY->points );
565 ptwXY->overflowLength = 0;
566 ptwXY->overflowAllocatedSize = 0;
567 ptwXY->overflowPoints = (ptwXYOverflowPoint *)
nfu_free( ptwXY->overflowPoints );
578 return( (ptwXYPoints *) NULL );
585 return( ptwXY->length );
592 return( ptwXY->length - ptwXY->overflowLength );
599 nfu_status status = nfu_Okay;
602 double const *d = xy;
605 if( length > ptwXY->allocatedSize ) {
607 if( status != nfu_Okay )
return( status );
609 for( i = 0, p = ptwXY->points; i < length; i++, p++ ) {
612 status = nfu_XNotAscending;
613 ptwXY->status = nfu_XNotAscending;
622 ptwXY->overflowHeader.next = &(ptwXY->overflowHeader);
623 ptwXY->overflowHeader.prior = &(ptwXY->overflowHeader);
624 ptwXY->overflowLength = 0;
625 ptwXY->length = length;
626 return( ptwXY->status = status );
638 if( ( status =
ptwXY_clear( ptwXY ) ) != nfu_Okay )
return( status );
639 if( length > ptwXY->allocatedSize ) {
642 for( i = 0, p = ptwXY->points; i < length; i++, p++, x++, y++ ) {
645 status = ptwXY->status = nfu_XNotAscending;
654 ptwXY->length = length;
662 int64_t
n = ptwXY->length - ( i2 - i1 );
665 if( ( i1 < 0 ) || ( i1 > i2 ) || ( i2 > ptwXY->length ) )
return( nfu_badIndex );
667 for( ; i2 < ptwXY->length; i1++, i2++ ) ptwXY->points[i1] = ptwXY->points[i2];
670 return( ptwXY->status );
677 if( ptwXY->status != nfu_Okay )
return( NULL );
678 if( ( index < 0 ) || ( index >= ptwXY->length ) )
return( NULL );
687 ptwXYOverflowPoint *overflowPoint;
689 for( overflowPoint = ptwXY->overflowHeader.next, i = 0; overflowPoint != &(ptwXY->overflowHeader); overflowPoint = overflowPoint->next, i++ ) {
690 if( overflowPoint->index == index )
return( &(overflowPoint->point) );
691 if( overflowPoint->index > index )
break;
693 return( &(ptwXY->points[index - i]) );
702 if( p == NULL )
return( nfu_badIndex );
710 ptwXY_lessEqualGreaterX
ptwXY_getPointsAroundX( ptwXYPoints *ptwXY,
double x, ptwXYOverflowPoint *lessThanEqualXPoint, ptwXYOverflowPoint *greaterThanXPoint ) {
713 ptwXYPoint *closePoint;
721 ptwXYOverflowPoint *greaterThanXPoint,
double eps,
int *closeIsEqual, ptwXYPoint **closePoint ) {
724 int64_t indexMin, indexMid, indexMax;
725 ptwXY_dataFrom xMinFrom, xMaxFrom;
727 ptwXYOverflowPoint *overflowPoint, *overflowHeader = &(ptwXY->overflowHeader);
728 ptwXY_lessEqualGreaterX status = ptwXY_lessEqualGreaterX_empty;
729 ptwXYPoint *lowerPoint = NULL, *upperPoint = NULL;
733 if( ptwXY->length != 0 ) {
735 status = ptwXY_lessEqualGreaterX_lessThan;
736 if( xMinFrom == ptwXY_dataFrom_Points ) {
737 greaterThanXPoint->prior = overflowHeader;
738 greaterThanXPoint->index = 0;
739 greaterThanXPoint->point = ptwXY->points[0];
740 *closePoint = &(ptwXY->points[0]); }
742 *greaterThanXPoint = *(overflowHeader->next);
743 *closePoint = &(overflowHeader->next->point);
745 else if( x > xMax ) {
746 status = ptwXY_lessEqualGreaterX_greater;
747 if( xMaxFrom == ptwXY_dataFrom_Points ) {
748 lessThanEqualXPoint->prior = overflowHeader->prior;
749 lessThanEqualXPoint->index = nonOverflowLength - 1;
750 lessThanEqualXPoint->point = ptwXY->points[lessThanEqualXPoint->index];
751 *closePoint = &(ptwXY->points[lessThanEqualXPoint->index]); }
753 *lessThanEqualXPoint = *(overflowHeader->prior);
754 *closePoint = &(overflowHeader->prior->point);
757 status = ptwXY_lessEqualGreaterX_between;
758 for( overflowPoint = overflowHeader->next, overflowIndex = 0; overflowPoint != overflowHeader;
759 overflowPoint = overflowPoint->next, overflowIndex++ )
if( overflowPoint->point.x > x )
break;
760 overflowPoint = overflowPoint->prior;
761 if( ( overflowPoint != overflowHeader ) && ( overflowPoint->point.x ==
x ) ) {
762 status = ptwXY_lessEqualGreaterX_equal;
763 *lessThanEqualXPoint = *overflowPoint; }
764 else if( ptwXY->length == 1 ) {
765 status = ptwXY_lessEqualGreaterX_equal;
766 lessThanEqualXPoint->index = 0;
767 lessThanEqualXPoint->point = ptwXY->points[0]; }
770 indexMax = nonOverflowLength - 1;
771 indexMid = ( indexMin + indexMax ) >> 1;
772 while( ( indexMin != indexMid ) && ( indexMid != indexMax ) ) {
773 if( ptwXY->points[indexMid].x > x ) {
774 indexMax = indexMid; }
778 indexMid = ( indexMin + indexMax ) >> 1;
780 if( ptwXY->points[indexMin].x == x ) {
781 status = ptwXY_lessEqualGreaterX_equal;
782 lessThanEqualXPoint->index = indexMin;
783 lessThanEqualXPoint->point = ptwXY->points[indexMin]; }
784 else if( ptwXY->points[indexMax].x == x ) {
785 status = ptwXY_lessEqualGreaterX_equal;
786 lessThanEqualXPoint->index = indexMax;
787 lessThanEqualXPoint->point = ptwXY->points[indexMax]; }
789 if( ptwXY->points[indexMin].x > x ) indexMax = 0;
790 if( ptwXY->points[indexMax].x < x ) indexMin = indexMax;
791 if( ( overflowPoint == overflowHeader ) ||
792 ( ( ptwXY->points[indexMin].x > overflowPoint->point.x ) && ( ptwXY->points[indexMin].x < x ) ) ) {
793 if( overflowPoint != overflowHeader ) lessThanEqualXPoint->prior = overflowPoint;
794 lowerPoint = &(ptwXY->points[indexMin]);
795 lessThanEqualXPoint->index = indexMin;
796 lessThanEqualXPoint->point = ptwXY->points[indexMin]; }
798 lowerPoint = &(overflowPoint->point);
799 *lessThanEqualXPoint = *overflowPoint;
801 if( ( overflowPoint->next == overflowHeader ) ||
802 ( ( ptwXY->points[indexMax].x < overflowPoint->next->point.x ) && ( ptwXY->points[indexMax].x > x ) ) ) {
803 upperPoint = &(ptwXY->points[indexMax]);
804 greaterThanXPoint->index = indexMax;
805 greaterThanXPoint->point = ptwXY->points[indexMax]; }
807 upperPoint = &(overflowPoint->next->point);
808 *greaterThanXPoint = *(overflowPoint->next);
817 double absX = std::fabs( x );
819 if( status == ptwXY_lessEqualGreaterX_lessThan ) {
820 if( absX < std::fabs( greaterThanXPoint->point.x ) ) absX = std::fabs( greaterThanXPoint->point.x );
821 if( ( greaterThanXPoint->point.x - x ) < eps * absX ) *closeIsEqual = 1; }
822 else if( status == ptwXY_lessEqualGreaterX_greater ) {
823 if( absX < std::fabs( lessThanEqualXPoint->point.x ) ) absX = std::fabs( lessThanEqualXPoint->point.x );
824 if( ( x - lessThanEqualXPoint->point.x ) < eps * absX ) *closeIsEqual = -1; }
825 else if( status == ptwXY_lessEqualGreaterX_between ) {
826 if( ( x - lessThanEqualXPoint->point.x ) < ( greaterThanXPoint->point.x - x ) ) {
827 *closePoint = lowerPoint;
828 if( absX < std::fabs( lessThanEqualXPoint->point.x ) ) absX = std::fabs( lessThanEqualXPoint->point.x );
829 if( ( x - lessThanEqualXPoint->point.x ) < eps * absX ) *closeIsEqual = -1; }
831 *closePoint = upperPoint;
832 if( absX < std::fabs( greaterThanXPoint->point.x ) ) absX = std::fabs( greaterThanXPoint->point.x );
833 if( ( greaterThanXPoint->point.x - x ) < eps * absX ) *closeIsEqual = 1;
835 else if( status == ptwXY_lessEqualGreaterX_equal ) {
846 nfu_status status = nfu_XOutsideDomain;
847 ptwXYOverflowPoint lessThanEqualXPoint, greaterThanXPoint;
848 ptwXY_lessEqualGreaterX legx =
ptwXY_getPointsAroundX( ptwXY, x, &lessThanEqualXPoint, &greaterThanXPoint );
851 if( ptwXY->status != nfu_Okay )
return( ptwXY->status );
853 case ptwXY_lessEqualGreaterX_empty :
854 case ptwXY_lessEqualGreaterX_lessThan :
855 case ptwXY_lessEqualGreaterX_greater :
857 case ptwXY_lessEqualGreaterX_equal :
859 *y = lessThanEqualXPoint.point.y;
861 case ptwXY_lessEqualGreaterX_between :
862 if( ptwXY->interpolationOtherInfo.getValueFunc != NULL ) {
863 status = ptwXY->interpolationOtherInfo.getValueFunc( ptwXY->interpolationOtherInfo.argList, x, y,
864 lessThanEqualXPoint.point.x, lessThanEqualXPoint.point.y, greaterThanXPoint.point.x, greaterThanXPoint.point.y ); }
866 status =
ptwXY_interpolatePoint( ptwXY->interpolation, x, y, lessThanEqualXPoint.point.x, lessThanEqualXPoint.point.y,
867 greaterThanXPoint.point.x, greaterThanXPoint.point.y );
887 nfu_status status = nfu_Okay;
888 ptwXY_lessEqualGreaterX legx;
889 ptwXYPoint *point = NULL, newPoint = {
x, y };
890 ptwXYOverflowPoint *overflowPoint, *p, *overflowHeader = &(ptwXY->overflowHeader);
891 ptwXYOverflowPoint lessThanEqualXPoint, greaterThanXPoint;
892 ptwXYPoint *closePoint;
894 if( ptwXY->status != nfu_Okay )
return( ptwXY->status );
898 case ptwXY_lessEqualGreaterX_lessThan :
899 case ptwXY_lessEqualGreaterX_greater :
900 case ptwXY_lessEqualGreaterX_between :
902 if( !
override )
return( status );
904 legx = ptwXY_lessEqualGreaterX_equal;
907 if( ( legx == ptwXY_lessEqualGreaterX_greater ) && ( nonOverflowLength < ptwXY->allocatedSize ) ) {
908 point = &(ptwXY->points[nonOverflowLength]); }
910 if( ptwXY->overflowLength == ptwXY->overflowAllocatedSize )
911 return(
ptwXY_coalescePoints( ptwXY, ptwXY->length + ptwXY->overflowAllocatedSize, &newPoint, 0 ) );
912 overflowPoint = &(ptwXY->overflowPoints[ptwXY->overflowLength]);
913 if( legx == ptwXY_lessEqualGreaterX_lessThan ) {
914 overflowPoint->prior = greaterThanXPoint.prior;
915 overflowPoint->index = 0; }
917 if( legx == ptwXY_lessEqualGreaterX_greater ) {
918 overflowPoint->prior = overflowHeader->prior;
919 overflowPoint->index = ptwXY->length; }
921 overflowPoint->prior = lessThanEqualXPoint.prior;
922 if( lessThanEqualXPoint.next != NULL ) {
923 if( lessThanEqualXPoint.point.x < x )
924 overflowPoint->prior = lessThanEqualXPoint.prior->next;
927 for( p = overflowHeader->next, i = 1; p != overflowHeader; p = p->next, i++ )
928 if( p->point.x > x )
break;
930 overflowPoint->index = lessThanEqualXPoint.index + i;
933 overflowPoint->next = overflowPoint->prior->next;
934 overflowPoint->prior->next = overflowPoint;
935 overflowPoint->next->prior = overflowPoint;
936 point = &(overflowPoint->point);
937 for( overflowPoint = overflowPoint->next; overflowPoint != overflowHeader; overflowPoint = overflowPoint->next ) {
938 overflowPoint->index++;
940 ptwXY->overflowLength++;
944 case ptwXY_lessEqualGreaterX_empty :
945 point = ptwXY->points;
947 case ptwXY_lessEqualGreaterX_equal :
948 if( closeIsEqual && !
override )
return( status );
949 if( lessThanEqualXPoint.next == NULL ) {
950 point = &(ptwXY->points[lessThanEqualXPoint.index]); }
952 point = &(lessThanEqualXPoint.prior->next->point);
956 if( status == nfu_Okay ) {
959 if( legx != ptwXY_lessEqualGreaterX_equal ) ptwXY->length++;
976 double *xs, *p1, *p2;
979 if( length < 0 )
return( nfu_badInput );
980 if( length == 0 )
return( nfu_Okay );
981 if( ( xs = (
double *)
nfu_malloc( length *
sizeof(
double ) ) ) == NULL )
return( nfu_mallocError );
982 for( i = 0, p1 = xs, p2 = xys; i < length; i++, p1++, p2 += 2 ) *p1 = *p2;
991 static nfu_status
ptwXY_mergeFrom( ptwXYPoints *ptwXY,
int ,
int length,
double *xs,
double *ys ) {
994 double *sortedXs, *p1, *p2;
996 ptwXYPoint *point1, *point2;
998 if( length < 0 )
return( nfu_badInput );
999 if( length == 0 )
return( nfu_Okay );
1004 if( ( sortedXs = (
double *)
nfu_malloc( length *
sizeof(
double ) ) ) == NULL )
return( nfu_mallocError );
1006 for( i = 0, p1 = sortedXs, p2 = xs; i < length; i++, p1++, p2++ ) *p1 = *p2;
1011 for( i = 0, p1 = sortedXs, j = 0, n = 0; i < length; i++, p1++, n++ ) {
1012 for( ; j < ptwXY->length; j++, n++ ) {
1013 if( *p1 <= ptwXY->points[j].
x )
break;
1015 if( j == ptwXY->length )
break;
1017 n += (int) ( ( length - i ) + ( ptwXY->length - j ) );
1020 point1 = &(ptwXY->points[n-1]);
1021 point2 = &(ptwXY->points[length-1]);
1022 for( i = 0, j = 0, p1 = &(sortedXs[length-1]); ( i < length ) && ( j < length ) && ( n > 0 ); n--, point1-- ) {
1023 if( *p1 >= point2->x ) {
1025 point1->y = ys[(int)(p1 - xs)];
1026 if( *p1 >= point2->x ) {
1038 for( ; i < length; i++, p1--, point1-- ) {
1040 point1->y = ys[(int)(p1 - xs)];
1042 for( ; j < length; j++, point1--, point2-- ) *point1 = *point2;
1053 double d1 = *((
double *) x1p),
d2 = *((
double *) x2p);
1055 if( d1 <
d2 )
return( -1 );
1056 if( d1 ==
d2 )
return( 0 );
1065 ptwXY_dataFrom dataFrom;
1067 if( ptwXY->length != 0 ) {
1069 if( xMax >= x )
return( nfu_XNotAscending );
1072 if( nonOverflowLength < ptwXY->allocatedSize ) {
1073 ptwXY->points[nonOverflowLength].x =
x;
1074 ptwXY->points[nonOverflowLength].y = y; }
1076 if( ptwXY->overflowLength == ptwXY->overflowAllocatedSize ) {
1077 ptwXYPoint newPoint = {
x, y };
1078 return(
ptwXY_coalescePoints( ptwXY, ptwXY->length + ptwXY->overflowAllocatedSize, &newPoint, 0 ) ); }
1080 ptwXYOverflowPoint *overflowPoint = &(ptwXY->overflowPoints[ptwXY->overflowLength]);
1082 overflowPoint->prior = ptwXY->overflowHeader.prior;
1083 overflowPoint->next = overflowPoint->prior->next;
1084 overflowPoint->index = ptwXY->length;
1085 overflowPoint->prior->next = overflowPoint;
1086 overflowPoint->next->prior = overflowPoint;
1087 overflowPoint->point.x =
x;
1088 overflowPoint->point.y = y;
1089 ptwXY->overflowLength++;
1101 ptwXYOverflowPoint *overflowPoint, *pm1, *pp1;
1103 if( ptwXY->status != nfu_Okay )
return( ptwXY->status );
1105 if( ( index < 0 ) || ( index >= ptwXY->length ) )
return( nfu_badIndex );
1106 for( overflowPoint = ptwXY->overflowHeader.next, i = 0; overflowPoint != &(ptwXY->overflowHeader); overflowPoint = overflowPoint->next, i++ ) {
1107 if( overflowPoint->index >= index )
break;
1110 pm1 = pp1 = overflowPoint;
1111 if( overflowPoint->index == index ) {
1112 pp1 = overflowPoint->next;
1115 if( ( pp1 != &(ptwXY->overflowHeader) ) && ( pp1->index == ( index + 1 ) ) ) {
1116 if( pp1->point.x <= x )
return( nfu_badIndexForX ); }
1118 if( ( ( index + 1 ) < ptwXY->length ) && ( ptwXY->points[index + 1 - ip1].x <= x ) )
return( nfu_badIndexForX );
1120 if( overflowPoint != &(ptwXY->overflowHeader) ) pm1 = overflowPoint->prior;
1121 if( ( pm1 != &(ptwXY->overflowHeader) ) && ( pm1->index == ( index - 1 ) ) ) {
1122 if( pm1->point.x >= x )
return( nfu_badIndexForX ); }
1124 if( ( ( index - 1 ) >= 0 ) && ( ptwXY->points[index - 1 - i].x >= x ) )
return( nfu_badIndexForX );
1126 if( ( overflowPoint != &(ptwXY->overflowHeader) ) && ( overflowPoint->index == index ) ) {
1127 overflowPoint->point.x =
x;
1128 overflowPoint->point.y = y; }
1131 ptwXY->points[index].x =
x;
1132 ptwXY->points[index].y = y;
1141 nfu_status status = nfu_Okay;
1142 ptwXYOverflowPoint lessThanEqualXPoint, greaterThanXPoint;
1143 ptwXY_lessEqualGreaterX legx =
ptwXY_getPointsAroundX( ptwXY, x, &lessThanEqualXPoint, &greaterThanXPoint );
1147 if( ( side !=
'-' ) && ( side !=
'+' ) )
return( nfu_badInput );
1150 case ptwXY_lessEqualGreaterX_empty :
1151 case ptwXY_lessEqualGreaterX_lessThan :
1152 case ptwXY_lessEqualGreaterX_greater :
1153 status = nfu_XOutsideDomain;
1155 case ptwXY_lessEqualGreaterX_between :
1156 *slope = ( greaterThanXPoint.point.y - lessThanEqualXPoint.point.y ) /
1157 ( greaterThanXPoint.point.x - lessThanEqualXPoint.point.x );
1159 case ptwXY_lessEqualGreaterX_equal :
1161 if( lessThanEqualXPoint.index == 0 ) {
1162 status = nfu_XOutsideDomain; }
1165 *slope = ( lessThanEqualXPoint.point.y - point->y ) / ( lessThanEqualXPoint.point.x - point->x );
1168 if( lessThanEqualXPoint.index == ( ptwXY->length - 1 ) ) {
1169 status = nfu_XOutsideDomain; }
1172 *slope = ( point->y - lessThanEqualXPoint.point.y ) / ( point->x - lessThanEqualXPoint.point.x );
1187 *dataFrom = ptwXY_dataFrom_Unknown;
1188 if( ptwXY->overflowLength > 0 ) {
1189 *dataFrom = ptwXY_dataFrom_Overflow;
1190 xMin = ptwXY->overflowHeader.next->point.x;
1191 if( nonOverflowLength >= 0 ) {
1192 if( xMin > ptwXY->points[0].x ) {
1193 *dataFrom = ptwXY_dataFrom_Points;
1194 xMin = ptwXY->points[0].x;
1197 else if( nonOverflowLength > 0 ) {
1198 *dataFrom = ptwXY_dataFrom_Points;
1199 xMin = ptwXY->points[0].x;
1208 ptwXY_dataFrom dataFrom;
1220 *dataFrom = ptwXY_dataFrom_Unknown;
1221 if( ptwXY->overflowLength > 0 ) {
1222 *dataFrom = ptwXY_dataFrom_Overflow;
1223 xMax = ptwXY->overflowHeader.prior->point.x;
1224 if( ( nonOverflowLength > 0 ) ) {
1225 if( xMax < ptwXY->points[nonOverflowLength-1].
x ) {
1226 *dataFrom = ptwXY_dataFrom_Points;
1227 xMax = ptwXY->points[nonOverflowLength-1].x;
1230 else if( ptwXY->length > 0 ) {
1231 *dataFrom = ptwXY_dataFrom_Points;
1232 xMax = ptwXY->points[nonOverflowLength-1].x;
1241 ptwXY_dataFrom dataFrom;
1251 ptwXYPoint *p = ptwXY->points;
1252 ptwXYOverflowPoint *overflowPoint = ptwXY->overflowHeader.next;
1255 if( ptwXY->length == 0 )
return( 0. );
1258 for( i = 1, p++; i <
n; i++, p++ ) yMin = ( ( yMin < p->y ) ? yMin : p->y ); }
1260 yMin = overflowPoint->point.y;
1262 for( ; overflowPoint != &(ptwXY->overflowHeader); overflowPoint = overflowPoint->next )
1263 yMin = ( ( yMin < overflowPoint->point.y ) ? yMin : overflowPoint->point.y );
1272 ptwXYPoint *p = ptwXY->points;
1273 ptwXYOverflowPoint *overflowPoint = ptwXY->overflowHeader.next;
1276 if( ptwXY->length == 0 )
return( 0. );
1279 for( i = 1, p++; i < n; i++, p++ ) yMax = ( ( yMax > p->y ) ? yMax : p->y ); }
1281 yMax = overflowPoint->point.y;
1283 for( ; overflowPoint != &(ptwXY->overflowHeader); overflowPoint = overflowPoint->next )
1284 yMax = ( ( yMax > overflowPoint->point.y ) ? yMax : overflowPoint->point.y );
1292 overflowPoint->prior = prior;
1293 overflowPoint->next = next;
1294 overflowPoint->index = -1;
1295 overflowPoint->point.x = 0.;
1296 overflowPoint->point.y = 0.;
1299 #if defined __cplusplus
void * nfu_realloc(size_t size, void *old)
static char const logLinInterpolationString[]
ptwXY_lessEqualGreaterX ptwXY_getPointsAroundX(ptwXYPoints *ptwXY, double x, ptwXYOverflowPoint *lessThanEqualXPoint, ptwXYOverflowPoint *greaterThanXPoint)
nfu_status ptwXY_setXYDataFromXsAndYs(ptwXYPoints *ptwXY, int64_t length, double const *x, double const *y)
double ptwXY_getXMinAndFrom(ptwXYPoints *ptwXY, ptwXY_dataFrom *dataFrom)
nfu_status ptwXY_release(ptwXYPoints *ptwXY)
ptwXY_interpolation ptwXY_getInterpolation(ptwXYPoints *ptwXY)
ptwXYPoints * ptwXY_xMaxSlice(ptwXYPoints *ptwXY, double xMax, int64_t secondarySize, int fill, nfu_status *status)
int64_t ptwXY_length(ptwXYPoints *ptwXY)
ptwXYPoints * ptwXY_create(ptwXY_interpolation interpolation, ptwXY_interpolationOtherInfo const *interpolationOtherInfo, double biSectionMax, double accuracy, int64_t primarySize, int64_t secondarySize, int64_t length, double const *xy, nfu_status *status, int userFlag)
int64_t ptwXY_getNonOverflowLength(ptwXYPoints const *ptwXY)
static char const logLogInterpolationString[]
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
double ptwXY_getXMaxAndFrom(ptwXYPoints *ptwXY, ptwXY_dataFrom *dataFrom)
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
static const G4double eps
nfu_status ptwXY_getSlopeAtX(ptwXYPoints *ptwXY, double x, const char side, double *slope)
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
void ptwXY_setUserFlag(ptwXYPoints *ptwXY, int userFlag)
static void ptwXY_initialOverflowPoint(ptwXYOverflowPoint *overflowPoint, ptwXYOverflowPoint *prior, ptwXYOverflowPoint *next)
double ptwXY_getBiSectionMax(ptwXYPoints *ptwXY)
char const * ptwXY_getInterpolationString(ptwXYPoints *ptwXY)
double ptwXY_getXMax(ptwXYPoints *ptwXY)
nfu_status ptwXY_mergeFromXYs(ptwXYPoints *ptwXY, int length, double *xys)
ptwXYPoints * ptwXY_slice(ptwXYPoints *ptwXY, int64_t index1, int64_t index2, int64_t secondarySize, nfu_status *status)
nfu_status ptwXY_setXYData(ptwXYPoints *ptwXY, int64_t length, double const *xy)
double ptwXY_getAccuracy(ptwXYPoints *ptwXY)
nfu_status ptwXY_setValueAtX_overrideIfClose(ptwXYPoints *ptwXY, double x, double y, double eps, int override)
ptwXY_lessEqualGreaterX ptwXY_getPointsAroundX_closeIsEqual(ptwXYPoints *ptwXY, double x, ptwXYOverflowPoint *lessThanEqualXPoint, ptwXYOverflowPoint *greaterThanXPoint, double eps, int *closeIsEqual, ptwXYPoint **closePoint)
nfu_status ptwXY_copy(ptwXYPoints *dest, ptwXYPoints *src)
static char const linLinInterpolationString[]
double ptwXY_getYMax(ptwXYPoints *ptwXY)
nfu_status ptwXY_mergeFromXsAndYs(ptwXYPoints *ptwXY, int length, double *xs, double *ys)
nfu_status ptwXY_coalescePoints(ptwXYPoints *ptwXY, int64_t size, ptwXYPoint *newPoint, int forceSmallerResize)
ptwXYPoints * ptwXY_xSlice(ptwXYPoints *ptwXY, double xMin, double xMax, int64_t secondarySize, int fill, nfu_status *status)
double ptwXY_getYMin(ptwXYPoints *ptwXY)
nfu_status ptwXY_getXYPairAtIndex(ptwXYPoints *ptwXY, int64_t index, double *x, double *y)
nfu_status ptwXY_reallocatePoints(ptwXYPoints *ptwXY, int64_t size, int forceSmallerResize)
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_mergeFrom(ptwXYPoints *ptwXY, int incY, int length, double *xs, double *ys)
int ptwXY_getUserFlag(ptwXYPoints *ptwXY)
nfu_status ptwXY_getValueAtX(ptwXYPoints *ptwXY, double x, double *y)
const G4double x[NPOINTSGL]
ptwXYPoints * ptwXY_createFrom_Xs_Ys(ptwXY_interpolation interpolation, ptwXY_interpolationOtherInfo const *interpolationOtherInfo, double biSectionMax, double accuracy, int64_t primarySize, int64_t secondarySize, int64_t length, double const *Xs, double const *Ys, nfu_status *status, int userFlag)
ptwXYPoint * ptwXY_getPointAtIndex(ptwXYPoints *ptwXY, int64_t index)
double ptwXY_getXMin(ptwXYPoints *ptwXY)
nfu_status ptwXY_reallocateOverflowPoints(ptwXYPoints *ptwXY, int64_t size)
ptwXYPoints * ptwXY_cloneToInterpolation(ptwXYPoints *ptwXY, ptwXY_interpolation interpolationTo, nfu_status *status)
double ptwXY_setAccuracy(ptwXYPoints *ptwXY, double accuracy)
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_appendXY(ptwXYPoints *ptwXY, double x, double y)
double ptwXY_setBiSectionMax(ptwXYPoints *ptwXY, double biSectionMax)
ptwXYPoints * ptwXY_xMinSlice(ptwXYPoints *ptwXY, double xMin, int64_t secondarySize, int fill, nfu_status *status)
static char const flatInterpolationString[]
nfu_status ptwXY_deletePoints(ptwXYPoints *ptwXY, int64_t i1, int64_t i2)
nfu_status ptwXY_interpolatePoint(ptwXY_interpolation interpolation, double x, double *y, double x1, double y1, double x2, double y2)
nfu_status ptwXY_setup(ptwXYPoints *ptwXY, ptwXY_interpolation interpolation, ptwXY_interpolationOtherInfo const *interpolationOtherInfo, double biSectionMax, double accuracy, int64_t primarySize, int64_t secondarySize, int userFlag)
nfu_status ptwXY_setXYPairAtIndex(ptwXYPoints *ptwXY, int64_t index, double x, double y)
void * nfu_calloc(size_t size, size_t n)
nfu_status ptwXY_getStatus(ptwXYPoints *ptwXY)
void * nfu_malloc(size_t size)
static char const linLogInterpolationString[]
static int ptwXY_mergeCompareFunction(void const *x1p, void const *x2p)