Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ptwXY_core.cc
Go to the documentation of this file.
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;
55  ptwXY->interpolationOtherInfo.argList = NULL;
56  switch( interpolation ) {
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 );
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 
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 */
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 */
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 ) {
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 ) {
170  return( dest->status = nfu_mallocError );
171  } }
172  else {
174  }
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 */
209 
210  return( ptwXY_slice( ptwXY, 0, ptwXY->length, ptwXY->overflowAllocatedSize, status ) );
211 }
212 /*
213 ************************************************************
214 */
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 ) {
225  n1->interpolation = interpolationTo;
226  switch( interpolationTo ) {
237  case ptwXY_interpolationOther : /* Does not happen, but needed to stop compilers from complaining. */
238  break;
239  }
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 */
338 
339  return( ptwXY->interpolation );
340 }
341 /*
342 ************************************************************
343 */
345 
347 }
348 /*
349 ************************************************************
350 */
352 
353  return( ptwXY->status );
354 }
355 /*
356 ************************************************************
357 */
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 */
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 */
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 */
530 
531  return( ptwXY_coalescePoints( ptwXY, ptwXY->length, NULL, 0 ) );
532 }
533 /*
534 ************************************************************
535 */
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 */
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 )
557  }
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;
568 
569  return( nfu_Okay );
570 }
571 /*
572 ************************************************************
573 */
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 */
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 */
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);
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 ) {
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 ) {
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 ) ) {
763  *lessThanEqualXPoint = *overflowPoint; }
764  else if( ptwXY->length == 1 ) { /* If here and length = 1, then ptwXY->points[0].x == x. */
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 ) {
782  lessThanEqualXPoint->index = indexMin;
783  lessThanEqualXPoint->point = ptwXY->points[indexMin]; }
784  else if( ptwXY->points[indexMax].x == x ) {
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 
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 ) {
856  break;
858  status = nfu_Okay;
859  *y = lessThanEqualXPoint.point.y;
860  break;
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;
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 ) {
901  if( closeIsEqual ) {
902  if( !override ) return( status );
903  point = closePoint;
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;
945  point = ptwXY->points; /* ptwXY_minimumSize must be > 0 so there is always space here. */
946  break;
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 ) {
1153  status = nfu_XOutsideDomain;
1154  break;
1156  *slope = ( greaterThanXPoint.point.y - lessThanEqualXPoint.point.y ) /
1157  ( greaterThanXPoint.point.x - lessThanEqualXPoint.point.x );
1158  break;
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 */
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
double minFractional_dx
Definition: ptwXY.h:92
#define ptwXY_minimumOverflowSize
Definition: ptwXY.h:21
static char const logLinInterpolationString[]
Definition: ptwXY_core.cc:19
ptwXYPoint * ptwXY_getPointAtIndex(ptwXYPoints *ptwXY, int64_t index)
Definition: ptwXY_core.cc:675
int64_t mallocFailedSize
Definition: ptwXY.h:97
nfu_status ptwXY_getValueAtX(ptwXYPoints *ptwXY, double x, double *y)
Definition: ptwXY_core.cc:844
double ptwXY_getAccuracy(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:372
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
nfu_status ptwXY_mergeFromXYs(ptwXYPoints *ptwXY, int length, double *xys)
Definition: ptwXY_core.cc:973
nfu_status ptwXY_coalescePoints(ptwXYPoints *ptwXY, int64_t size, ptwXYPoint *newPoint, int forceSmallerResize)
Definition: ptwXY_core.cc:469
ptwXYPoint point
Definition: ptwXY.h:80
nfu_status ptwXY_interpolatePoint(ptwXY_interpolation interpolation, double x, double *y, double x1, double y1, double x2, double y2)
int userFlag
Definition: ptwXY.h:89
ptwXYPoint * points
Definition: ptwXY.h:99
double ptwXY_setAccuracy(ptwXYPoints *ptwXY, double accuracy)
Definition: ptwXY_core.cc:379
double ptwXY_getXMaxAndFrom(ptwXYPoints *ptwXY, ptwXY_dataFrom *dataFrom)
Definition: ptwXY_core.cc:1215
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)
Definition: ptwXY_core.cc:126
static char const logLogInterpolationString[]
Definition: ptwXY_core.cc:20
const char * p
Definition: xmltok.h:285
static const G4double d2
ptwXYPoints * ptwXY_clone(ptwXYPoints *ptwXY, nfu_status *status)
Definition: ptwXY_core.cc:208
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)
Definition: ptwXY_core.cc:108
nfu_status ptwXY_deletePoints(ptwXYPoints *ptwXY, int64_t i1, int64_t i2)
Definition: ptwXY_core.cc:660
nfu_status ptwXY_reallocateOverflowPoints(ptwXYPoints *ptwXY, int64_t size)
Definition: ptwXY_core.cc:439
ptwXYPoints * ptwXY_cloneToInterpolation(ptwXYPoints *ptwXY, ptwXY_interpolation interpolationTo, nfu_status *status)
Definition: ptwXY_core.cc:215
nfu_status ptwXY_mergeFromXsAndYs(ptwXYPoints *ptwXY, int length, double *xs, double *ys)
Definition: ptwXY_core.cc:966
ptwXY_getValue_callback getValueFunc
Definition: ptwXY.h:71
static const G4double eps
enum ptwXY_lessEqualGreaterX_e ptwXY_lessEqualGreaterX
nfu_status ptwXY_appendXY(ptwXYPoints *ptwXY, double x, double y)
Definition: ptwXY_core.cc:1062
int ptwXY_getUserFlag(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:358
tuple x
Definition: test.py:50
ptwXYPoint * ptwXY_getPointAtIndex_Unsafely(ptwXYPoints *ptwXY, int64_t index)
Definition: ptwXY_core.cc:684
nfu_status ptwXY_copy(ptwXYPoints *dest, ptwXYPoints *src)
Definition: ptwXY_core.cc:148
int64_t ptwXY_length(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:583
struct ptwXYOverflowPoint_s * prior
Definition: ptwXY.h:77
double ptwXY_getYMin(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:1248
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
static void ptwXY_initialOverflowPoint(ptwXYOverflowPoint *overflowPoint, ptwXYOverflowPoint *prior, ptwXYOverflowPoint *next)
Definition: ptwXY_core.cc:1290
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
Definition: ptwXY_core.cc:876
double biSectionMax
Definition: ptwXY.h:90
int64_t length
Definition: ptwXY.h:93
void ptwXY_setUserFlag(ptwXYPoints *ptwXY, int userFlag)
Definition: ptwXY_core.cc:365
nfu_status ptwXY_setXYDataFromXsAndYs(ptwXYPoints *ptwXY, int64_t length, double const *x, double const *y)
Definition: ptwXY_core.cc:631
double ptwXY_getYMax(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:1269
void * nfu_calloc(size_t size, size_t n)
nfu_status ptwXY_release(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:549
double ptwXY_setBiSectionMax(ptwXYPoints *ptwXY, double biSectionMax)
Definition: ptwXY_core.cc:397
int64_t ptwXY_getNonOverflowLength(ptwXYPoints const *ptwXY)
Definition: ptwXY_core.cc:590
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)
Definition: ptwXY_core.cc:46
ptwXY_lessEqualGreaterX ptwXY_getPointsAroundX_closeIsEqual(ptwXYPoints *ptwXY, double x, ptwXYOverflowPoint *lessThanEqualXPoint, ptwXYOverflowPoint *greaterThanXPoint, double eps, int *closeIsEqual, ptwXYPoint **closePoint)
Definition: ptwXY_core.cc:720
nfu_status ptwXY_getStatus(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:351
enum ptwXY_dataFrom_e ptwXY_dataFrom
static char const linLinInterpolationString[]
Definition: ptwXY_core.cc:17
nfu_status ptwXY_setValueAtX_overrideIfClose(ptwXYPoints *ptwXY, double x, double y, double eps, int override)
Definition: ptwXY_core.cc:883
double nfu_getNAN(void)
Definition: nf_utilities.cc:54
nfu_status ptwXY_setXYPairAtIndex(ptwXYPoints *ptwXY, int64_t index, double x, double y)
Definition: ptwXY_core.cc:1098
enum nfu_status_e nfu_status
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
double ptwXY_getXMin(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:1206
const G4int n
void * nfu_free(void *p)
double ptwXY_getBiSectionMax(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:390
static const G4double d1
ptwXYPoints * ptwXY_xMinSlice(ptwXYPoints *ptwXY, double xMin, int64_t secondarySize, int fill, nfu_status *status)
Definition: ptwXY_core.cc:315
ptwXY_interpolationOtherInfo interpolationOtherInfo
Definition: ptwXY.h:88
enum ptwXY_interpolation_e ptwXY_interpolation
int64_t allocatedSize
Definition: ptwXY.h:94
nfu_status ptwXY_clear(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:536
ptwXY_sigma typeY
Definition: ptwXY.h:86
double y
Definition: ptwXY.h:62
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)
Definition: ptwXY_core.cc:29
struct ptwXYOverflowPoint_s * next
Definition: ptwXY.h:78
ptwXYPoints * ptwXY_xSlice(ptwXYPoints *ptwXY, double xMin, double xMax, int64_t secondarySize, int fill, nfu_status *status)
Definition: ptwXY_core.cc:274
static nfu_status ptwXY_mergeFrom(ptwXYPoints *ptwXY, int incY, int length, double *xs, double *ys)
Definition: ptwXY_core.cc:991
double x
Definition: ptwXY.h:62
int64_t overflowAllocatedSize
Definition: ptwXY.h:96
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
ptwXYPoints * ptwXY_slice(ptwXYPoints *ptwXY, int64_t index1, int64_t index2, int64_t secondarySize, nfu_status *status)
Definition: ptwXY_core.cc:248
int64_t overflowLength
Definition: ptwXY.h:95
nfu_status status
Definition: ptwXY.h:85
#define ptwXY_maxBiSectionMax
Definition: ptwXY.h:22
char const * ptwXY_getInterpolationString(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:344
nfu_status ptwXY_setXYData(ptwXYPoints *ptwXY, int64_t length, double const *xy)
Definition: ptwXY_core.cc:597
nfu_status ptwXY_getSlopeAtX(ptwXYPoints *ptwXY, double x, const char side, double *slope)
Definition: ptwXY_core.cc:1139
double ptwXY_getXMinAndFrom(ptwXYPoints *ptwXY, ptwXY_dataFrom *dataFrom)
Definition: ptwXY_core.cc:1182
void * nfu_realloc(size_t size, void *old)
char const * interpolationString
Definition: ptwXY.h:70
ptwXYOverflowPoint overflowHeader
Definition: ptwXY.h:98
ptwXY_interpolation ptwXY_getInterpolation(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:337
int64_t index
Definition: ptwXY.h:79
struct ptwXYOverflowPoint_s ptwXYOverflowPoint
double accuracy
Definition: ptwXY.h:91
#define ptwXY_minimumSize
Definition: ptwXY.h:20
static char const flatInterpolationString[]
Definition: ptwXY_core.cc:21
void * nfu_malloc(size_t size)
double ptwXY_getXMax(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:1239
nfu_status ptwXY_reallocatePoints(ptwXYPoints *ptwXY, int64_t size, int forceSmallerResize)
Definition: ptwXY_core.cc:410
ptwXYPoints * ptwXY_xMaxSlice(ptwXYPoints *ptwXY, double xMax, int64_t secondarySize, int fill, nfu_status *status)
Definition: ptwXY_core.cc:326
struct ptwXYPoint_s ptwXYPoint
ptwXY_lessEqualGreaterX ptwXY_getPointsAroundX(ptwXYPoints *ptwXY, double x, ptwXYOverflowPoint *lessThanEqualXPoint, ptwXYOverflowPoint *greaterThanXPoint)
Definition: ptwXY_core.cc:710
#define ptwXY_minAccuracy
Definition: ptwXY.h:23
ptwXY_sigma typeX
Definition: ptwXY.h:86
ptwXYOverflowPoint * overflowPoints
Definition: ptwXY.h:100
nfu_status ptwXY_getXYPairAtIndex(ptwXYPoints *ptwXY, int64_t index, double *x, double *y)
Definition: ptwXY_core.cc:698
static char const linLogInterpolationString[]
Definition: ptwXY_core.cc:18
static int ptwXY_mergeCompareFunction(void const *x1p, void const *x2p)
Definition: ptwXY_core.cc:1051