| Geant4 Cross Reference |
1 /*
2 # <<BEGIN-copyright>>
3 # <<END-copyright>>
4 */
5
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <cmath>
9
10 #include "ptwXY.h"
11
12 #if defined __cplusplus
13 namespace GIDI {
14 using namespace GIDI;
15 #endif
16
17 static char const linLinInterpolationString[] = "linear,linear";
18 static char const linLogInterpolationString[] = "linear,log";
19 static char const logLinInterpolationString[] = "log,linear";
20 static char const logLogInterpolationString[] = "log,log";
21 static char const flatInterpolationString[] = "flat";
22
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 );
25 static int ptwXY_mergeCompareFunction( void const *x1p, void const *x2p );
26 /*
27 ************************************************************
28 */
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 ) {
31
32 ptwXYPoints *ptwXY = (ptwXYPoints *) nfu_calloc( sizeof( ptwXYPoints ), 1 );
33
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 );
40 }
41 return( ptwXY );
42 }
43 /*
44 ************************************************************
45 */
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 ) {
48
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 :
58 ptwXY->interpolationOtherInfo.interpolationString = linLinInterpolationString; break;
59 case ptwXY_interpolationLinLog :
60 ptwXY->interpolationOtherInfo.interpolationString = linLogInterpolationString; break;
61 case ptwXY_interpolationLogLin :
62 ptwXY->interpolationOtherInfo.interpolationString = logLinInterpolationString; break;
63 case ptwXY_interpolationLogLog :
64 ptwXY->interpolationOtherInfo.interpolationString = logLogInterpolationString; break;
65 case ptwXY_interpolationFlat :
66 ptwXY->interpolationOtherInfo.interpolationString = flatInterpolationString; break;
67 case ptwXY_interpolationOther : /* For ptwXY_interpolationOther, interpolationOtherInfo and interpolationString must be defined. */
68 if( interpolationOtherInfo == NULL ) {
69 ptwXY->status = nfu_otherInterpolation; }
70 else {
71 if( interpolationOtherInfo->interpolationString == NULL ) {
72 ptwXY->status = nfu_otherInterpolation; }
73 else {
74 if( ( ptwXY->interpolationOtherInfo.interpolationString = strdup( interpolationOtherInfo->interpolationString ) ) == NULL ) {
75 ptwXY->status = nfu_mallocError;
76 }
77 }
78 ptwXY->interpolationOtherInfo.getValueFunc = interpolationOtherInfo->getValueFunc;
79 ptwXY->interpolationOtherInfo.argList = interpolationOtherInfo->argList;
80 }
81 }
82 ptwXY->userFlag = 0;
83 ptwXY_setUserFlag( ptwXY, userFlag );
84 ptwXY->biSectionMax = ptwXY_maxBiSectionMax;
85 ptwXY_setBiSectionMax( ptwXY, biSectionMax );
86 ptwXY->accuracy = ptwXY_minAccuracy;
87 ptwXY_setAccuracy( ptwXY, accuracy );
88
89 ptwXY->length = 0;
90 ptwXY->allocatedSize = 0;
91 ptwXY->overflowLength = 0;
92 ptwXY->overflowAllocatedSize = 0;
93 ptwXY->mallocFailedSize = 0;
94
95 ptwXY_initialOverflowPoint( &(ptwXY->overflowHeader), &(ptwXY->overflowHeader), &(ptwXY->overflowHeader) );
96
97 ptwXY->points = NULL;
98 ptwXY->overflowPoints = NULL;
99
100 ptwXY_reallocatePoints( ptwXY, primarySize, 0 );
101 ptwXY_reallocateOverflowPoints( ptwXY, secondarySize );
102 if( ptwXY->status != nfu_Okay ) ptwXY_release( ptwXY );
103 return( ptwXY->status );
104 }
105 /*
106 ************************************************************
107 */
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 ) {
111
112 ptwXYPoints *ptwXY;
113
114 if( primarySize < length ) primarySize = length;
115 if( ( ptwXY = ptwXY_new( interpolation, interpolationOtherInfo, biSectionMax, accuracy, primarySize,
116 secondarySize, status, userFlag ) ) != NULL ) {
117 if( ( *status = ptwXY_setXYData( ptwXY, length, xy ) ) != nfu_Okay ) {
118 ptwXY = ptwXY_free( ptwXY );
119 }
120 }
121 return( ptwXY );
122 }
123 /*
124 ************************************************************
125 */
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 ) {
129
130 int i;
131 ptwXYPoints *ptwXY;
132
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];
139 }
140 ptwXY->length = length;
141 }
142
143 return( ptwXY );
144 }
145 /*
146 ************************************************************
147 */
148 nfu_status ptwXY_copy( ptwXYPoints *dest, ptwXYPoints *src ) {
149
150 int64_t i, nonOverflowLength = ptwXY_getNonOverflowLength( src );
151 ptwXYPoint *pointFrom, *pointTo;
152 ptwXYOverflowPoint *o, *overflowHeader = &(src->overflowHeader);
153
154 if( dest->status != nfu_Okay ) return( dest->status );
155 if( src->status != nfu_Okay ) return( src->status );
156
157 ptwXY_clear( dest );
158 if( dest->interpolation == ptwXY_interpolationOther ) {
159 if( dest->interpolationOtherInfo.interpolationString != NULL ) {
160 dest->interpolationOtherInfo.interpolationString = (char const *) nfu_free( (void *) dest->interpolationOtherInfo.interpolationString );
161 }
162 }
163 dest->interpolation = ptwXY_interpolationLinLin; /* This and prior lines are in case interpolation is 'other' and ptwXY_reallocatePoints fails. */
164 if( dest->allocatedSize < src->length ) ptwXY_reallocatePoints( dest, src->length, 0 );
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 );
171 } }
172 else {
173 dest->interpolationOtherInfo.interpolationString = src->interpolationOtherInfo.interpolationString;
174 }
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;
184 i = 0;
185 while( o != overflowHeader ) {
186 if( i < nonOverflowLength ) {
187 if( pointFrom->x < o->point.x ) {
188 *pointTo = *pointFrom;
189 i++;
190 pointFrom++; }
191 else {
192 *pointTo = o->point;
193 o = o->next;
194 } }
195 else {
196 *pointTo = o->point;
197 o = o->next;
198 }
199 pointTo++;
200 } // Loop checking, 11.06.2015, T. Koi
201 for( ; i < nonOverflowLength; i++, pointFrom++, pointTo++ ) *pointTo = *pointFrom;
202 dest->length = src->length;
203 return( dest->status );
204 }
205 /*
206 ************************************************************
207 */
208 ptwXYPoints *ptwXY_clone( ptwXYPoints *ptwXY, nfu_status *status ) {
209
210 return( ptwXY_slice( ptwXY, 0, ptwXY->length, ptwXY->overflowAllocatedSize, status ) );
211 }
212 /*
213 ************************************************************
214 */
215 ptwXYPoints *ptwXY_cloneToInterpolation( ptwXYPoints *ptwXY, ptwXY_interpolation interpolationTo, nfu_status *status ) {
216
217 ptwXYPoints *n1;
218
219 if( interpolationTo == ptwXY_interpolationOther ) {
220 *status = nfu_otherInterpolation;
221 return( NULL );
222 }
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 :
228 n1->interpolationOtherInfo.interpolationString = linLinInterpolationString; break;
229 case ptwXY_interpolationLinLog :
230 n1->interpolationOtherInfo.interpolationString = linLogInterpolationString; break;
231 case ptwXY_interpolationLogLin :
232 n1->interpolationOtherInfo.interpolationString = logLinInterpolationString; break;
233 case ptwXY_interpolationLogLog :
234 n1->interpolationOtherInfo.interpolationString = logLogInterpolationString; break;
235 case ptwXY_interpolationFlat :
236 n1->interpolationOtherInfo.interpolationString = flatInterpolationString; break;
237 case ptwXY_interpolationOther : /* Does not happen, but needed to stop compilers from complaining. */
238 break;
239 }
240 n1->interpolationOtherInfo.getValueFunc = NULL;
241 n1->interpolationOtherInfo.argList = NULL;
242 }
243 return( n1 );
244 }
245 /*
246 ************************************************************
247 */
248 ptwXYPoints *ptwXY_slice( ptwXYPoints *ptwXY, int64_t index1, int64_t index2, int64_t secondarySize, nfu_status *status ) {
249
250 int64_t i, length;
251 ptwXYPoints *n;
252
253 *status = nfu_badSelf;
254 if( ptwXY->status != nfu_Okay ) return( NULL );
255
256 *status = nfu_badIndex;
257 if( index2 < index1 ) return( NULL );
258 if( index1 < 0 ) index1 = 0;
259 if( index2 > ptwXY->length ) index2 = ptwXY->length;
260
261 length = index2 - index1;
262 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( NULL );
263 if( ( n = ptwXY_new( ptwXY->interpolation, &(ptwXY->interpolationOtherInfo), ptwXY->biSectionMax,
264 ptwXY->accuracy, length, secondarySize, status, ptwXY->userFlag ) ) == NULL ) return( NULL );
265
266 *status = n->status = ptwXY->status;
267 for( i = index1; i < index2; i++ ) n->points[i - index1] = ptwXY->points[i];
268 n->length = length;
269 return( n );
270 }
271 /*
272 ************************************************************
273 */
274 ptwXYPoints *ptwXY_xSlice( ptwXYPoints *ptwXY, double xMin, double xMax, int64_t secondarySize, int fill, nfu_status *status ) {
275
276 int64_t i, i1, i2;
277 double y;
278 ptwXYPoints *n = NULL;
279
280 if( ( *status = ptwXY->status ) != nfu_Okay ) return( NULL );
281
282 if( ( ptwXY->length == 0 ) || ( ptwXY_getXMin( ptwXY ) >= xMax ) || ( ptwXY_getXMax( ptwXY ) <= xMin ) ) {
283 n = ptwXY_new( ptwXY->interpolation, &(ptwXY->interpolationOtherInfo), ptwXY->biSectionMax,
284 ptwXY->accuracy, 0, secondarySize, status, ptwXY->userFlag ); }
285 else {
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 ) ) {
289 if( ( *status = ptwXY_getValueAtX( n, xMax, &y ) ) != nfu_Okay ) goto Err;
290 if( ( *status = ptwXY_setValueAtX( n, xMax, y ) ) != nfu_Okay ) goto Err;
291 }
292 if( fill && ( n->points[0].x < xMin ) ) {
293 if( ( *status = ptwXY_getValueAtX( n, xMin, &y ) ) != nfu_Okay ) goto Err;
294 if( ( *status = ptwXY_setValueAtX( n, xMin, y ) ) != nfu_Okay ) goto Err;
295 }
296 ptwXY_coalescePoints( n, n->length + n->overflowAllocatedSize, NULL, 0 );
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;
299 i2++;
300 if( i1 > 0 ) {
301 for( i = i1; i < i2; i++ ) n->points[i- i1] = n->points[i];
302 }
303 n->length = i2 - i1;
304 }
305 }
306 return( n );
307
308 Err:
309 if( n != NULL ) ptwXY_free( n );
310 return( NULL );
311 }
312 /*
313 ************************************************************
314 */
315 ptwXYPoints *ptwXY_xMinSlice( ptwXYPoints *ptwXY, double xMin, int64_t secondarySize, int fill, nfu_status *status ) {
316
317 double xMax = 1.1 * xMin + 1;
318
319 if( xMin < 0 ) xMax = 0.9 * xMin + 1;
320 if( ptwXY->length > 0 ) xMax = ptwXY_getXMax( ptwXY );
321 return( ptwXY_xSlice( ptwXY, xMin, xMax, secondarySize, fill, status ) );
322 }
323 /*
324 ************************************************************
325 */
326 ptwXYPoints *ptwXY_xMaxSlice( ptwXYPoints *ptwXY, double xMax, int64_t secondarySize, int fill, nfu_status *status ) {
327
328 double xMin = 0.9 * xMax - 1;
329
330 if( xMax < 0 ) xMin = 1.1 * xMax - 1;
331 if( ptwXY->length > 0 ) xMin = ptwXY_getXMin( ptwXY );
332 return( ptwXY_xSlice( ptwXY, xMin, xMax, secondarySize, fill, status ) );
333 }
334 /*
335 ************************************************************
336 */
337 ptwXY_interpolation ptwXY_getInterpolation( ptwXYPoints *ptwXY ) {
338
339 return( ptwXY->interpolation );
340 }
341 /*
342 ************************************************************
343 */
344 char const *ptwXY_getInterpolationString( ptwXYPoints *ptwXY ) {
345
346 return( ptwXY->interpolationOtherInfo.interpolationString );
347 }
348 /*
349 ************************************************************
350 */
351 nfu_status ptwXY_getStatus( ptwXYPoints *ptwXY ) {
352
353 return( ptwXY->status );
354 }
355 /*
356 ************************************************************
357 */
358 int ptwXY_getUserFlag( ptwXYPoints *ptwXY ) {
359
360 return( ptwXY->userFlag );
361 }
362 /*
363 ************************************************************
364 */
365 void ptwXY_setUserFlag( ptwXYPoints *ptwXY, int userFlag ) {
366
367 ptwXY->userFlag = userFlag;
368 }
369 /*
370 ************************************************************
371 */
372 double ptwXY_getAccuracy( ptwXYPoints *ptwXY ) {
373
374 return( ptwXY->accuracy );
375 }
376 /*
377 ************************************************************
378 */
379 double ptwXY_setAccuracy( ptwXYPoints *ptwXY, double accuracy ) {
380
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 );
386 }
387 /*
388 ************************************************************
389 */
390 double ptwXY_getBiSectionMax( ptwXYPoints *ptwXY ) {
391
392 return( ptwXY->biSectionMax );
393 }
394 /*
395 ************************************************************
396 */
397 double ptwXY_setBiSectionMax( ptwXYPoints *ptwXY, double biSectionMax ) {
398
399 if( biSectionMax < 0 ) {
400 biSectionMax = 0; }
401 else if( biSectionMax > ptwXY_maxBiSectionMax ) {
402 biSectionMax = ptwXY_maxBiSectionMax;
403 }
404 ptwXY->biSectionMax = biSectionMax;
405 return( ptwXY->biSectionMax );
406 }
407 /*
408 ************************************************************
409 */
410 nfu_status ptwXY_reallocatePoints( ptwXYPoints *ptwXY, int64_t size, int forceSmallerResize ) {
411 /*
412 * This is for allocating/reallocating the primary data memory.
413 */
414 if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
415
416 if( size < ptwXY_minimumSize ) size = ptwXY_minimumSize; /* ptwXY_minimumSize must be > 0. */
417 if( size < ptwXY->length ) size = ptwXY->length;
418 if( size != ptwXY->allocatedSize ) {
419 if( size > ptwXY->allocatedSize ) { /* Increase size of allocated points. */
420 ptwXY->points = (ptwXYPoint *) nfu_realloc( (size_t) size * sizeof( ptwXYPoint ), ptwXY->points ); }
421 else if( ( ptwXY->allocatedSize > 2 * size ) || forceSmallerResize ) { /* Decrease size, if at least 1/2 size reduction or if forced to. */
422 ptwXY->points = (ptwXYPoint *) nfu_realloc( (size_t) size * sizeof( ptwXYPoint ), ptwXY->points ); }
423 else {
424 size = ptwXY->allocatedSize; /* Size is < ptwXY->allocatedSize, but realloc not called. */
425 }
426 if( ptwXY->points == NULL ) {
427 ptwXY->length = 0;
428 ptwXY->mallocFailedSize = size;
429 size = 0;
430 ptwXY->status = nfu_mallocError;
431 }
432 ptwXY->allocatedSize = size;
433 }
434 return( ptwXY->status );
435 }
436 /*
437 ************************************************************
438 */
439 nfu_status ptwXY_reallocateOverflowPoints( ptwXYPoints *ptwXY, int64_t size ) {
440 /*
441 * This is for allocating/reallocating the secondary data memory.
442 */
443 nfu_status status = nfu_Okay;
444
445 if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
446
447 if( size < ptwXY_minimumOverflowSize ) size = ptwXY_minimumOverflowSize; /* ptwXY_minimumOverflowSize must be > 0. */
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 ) {
453 ptwXY->length = 0;
454 ptwXY->overflowLength = 0;
455 ptwXY->mallocFailedSize = size;
456 size = 0;
457 ptwXY->status = nfu_mallocError;
458 }
459 }
460 ptwXY->overflowAllocatedSize = size; }
461 else {
462 ptwXY->status = status;
463 }
464 return( ptwXY->status );
465 }
466 /*
467 ************************************************************
468 */
469 nfu_status ptwXY_coalescePoints( ptwXYPoints *ptwXY, int64_t size, ptwXYPoint *newPoint, int forceSmallerResize ) {
470
471 int addNewPoint;
472 int64_t length = ptwXY->length + ( ( newPoint != NULL ) ? 1 : 0 );
473 ptwXYOverflowPoint *last = ptwXY->overflowHeader.prior;
474 ptwXYPoint *pointsFrom, *pointsTo;
475
476 if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
477 if( ptwXY->overflowLength == 0 ) return( nfu_Okay );
478
479 if( size < length ) size = length;
480 if( size > ptwXY->allocatedSize ) {
481 if( ptwXY_reallocatePoints( ptwXY, size, forceSmallerResize ) != nfu_Okay ) return( ptwXY->status );
482 }
483 pointsFrom = &(ptwXY->points[ptwXY_getNonOverflowLength( ptwXY ) - 1]);
484 pointsTo = &(ptwXY->points[length - 1]);
485 while( last != &(ptwXY->overflowHeader) ) {
486 addNewPoint = 0;
487 if( newPoint != NULL ) {
488 if( ( pointsFrom >= ptwXY->points ) && ( pointsFrom->x > last->point.x ) ) {
489 if( newPoint->x > pointsFrom->x ) addNewPoint = 1; }
490 else {
491 if( newPoint->x > last->point.x ) addNewPoint = 1;
492 }
493 if( addNewPoint == 1 ) {
494 *pointsTo = *newPoint;
495 newPoint = NULL;
496 }
497 }
498 if( addNewPoint == 0 ) {
499 if( ( pointsFrom >= ptwXY->points ) && ( pointsFrom->x > last->point.x ) ) {
500 *pointsTo = *pointsFrom;
501 pointsFrom--; }
502 else {
503 *pointsTo = last->point;
504 last = last->prior;
505 }
506 }
507 pointsTo--;
508 } // Loop checking, 11.06.2015, T. Koi
509 while( ( newPoint != NULL ) && ( pointsFrom >= ptwXY->points ) ) {
510 if( newPoint->x > pointsFrom->x ) {
511 *pointsTo = *newPoint;
512 newPoint = NULL; }
513 else {
514 *pointsTo = *pointsFrom;
515 pointsFrom--;
516 }
517 pointsTo--;
518 } // Loop checking, 11.06.2015, T. Koi
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;
524 return( nfu_Okay );
525 }
526 /*
527 ************************************************************
528 */
529 nfu_status ptwXY_simpleCoalescePoints( ptwXYPoints *ptwXY ) {
530
531 return( ptwXY_coalescePoints( ptwXY, ptwXY->length, NULL, 0 ) );
532 }
533 /*
534 ************************************************************
535 */
536 nfu_status ptwXY_clear( ptwXYPoints *ptwXY ) {
537
538 if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
539
540 ptwXY->length = 0;
541 ptwXY->overflowLength = 0;
542 ptwXY->overflowHeader.prior = &(ptwXY->overflowHeader);
543 ptwXY->overflowHeader.next = &(ptwXY->overflowHeader);
544 return( nfu_Okay );
545 }
546 /*
547 ************************************************************
548 */
549 nfu_status ptwXY_release( ptwXYPoints *ptwXY ) {
550 /*
551 * Note, this routine does not free ptwXY (i.e., it does not undo all of ptwXY_new).
552 */
553
554 if( ptwXY->interpolation == ptwXY_interpolationOther ) {
555 if( ptwXY->interpolationOtherInfo.interpolationString != NULL )
556 ptwXY->interpolationOtherInfo.interpolationString = (char const *) nfu_free( (void *) ptwXY->interpolationOtherInfo.interpolationString );
557 }
558 ptwXY->interpolation = ptwXY_interpolationLinLin;
559 ptwXY->interpolationOtherInfo.getValueFunc = NULL;
560 ptwXY->interpolationOtherInfo.argList = NULL;
561 ptwXY->length = 0;
562 ptwXY->allocatedSize = 0;
563 ptwXY->points = (ptwXYPoint *) nfu_free( ptwXY->points );
564
565 ptwXY->overflowLength = 0;
566 ptwXY->overflowAllocatedSize = 0;
567 ptwXY->overflowPoints = (ptwXYOverflowPoint *) nfu_free( ptwXY->overflowPoints );
568
569 return( nfu_Okay );
570 }
571 /*
572 ************************************************************
573 */
574 ptwXYPoints *ptwXY_free( ptwXYPoints *ptwXY ) {
575
576 if( ptwXY != NULL ) ptwXY_release( ptwXY );
577 nfu_free( ptwXY );
578 return( (ptwXYPoints *) NULL );
579 }
580 /*
581 ************************************************************
582 */
583 int64_t ptwXY_length( ptwXYPoints *ptwXY ) {
584
585 return( ptwXY->length );
586 }
587 /*
588 ************************************************************
589 */
590 int64_t ptwXY_getNonOverflowLength( ptwXYPoints const *ptwXY ) {
591
592 return( ptwXY->length - ptwXY->overflowLength );
593 }
594 /*
595 ************************************************************
596 */
597 nfu_status ptwXY_setXYData( ptwXYPoints *ptwXY, int64_t length, double const *xy ) {
598
599 nfu_status status = nfu_Okay;
600 int64_t i;
601 ptwXYPoint *p;
602 double const *d = xy;
603 double xOld = 0.;
604
605 if( length > ptwXY->allocatedSize ) {
606 status = ptwXY_reallocatePoints( ptwXY, length, 0 );
607 if( status != nfu_Okay ) return( status );
608 }
609 for( i = 0, p = ptwXY->points; i < length; i++, p++ ) {
610 if( i != 0 ) {
611 if( *d <= xOld ) {
612 status = nfu_XNotAscending;
613 ptwXY->status = nfu_XNotAscending;
614 length = 0;
615 break;
616 }
617 }
618 xOld = *d;
619 p->x = *(d++);
620 p->y = *(d++);
621 }
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 );
627 }
628 /*
629 ************************************************************
630 */
631 nfu_status ptwXY_setXYDataFromXsAndYs( ptwXYPoints *ptwXY, int64_t length, double const *x, double const *y ) {
632
633 nfu_status status;
634 int64_t i;
635 ptwXYPoint *p;
636 double xOld = 0.;
637
638 if( ( status = ptwXY_clear( ptwXY ) ) != nfu_Okay ) return( status );
639 if( length > ptwXY->allocatedSize ) {
640 if( ( status = ptwXY_reallocatePoints( ptwXY, length, 0 ) ) != nfu_Okay ) return( status );
641 }
642 for( i = 0, p = ptwXY->points; i < length; i++, p++, x++, y++ ) {
643 if( i != 0 ) {
644 if( *x <= xOld ) {
645 status = ptwXY->status = nfu_XNotAscending;
646 length = 0;
647 break;
648 }
649 }
650 xOld = *x;
651 p->x = *x;
652 p->y = *y;
653 }
654 ptwXY->length = length;
655 return( status );
656 }
657 /*
658 ************************************************************
659 */
660 nfu_status ptwXY_deletePoints( ptwXYPoints *ptwXY, int64_t i1, int64_t i2 ) {
661
662 int64_t n = ptwXY->length - ( i2 - i1 );
663
664 if( ( ptwXY->status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( ptwXY->status );
665 if( ( i1 < 0 ) || ( i1 > i2 ) || ( i2 > ptwXY->length ) ) return( nfu_badIndex );
666 if( i1 != i2 ) {
667 for( ; i2 < ptwXY->length; i1++, i2++ ) ptwXY->points[i1] = ptwXY->points[i2];
668 ptwXY->length = n;
669 }
670 return( ptwXY->status );
671 }
672 /*
673 ************************************************************
674 */
675 ptwXYPoint *ptwXY_getPointAtIndex( ptwXYPoints *ptwXY, int64_t index ) {
676
677 if( ptwXY->status != nfu_Okay ) return( NULL );
678 if( ( index < 0 ) || ( index >= ptwXY->length ) ) return( NULL );
679 return( ptwXY_getPointAtIndex_Unsafely( ptwXY, index ) );
680 }
681 /*
682 ************************************************************
683 */
684 ptwXYPoint *ptwXY_getPointAtIndex_Unsafely( ptwXYPoints *ptwXY, int64_t index ) {
685
686 int64_t i;
687 ptwXYOverflowPoint *overflowPoint;
688
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;
692 }
693 return( &(ptwXY->points[index - i]) );
694 }
695 /*
696 ************************************************************
697 */
698 nfu_status ptwXY_getXYPairAtIndex( ptwXYPoints *ptwXY, int64_t index, double *x, double *y ) {
699
700 ptwXYPoint *p = ptwXY_getPointAtIndex( ptwXY, index );
701
702 if( p == NULL ) return( nfu_badIndex );
703 *x = p->x;
704 *y = p->y;
705 return( nfu_Okay );
706 }
707 /*
708 ************************************************************
709 */
710 ptwXY_lessEqualGreaterX ptwXY_getPointsAroundX( ptwXYPoints *ptwXY, double x, ptwXYOverflowPoint *lessThanEqualXPoint, ptwXYOverflowPoint *greaterThanXPoint ) {
711
712 int closeIsEqual;
713 ptwXYPoint *closePoint;
714
715 return( ptwXY_getPointsAroundX_closeIsEqual( ptwXY, x, lessThanEqualXPoint, greaterThanXPoint, 0, &closeIsEqual, &closePoint ) );
716 }
717 /*
718 ************************************************************
719 */
720 ptwXY_lessEqualGreaterX ptwXY_getPointsAroundX_closeIsEqual( ptwXYPoints *ptwXY, double x, ptwXYOverflowPoint *lessThanEqualXPoint,
721 ptwXYOverflowPoint *greaterThanXPoint, double eps, int *closeIsEqual, ptwXYPoint **closePoint ) {
722
723 int64_t overflowIndex, nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
724 int64_t indexMin, indexMid, indexMax;
725 ptwXY_dataFrom xMinFrom, xMaxFrom;
726 double xMin = ptwXY_getXMinAndFrom( ptwXY, &xMinFrom ), xMax = ptwXY_getXMaxAndFrom( ptwXY, &xMaxFrom );
727 ptwXYOverflowPoint *overflowPoint, *overflowHeader = &(ptwXY->overflowHeader);
728 ptwXY_lessEqualGreaterX status = ptwXY_lessEqualGreaterX_empty;
729 ptwXYPoint *lowerPoint = NULL, *upperPoint = NULL;
730
731 ptwXY_initialOverflowPoint( lessThanEqualXPoint, overflowHeader, NULL );
732 ptwXY_initialOverflowPoint( greaterThanXPoint, overflowHeader, NULL );
733 if( ptwXY->length != 0 ) {
734 if( x < xMin ) {
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]); }
741 else {
742 *greaterThanXPoint = *(overflowHeader->next);
743 *closePoint = &(overflowHeader->next->point);
744 } }
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]); }
752 else {
753 *lessThanEqualXPoint = *(overflowHeader->prior);
754 *closePoint = &(overflowHeader->prior->point);
755 } }
756 else { /* xMin <= x <= xMax */
757 status = ptwXY_lessEqualGreaterX_between; /* Default for this condition, can only be between or equal. */
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 ) { /* If here and length = 1, then ptwXY->points[0].x == x. */
765 status = ptwXY_lessEqualGreaterX_equal;
766 lessThanEqualXPoint->index = 0;
767 lessThanEqualXPoint->point = ptwXY->points[0]; }
768 else { /* ptwXY->length > 1 */
769 indexMin = 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; }
775 else {
776 indexMin = indexMid;
777 }
778 indexMid = ( indexMin + indexMax ) >> 1;
779 } // Loop checking, 11.06.2015, T. Koi
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]; }
788 else {
789 if( ptwXY->points[indexMin].x > x ) indexMax = 0;
790 if( ptwXY->points[indexMax].x < x ) indexMin = indexMax;
791 if( ( overflowPoint == overflowHeader ) || /* x < xMin of overflow points. */
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]; }
797 else {
798 lowerPoint = &(overflowPoint->point);
799 *lessThanEqualXPoint = *overflowPoint;
800 }
801 if( ( overflowPoint->next == overflowHeader ) || /* x > xMax of overflow points. */
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]; }
806 else {
807 upperPoint = &(overflowPoint->next->point);
808 *greaterThanXPoint = *(overflowPoint->next);
809 }
810 }
811 }
812 }
813 }
814
815 *closeIsEqual = 0;
816 if( eps > 0 ) {
817 double absX = std::fabs( x );
818
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 ) ) { /* x is closer to lower point. */
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; }
830 else { /* x is closer to upper point. */
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;
834 } }
835 else if( status == ptwXY_lessEqualGreaterX_equal ) {
836 *closeIsEqual = 1;
837 }
838 }
839 return( status );
840 }
841 /*
842 ************************************************************
843 */
844 nfu_status ptwXY_getValueAtX( ptwXYPoints *ptwXY, double x, double *y ) {
845
846 nfu_status status = nfu_XOutsideDomain;
847 ptwXYOverflowPoint lessThanEqualXPoint, greaterThanXPoint;
848 ptwXY_lessEqualGreaterX legx = ptwXY_getPointsAroundX( ptwXY, x, &lessThanEqualXPoint, &greaterThanXPoint );
849
850 *y = 0.;
851 if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
852 switch( legx ) {
853 case ptwXY_lessEqualGreaterX_empty :
854 case ptwXY_lessEqualGreaterX_lessThan :
855 case ptwXY_lessEqualGreaterX_greater :
856 break;
857 case ptwXY_lessEqualGreaterX_equal :
858 status = nfu_Okay;
859 *y = lessThanEqualXPoint.point.y;
860 break;
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 ); }
865 else {
866 status = ptwXY_interpolatePoint( ptwXY->interpolation, x, y, lessThanEqualXPoint.point.x, lessThanEqualXPoint.point.y,
867 greaterThanXPoint.point.x, greaterThanXPoint.point.y );
868 }
869 break;
870 }
871 return( status );
872 }
873 /*
874 ************************************************************
875 */
876 nfu_status ptwXY_setValueAtX( ptwXYPoints *ptwXY, double x, double y ) {
877
878 return( ptwXY_setValueAtX_overrideIfClose( ptwXY, x, y, 0., 0 ) );
879 }
880 /*
881 ************************************************************
882 */
883 nfu_status ptwXY_setValueAtX_overrideIfClose( ptwXYPoints *ptwXY, double x, double y, double eps, int override ) {
884
885 int closeIsEqual;
886 int64_t nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY ), i;
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 = NULL;
893
894 if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
895
896 legx = ptwXY_getPointsAroundX_closeIsEqual( ptwXY, x, &lessThanEqualXPoint, &greaterThanXPoint, eps, &closeIsEqual, &closePoint );
897 switch( legx ) {
898 case ptwXY_lessEqualGreaterX_lessThan :
899 case ptwXY_lessEqualGreaterX_greater :
900 case ptwXY_lessEqualGreaterX_between :
901 if( closeIsEqual ) {
902 if( !override ) return( status );
903 point = closePoint;
904 legx = ptwXY_lessEqualGreaterX_equal;
905 x = point->x; }
906 else {
907 if( ( legx == ptwXY_lessEqualGreaterX_greater ) && ( nonOverflowLength < ptwXY->allocatedSize ) ) {
908 point = &(ptwXY->points[nonOverflowLength]); }
909 else {
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; }
916 else { /* Between or greater and must go in overflow area. */
917 if( legx == ptwXY_lessEqualGreaterX_greater ) {
918 overflowPoint->prior = overflowHeader->prior;
919 overflowPoint->index = ptwXY->length; }
920 else {
921 overflowPoint->prior = lessThanEqualXPoint.prior;
922 if( lessThanEqualXPoint.next != NULL ) {
923 if( lessThanEqualXPoint.point.x < x )
924 overflowPoint->prior = lessThanEqualXPoint.prior->next;
925 i = 1; }
926 else {
927 for( p = overflowHeader->next, i = 1; p != overflowHeader; p = p->next, i++ )
928 if( p->point.x > x ) break;
929 }
930 overflowPoint->index = lessThanEqualXPoint.index + i;
931 }
932 }
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++;
939 }
940 ptwXY->overflowLength++;
941 }
942 }
943 break;
944 case ptwXY_lessEqualGreaterX_empty :
945 point = ptwXY->points; /* ptwXY_minimumSize must be > 0 so there is always space here. */
946 break;
947 case ptwXY_lessEqualGreaterX_equal :
948 if( closeIsEqual && !override ) return( status );
949 if( lessThanEqualXPoint.next == NULL ) {
950 point = &(ptwXY->points[lessThanEqualXPoint.index]); }
951 else {
952 point = &(lessThanEqualXPoint.prior->next->point);
953 }
954 break;
955 }
956 if( status == nfu_Okay ) {
957 point->x = x;
958 point->y = y;
959 if( legx != ptwXY_lessEqualGreaterX_equal ) ptwXY->length++;
960 }
961 return( status );
962 }
963 /*
964 ************************************************************
965 */
966 nfu_status ptwXY_mergeFromXsAndYs( ptwXYPoints *ptwXY, int length, double *xs, double *ys ) {
967
968 return( ptwXY_mergeFrom( ptwXY, 1, length, xs, ys ) );
969 }
970 /*
971 ************************************************************
972 */
973 nfu_status ptwXY_mergeFromXYs( ptwXYPoints *ptwXY, int length, double *xys ) {
974
975 int i;
976 double *xs, *p1, *p2;
977 nfu_status status;
978
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;
983 status = ptwXY_mergeFrom( ptwXY, 2, length, xs, xys );
984 nfu_free( xs );
985
986 return( status );
987 }
988 /*
989 ************************************************************
990 */
991 static nfu_status ptwXY_mergeFrom( ptwXYPoints *ptwXY, int /*incY*/, int length, double *xs, double *ys ) {
992
993 int i, j, n;
994 double *sortedXs, *p1, *p2;
995 nfu_status status;
996 ptwXYPoint *point1, *point2;
997
998 if( length < 0 ) return( nfu_badInput );
999 if( length == 0 ) return( nfu_Okay );
1000 if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status );
1001
1002 //if( ( sortedXs = (double *) nfu_malloc( length * sizeof( double * ) ) ) == NULL ) return( nfu_mallocError );
1003 //TK fixed for Coverity 63081
1004 if( ( sortedXs = (double *) nfu_malloc( length * sizeof( double ) ) ) == NULL ) return( nfu_mallocError );
1005
1006 for( i = 0, p1 = sortedXs, p2 = xs; i < length; i++, p1++, p2++ ) *p1 = *p2;
1007 //qsort( sortedXs, length, sizeof( double * ), ptwXY_mergeCompareFunction );
1008 //TK fixed for Coverity 63079
1009 qsort( sortedXs, length, sizeof( double ), ptwXY_mergeCompareFunction );
1010
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;
1014 }
1015 if( j == ptwXY->length ) break; /* Completed all ptwXY points. */
1016 }
1017 n += (int) ( ( length - i ) + ( ptwXY->length - j ) );
1018
1019 if( ( status = ptwXY_reallocatePoints( ptwXY, n, 0 ) ) == nfu_Okay ) {
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 ) {
1024 point1->x = *p1;
1025 point1->y = ys[(int)(p1 - xs)];
1026 if( *p1 >= point2->x ) {
1027 point2++;
1028 j++;
1029 }
1030 p1--;
1031 i--; }
1032 else {
1033 *point1 = *point2;
1034 point2++;
1035 j++;
1036 }
1037 }
1038 for( ; i < length; i++, p1--, point1-- ) {
1039 point1->x = *p1;
1040 point1->y = ys[(int)(p1 - xs)];
1041 }
1042 for( ; j < length; j++, point1--, point2-- ) *point1 = *point2;
1043 }
1044 nfu_free( sortedXs );
1045
1046 return( status );
1047 }
1048 /*
1049 ************************************************************
1050 */
1051 static int ptwXY_mergeCompareFunction( void const *x1p, void const *x2p ) {
1052
1053 double d1 = *((double *) x1p), d2 = *((double *) x2p);
1054
1055 if( d1 < d2 ) return( -1 );
1056 if( d1 == d2 ) return( 0 );
1057 return( 1 );
1058 }
1059 /*
1060 ************************************************************
1061 */
1062 nfu_status ptwXY_appendXY( ptwXYPoints *ptwXY, double x, double y ) {
1063
1064 int64_t nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
1065 ptwXY_dataFrom dataFrom;
1066
1067 if( ptwXY->length != 0 ) {
1068 double xMax = ptwXY_getXMaxAndFrom( ptwXY, &dataFrom );
1069 if( xMax >= x ) return( nfu_XNotAscending );
1070 }
1071
1072 if( nonOverflowLength < ptwXY->allocatedSize ) { /* Room at end of points. Also handles the case when length = 0. */
1073 ptwXY->points[nonOverflowLength].x = x;
1074 ptwXY->points[nonOverflowLength].y = y; }
1075 else {
1076 if( ptwXY->overflowLength == ptwXY->overflowAllocatedSize ) {
1077 ptwXYPoint newPoint = { x, y };
1078 return( ptwXY_coalescePoints( ptwXY, ptwXY->length + ptwXY->overflowAllocatedSize, &newPoint, 0 ) ); }
1079 else { /* Add to end of overflow. */
1080 ptwXYOverflowPoint *overflowPoint = &(ptwXY->overflowPoints[ptwXY->overflowLength]);
1081
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++;
1090 }
1091 }
1092 ptwXY->length++;
1093 return( nfu_Okay );
1094 }
1095 /*
1096 ************************************************************
1097 */
1098 nfu_status ptwXY_setXYPairAtIndex( ptwXYPoints *ptwXY, int64_t index, double x, double y ) {
1099
1100 int64_t i, ip1;
1101 ptwXYOverflowPoint *overflowPoint, *pm1, *pp1;
1102
1103 if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
1104
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;
1108 }
1109 ip1 = i;
1110 pm1 = pp1 = overflowPoint;
1111 if( overflowPoint->index == index ) { /* Note, if overflowPoint is header, then its index = -1. */
1112 pp1 = overflowPoint->next;
1113 ip1++;
1114 }
1115 if( ( pp1 != &(ptwXY->overflowHeader) ) && ( pp1->index == ( index + 1 ) ) ) { /* This if and else check that x < element[index+1]'s x values. */
1116 if( pp1->point.x <= x ) return( nfu_badIndexForX ); }
1117 else {
1118 if( ( ( index + 1 ) < ptwXY->length ) && ( ptwXY->points[index + 1 - ip1].x <= x ) ) return( nfu_badIndexForX );
1119 }
1120 if( overflowPoint != &(ptwXY->overflowHeader) ) pm1 = overflowPoint->prior;
1121 if( ( pm1 != &(ptwXY->overflowHeader) ) && ( pm1->index == ( index - 1 ) ) ) { /* This if and else check that x > element[index-1]'s x values. */
1122 if( pm1->point.x >= x ) return( nfu_badIndexForX ); }
1123 else {
1124 if( ( ( index - 1 ) >= 0 ) && ( ptwXY->points[index - 1 - i].x >= x ) ) return( nfu_badIndexForX );
1125 }
1126 if( ( overflowPoint != &(ptwXY->overflowHeader) ) && ( overflowPoint->index == index ) ) {
1127 overflowPoint->point.x = x;
1128 overflowPoint->point.y = y; }
1129 else {
1130 index -= i;
1131 ptwXY->points[index].x = x;
1132 ptwXY->points[index].y = y;
1133 }
1134 return( nfu_Okay );
1135 }
1136 /*
1137 ************************************************************
1138 */
1139 nfu_status ptwXY_getSlopeAtX( ptwXYPoints *ptwXY, double x, const char side, double *slope ) {
1140
1141 nfu_status status = nfu_Okay;
1142 ptwXYOverflowPoint lessThanEqualXPoint, greaterThanXPoint;
1143 ptwXY_lessEqualGreaterX legx = ptwXY_getPointsAroundX( ptwXY, x, &lessThanEqualXPoint, &greaterThanXPoint );
1144 ptwXYPoint *point;
1145
1146 *slope = 0.;
1147 if( ( side != '-' ) && ( side != '+' ) ) return( nfu_badInput );
1148
1149 switch( legx ) {
1150 case ptwXY_lessEqualGreaterX_empty :
1151 case ptwXY_lessEqualGreaterX_lessThan :
1152 case ptwXY_lessEqualGreaterX_greater :
1153 status = nfu_XOutsideDomain;
1154 break;
1155 case ptwXY_lessEqualGreaterX_between :
1156 *slope = ( greaterThanXPoint.point.y - lessThanEqualXPoint.point.y ) /
1157 ( greaterThanXPoint.point.x - lessThanEqualXPoint.point.x );
1158 break;
1159 case ptwXY_lessEqualGreaterX_equal :
1160 if( side == '-' ) {
1161 if( lessThanEqualXPoint.index == 0 ) {
1162 status = nfu_XOutsideDomain; }
1163 else {
1164 point = ptwXY_getPointAtIndex_Unsafely( ptwXY, lessThanEqualXPoint.index - 1 );
1165 *slope = ( lessThanEqualXPoint.point.y - point->y ) / ( lessThanEqualXPoint.point.x - point->x );
1166 } }
1167 else {
1168 if( lessThanEqualXPoint.index == ( ptwXY->length - 1 ) ) {
1169 status = nfu_XOutsideDomain; }
1170 else {
1171 point = ptwXY_getPointAtIndex_Unsafely( ptwXY, lessThanEqualXPoint.index + 1 );
1172 *slope = ( point->y - lessThanEqualXPoint.point.y ) / ( point->x - lessThanEqualXPoint.point.x );
1173 }
1174 }
1175 }
1176
1177 return( status );
1178 }
1179 /*
1180 ************************************************************
1181 */
1182 double ptwXY_getXMinAndFrom( ptwXYPoints *ptwXY, ptwXY_dataFrom *dataFrom ) {
1183
1184 int64_t nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
1185 double xMin = nfu_getNAN( );
1186
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;
1195 }
1196 } }
1197 else if( nonOverflowLength > 0 ) {
1198 *dataFrom = ptwXY_dataFrom_Points;
1199 xMin = ptwXY->points[0].x;
1200 }
1201 return( xMin );
1202 }
1203 /*
1204 ************************************************************
1205 */
1206 double ptwXY_getXMin( ptwXYPoints *ptwXY ) {
1207
1208 ptwXY_dataFrom dataFrom;
1209
1210 return( ptwXY_getXMinAndFrom( ptwXY, &dataFrom ) );
1211 }
1212 /*
1213 ************************************************************
1214 */
1215 double ptwXY_getXMaxAndFrom( ptwXYPoints *ptwXY, ptwXY_dataFrom *dataFrom ) {
1216
1217 int64_t nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
1218 double xMax = nfu_getNAN( );
1219
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;
1228 }
1229 } }
1230 else if( ptwXY->length > 0 ) {
1231 *dataFrom = ptwXY_dataFrom_Points;
1232 xMax = ptwXY->points[nonOverflowLength-1].x;
1233 }
1234 return( xMax );
1235 }
1236 /*
1237 ************************************************************
1238 */
1239 double ptwXY_getXMax( ptwXYPoints *ptwXY ) {
1240
1241 ptwXY_dataFrom dataFrom;
1242
1243 return( ptwXY_getXMaxAndFrom( ptwXY, &dataFrom ) );
1244 }
1245 /*
1246 ************************************************************
1247 */
1248 double ptwXY_getYMin( ptwXYPoints *ptwXY ) {
1249
1250 int64_t i, n = ptwXY_getNonOverflowLength( ptwXY );
1251 ptwXYPoint *p = ptwXY->points;
1252 ptwXYOverflowPoint *overflowPoint = ptwXY->overflowHeader.next;
1253 double yMin;
1254
1255 if( ptwXY->length == 0 ) return( 0. );
1256 if( n > 0 ) {
1257 yMin = p->y;
1258 for( i = 1, p++; i < n; i++, p++ ) yMin = ( ( yMin < p->y ) ? yMin : p->y ); }
1259 else {
1260 yMin = overflowPoint->point.y;
1261 }
1262 for( ; overflowPoint != &(ptwXY->overflowHeader); overflowPoint = overflowPoint->next )
1263 yMin = ( ( yMin < overflowPoint->point.y ) ? yMin : overflowPoint->point.y );
1264 return( yMin );
1265 }
1266 /*
1267 ************************************************************
1268 */
1269 double ptwXY_getYMax( ptwXYPoints *ptwXY ) {
1270
1271 int64_t i, n = ptwXY_getNonOverflowLength( ptwXY );
1272 ptwXYPoint *p = ptwXY->points;
1273 ptwXYOverflowPoint *overflowPoint = ptwXY->overflowHeader.next;
1274 double yMax;
1275
1276 if( ptwXY->length == 0 ) return( 0. );
1277 if( n > 0 ) {
1278 yMax = p->y;
1279 for( i = 1, p++; i < n; i++, p++ ) yMax = ( ( yMax > p->y ) ? yMax : p->y ); }
1280 else {
1281 yMax = overflowPoint->point.y;
1282 }
1283 for( ; overflowPoint != &(ptwXY->overflowHeader); overflowPoint = overflowPoint->next )
1284 yMax = ( ( yMax > overflowPoint->point.y ) ? yMax : overflowPoint->point.y );
1285 return( yMax );
1286 }
1287 /*
1288 ************************************************************
1289 */
1290 static void ptwXY_initialOverflowPoint( ptwXYOverflowPoint *overflowPoint, ptwXYOverflowPoint *prior, ptwXYOverflowPoint *next ) {
1291
1292 overflowPoint->prior = prior;
1293 overflowPoint->next = next;
1294 overflowPoint->index = -1;
1295 overflowPoint->point.x = 0.;
1296 overflowPoint->point.y = 0.;
1297 }
1298
1299 #if defined __cplusplus
1300 }
1301 #endif
1302