| Geant4 Cross Reference |
1 /*
2 # <<BEGIN-copyright>>
3 # <<END-copyright>>
4 */
5
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <cmath>
9 #include <cstdlib> //for std::abs(int)
10 #include <float.h>
11
12 #include "ptwX.h"
13
14 #if defined __cplusplus
15 namespace GIDI {
16 using namespace GIDI;
17 #endif
18
19 static int ptwX_sort_descending( void const *p1, void const *p2 );
20 static int ptwX_sort_ascending( void const *p1, void const *p2 );
21 /*
22 ************************************************************
23 */
24 ptwXPoints *ptwX_new( int64_t size, nfu_status *status ) {
25
26 ptwXPoints *ptwX = (ptwXPoints *) nfu_calloc( sizeof( ptwXPoints ), 1 );
27
28 *status = nfu_mallocError;
29 if( ptwX == NULL ) return( NULL );
30 ptwX_setup( ptwX, size );
31 if( ( *status = ptwX->status ) != nfu_Okay ) ptwX = (ptwXPoints *) nfu_free( ptwX );
32 return( ptwX );
33 }
34 /*
35 ************************************************************
36 */
37 nfu_status ptwX_setup( ptwXPoints *ptwX, int64_t size ) {
38
39 ptwX->status = nfu_Okay;
40 ptwX->length = 0;
41 ptwX->allocatedSize = 0;
42 ptwX->mallocFailedSize = 0;
43 ptwX->points = NULL;
44 ptwX_reallocatePoints( ptwX, size, 0 );
45 return( ptwX->status );
46 }
47 /*
48 ************************************************************
49 */
50 ptwXPoints *ptwX_create( int64_t size, int64_t length, double const *xs, nfu_status *status ) {
51
52 ptwXPoints *ptwX = ptwX_new( size, status );
53
54 if( ptwX != NULL ) {
55 if( ( *status = ptwX_setData( ptwX, length, xs ) ) != nfu_Okay ) ptwX = ptwX_free( ptwX );
56 }
57 return( ptwX );
58 }
59 /*
60 ************************************************************
61 */
62 ptwXPoints *ptwX_createLine( int64_t size, int64_t length, double slope, double offset, nfu_status *status ) {
63
64 int64_t i1;
65 double *p1;
66 ptwXPoints *ptwX;
67
68 if( size < length ) size = length;
69 if( ( ptwX = ptwX_new( size, status ) ) != NULL ) {
70 for( i1 = 0, p1 = ptwX->points; i1 < length; i1++, p1++ ) *p1 = slope * i1 + offset;
71 ptwX->length = length;
72 }
73 return( ptwX );
74 }
75 /*
76 ************************************************************
77 */
78 nfu_status ptwX_copy( ptwXPoints *dest, ptwXPoints *src ) {
79
80 if( dest->status == nfu_Okay ) return( dest->status );
81 if( src->status == nfu_Okay ) return( src->status );
82 ptwX_clear( dest );
83 return( ptwX_setData( dest, src->length, src->points ) );
84 }
85 /*
86 ************************************************************
87 */
88 ptwXPoints *ptwX_clone( ptwXPoints *ptwX, nfu_status *status ) {
89
90 return( ptwX_slice( ptwX, 0, ptwX->length, status ) );
91 }
92 /*
93 ************************************************************
94 */
95 ptwXPoints *ptwX_slice( ptwXPoints *ptwX, int64_t index1, int64_t index2, nfu_status *status ) {
96
97 int64_t i, j, length;
98 ptwXPoints *n;
99
100 *status = nfu_badSelf;
101 if( ptwX->status != nfu_Okay ) return( NULL );
102 *status = nfu_badIndex;
103 if( index1 < 0 ) return( NULL );
104 if( index2 < index1 ) return( NULL );
105 if( index2 > ptwX->length ) return( NULL );
106 length = ( index2 - index1 );
107 if( ( n = ptwX_new( length, status ) ) == NULL ) return( n );
108 *status = n->status;
109 for( j = 0, i = index1; i < index2; i++, j++ ) n->points[j] = ptwX->points[i];
110 n->length = length;
111 return( n );
112 }
113 /*
114 ************************************************************
115 */
116 nfu_status ptwX_reallocatePoints( ptwXPoints *ptwX, int64_t size, int forceSmallerResize ) {
117
118 if( size < ptwX_minimumSize ) size = ptwX_minimumSize; /* ptwX_minimumSize must be > 0 for other routines to work properly. */
119 if( size < ptwX->length ) size = ptwX->length;
120 if( size != ptwX->allocatedSize ) {
121 if( size > ptwX->allocatedSize ) { /* Increase size of allocated points. */
122 ptwX->points = (double *) nfu_realloc( (size_t) size * sizeof( double ), ptwX->points ); }
123 else if( ( ptwX->allocatedSize > 2 * size ) || forceSmallerResize ) { /* Decrease size, if at least 1/2 size reduction or if forced to. */
124 ptwX->points = (double *) nfu_realloc( (size_t) size * sizeof( double ), ptwX->points );
125 }
126 if( ptwX->points == NULL ) {
127 ptwX->mallocFailedSize = size;
128 size = 0;
129 ptwX->status = nfu_mallocError;
130 }
131 ptwX->allocatedSize = size;
132 }
133
134 return( ptwX->status );
135 }
136 /*
137 ************************************************************
138 */
139 nfu_status ptwX_clear( ptwXPoints *ptwX ) {
140
141 ptwX->length = 0;
142 return( ptwX->status );
143 }
144 /*
145 ************************************************************
146 */
147 nfu_status ptwX_release( ptwXPoints *ptwX ) {
148
149 ptwX->length = 0;
150 ptwX->allocatedSize = 0;
151 ptwX->points = (double *) nfu_free( ptwX->points );
152
153 return( nfu_Okay );
154 }
155 /*
156 ************************************************************
157 */
158 ptwXPoints *ptwX_free( ptwXPoints *ptwX ) {
159
160 if( ptwX != NULL ) ptwX_release( ptwX );
161 return( (ptwXPoints *) nfu_free( ptwX ) );
162 }
163 /*
164 ************************************************************
165 */
166 int64_t ptwX_length( ptwXPoints *ptwX ) {
167
168 return( ptwX->length );
169 }
170 /*
171 ************************************************************
172 */
173 nfu_status ptwX_setData( ptwXPoints *ptwX, int64_t length, double const *xs ) {
174
175 int64_t i;
176
177 if( ptwX->status != nfu_Okay ) return( ptwX->status );
178
179 if( length > ptwX->allocatedSize ) {
180 ptwX_reallocatePoints( ptwX, length, 0 );
181 if( ptwX->status != nfu_Okay ) return( ptwX->status );
182 }
183 for( i = 0; i < length; i++ ) ptwX->points[i] = xs[i];
184 ptwX->length = length;
185
186 return( ptwX->status );
187 }
188 /*
189 ************************************************************
190 */
191 nfu_status ptwX_deletePoints( ptwXPoints *ptwX, int64_t i1, int64_t i2 ) {
192
193 int64_t n = ptwX->length - ( i2 - i1 );
194
195 if( ptwX->status != nfu_Okay ) return( ptwX->status );
196 if( ( i1 < 0 ) || ( i1 > i2 ) || ( i2 > ptwX->length ) ) return( nfu_badIndex );
197 if( i1 != i2 ) {
198 for( ; i2 < ptwX->length; i1++, i2++ ) ptwX->points[i1] = ptwX->points[i2];
199 ptwX->length = n;
200 }
201 return( ptwX->status );
202 }
203 /*
204 ************************************************************
205 */
206 double *ptwX_getPointAtIndex( ptwXPoints *ptwX, int64_t index ) {
207
208 if( ptwX->status != nfu_Okay ) return( NULL );
209 if( ( index < 0 ) || ( index >= ptwX->length ) ) return( NULL );
210 return( &(ptwX->points[index]) );
211 }
212 /*
213 ************************************************************
214 */
215 double ptwX_getPointAtIndex_Unsafely( ptwXPoints *ptwX, int64_t index ) {
216
217 return( ptwX->points[index] );
218 }
219 /*
220 ************************************************************
221 */
222 nfu_status ptwX_setPointAtIndex( ptwXPoints *ptwX, int64_t index, double x ) {
223
224 nfu_status status;
225
226 if( ptwX->status != nfu_Okay ) return( ptwX->status );
227 if( ( index < 0 ) || ( index > ptwX->length ) ) return( nfu_badIndex );
228 if( index == ptwX->allocatedSize ) {
229 if( ( status = ptwX_reallocatePoints( ptwX, ptwX->allocatedSize + 10, 0 ) ) != nfu_Okay ) return( status );
230 }
231 ptwX->points[index] = x;
232 if( index == ptwX->length ) ptwX->length++;
233 return( nfu_Okay );
234 }
235 /*
236 ************************************************************
237 */
238 nfu_status ptwX_insertPointsAtIndex( ptwXPoints *ptwX, int64_t index, int64_t n1, double const *xs ) {
239
240 nfu_status status;
241 int64_t i1, i2, n1p, size = n1 + ptwX->length;
242
243 if( ptwX->status != nfu_Okay ) return( ptwX->status );
244 if( n1 < 1 ) return( nfu_Okay );
245 if( ( index < 0 ) || ( index > ptwX->length ) ) return( nfu_badIndex );
246 if( size > ptwX->allocatedSize ) {
247 if( ( status = ptwX_reallocatePoints( ptwX, size, 0 ) ) != nfu_Okay ) return( status );
248 }
249 for( i1 = ptwX->length - 1, i2 = size - 1, n1p = ptwX->length - index + 1; n1p > 0; i1--, i2--, n1p-- ) ptwX->points[i2] = ptwX->points[i1];
250 for( i1 = 0, i2 = index; i1 < n1; i1++, i2++ ) ptwX->points[i2] = xs[i1];
251 ptwX->length += n1;
252 return( nfu_Okay );
253 }
254 /*
255 ************************************************************
256 */
257 int ptwX_ascendingOrder( ptwXPoints *ptwX ) {
258 /*
259 * Returns -1 list is descending, 1 if ascending and 0 otherwise (i.e., mixed).
260 */
261 int order = 1;
262 int64_t i;
263 double x1, x2;
264
265 if( ptwX->length < 2 ) return( 0 );
266
267 if( ( x1 = ptwX->points[0] ) < ( x2 = ptwX->points[1] ) ) { /* Check for ascending order. */
268 for( i = 2; i < ptwX->length; i++ ) {
269 x1 = x2;
270 x2 = ptwX->points[i];
271 if( x2 <= x1 ) return( 0 );
272 } }
273 else {
274 if( x1 == x2 ) return( 0 );
275 order = -1; /* Check for descending order. */
276 for( i = 2; i < ptwX->length; i++ ) {
277 x1 = x2;
278 x2 = ptwX->points[i];
279 if( x1 <= x2 ) return( 0 );
280 }
281 }
282 return( order );
283 }
284 /*
285 ************************************************************
286 */
287 ptwXPoints *ptwX_fromString( char const *str, char **endCharacter, nfu_status *status ) {
288
289 int64_t numberConverted;
290 double *doublePtr;
291 ptwXPoints *ptwX = NULL;
292
293 if( ( *status = nfu_stringToListOfDoubles( str, &numberConverted, &doublePtr, endCharacter ) ) != nfu_Okay ) return( NULL );
294 ptwX = ptwX_create( numberConverted, numberConverted, doublePtr, status );
295 nfu_free( doublePtr );
296 return( ptwX );
297 }
298 /*
299 ************************************************************
300 */
301 nfu_status ptwX_countOccurrences( ptwXPoints *ptwX, double value, int *count ) {
302
303 int64_t i1;
304
305 *count = 0;
306 for( i1 = 0; i1 < ptwX->length; i1++ ) {
307 if( ptwX->points[i1] == value ) (*count)++;
308 }
309 return( nfu_Okay );
310 }
311 /*
312 ************************************************************
313 */
314 nfu_status ptwX_reverse( ptwXPoints *ptwX ) {
315
316 int64_t i1, i2 = ptwX->length - 1, n1 = ptwX->length / 2;
317 double tmp;
318
319 for( i1 = 0; i1 < n1; i1++, i2-- ) {
320 tmp = ptwX->points[i1];
321 ptwX->points[i1] = ptwX->points[i2];
322 ptwX->points[i2] = tmp;
323 }
324 return( nfu_Okay );
325 }
326 /*
327 ************************************************************
328 */
329 nfu_status ptwX_sort( ptwXPoints *ptwX, enum ptwX_sort_order order ) {
330
331 int (*cmp)( void const *, void const * ) = ptwX_sort_descending;
332
333 if( order == ptwX_sort_order_ascending ) cmp = ptwX_sort_ascending;
334 qsort( ptwX->points, (size_t) ptwX->length, sizeof( ptwX->points[0] ), cmp );
335 return( nfu_Okay );
336 }
337 /*
338 ************************************************************
339 */
340 static int ptwX_sort_descending( void const *p1, void const *p2 ) { return( -ptwX_sort_ascending( p1, p2 ) ); }
341 static int ptwX_sort_ascending( void const *p1, void const *p2 ) {
342
343 double *d1 = (double *) p1, *d2 = (double *) p2;
344
345 if( *d1 < *d2 ) return( -1 );
346 if( *d1 == *d2 ) return( 0 );
347 return( 1 );
348 }
349 /*
350 ************************************************************
351 */
352 nfu_status ptwX_closesDifference( ptwXPoints *ptwX, double value, int64_t *index, double *difference ) {
353
354 return( ptwX_closesDifferenceInRange( ptwX, 0, ptwX->length, value, index, difference ) );
355 }
356 /*
357 ************************************************************
358 */
359 nfu_status ptwX_closesDifferenceInRange( ptwXPoints *ptwX, int64_t i1, int64_t i2, double value, int64_t *index, double *difference ) {
360 /*
361 * Finds the closes datum to value. If *difference is zero, datum is same as value.
362 */
363 double d1;
364
365 *index = -1;
366 *difference = -1;
367 if( ptwX->status != nfu_Okay ) return( ptwX->status );
368 if( i1 < 0 ) i1 = 0;
369 if( i2 > ptwX->length ) i2 = ptwX->length;
370 if( i1 >= i2 ) return( nfu_Okay );
371 *index = i1;
372 *difference = value - ptwX->points[i1];
373 for( i1++; i1 < i2; i1++ ) {
374 d1 = value - ptwX->points[i1];
375 if( std::fabs( *difference ) > std::fabs( d1 ) ) {
376 *index = i1;
377 *difference = d1;
378 }
379 }
380 return( nfu_Okay );
381 }
382 /*
383 ************************************************************
384 */
385 ptwXPoints *ptwX_unique( ptwXPoints *ptwX, int order, nfu_status *status ) {
386 /*
387 * If order < 0 order is descending, if order > 0 order is ascending, otherwise, order is the same as ptwX.
388 */
389 int64_t i1, i2, n1 = 0;
390 double x1, *p2;
391 ptwXPoints *ptwX2 = NULL;
392
393 if( order == 0 ) {
394 if( ( ptwX2 = ptwX_new( ptwX->length, status ) ) == NULL ) return( NULL );
395 for( i1 = 0; i1 < ptwX->length; i1++ ) {
396 x1 = ptwX->points[i1];
397 for( i2 = 0, p2 = ptwX2->points; i2 < ptwX2->length; i2++, p2++ ) {
398 if( *p2 == x1 ) break;
399 }
400 if( i2 == ptwX2->length ) {
401 ptwX2->points[ptwX2->length] = x1;
402 ptwX2->length++;
403 }
404 } }
405 else {
406 if( ( ptwX2 = ptwX_clone( ptwX, status ) ) == NULL ) return( NULL );
407 if( ( *status = ptwX_sort( ptwX2, ptwX_sort_order_ascending ) ) != nfu_Okay ) goto err;
408
409 if( ptwX2->length > 1 ) {
410 x1 = ptwX2->points[n1];
411 n1++;
412 for( i1 = 1; i1 < ptwX2->length; i1++ ) {
413 if( x1 != ptwX2->points[i1] ) {
414 x1 = ptwX2->points[i1];
415 ptwX2->points[n1] = x1;
416 n1++;
417 }
418 }
419 ptwX2->length = n1;
420 if( order < 0 ) {
421 if( ( *status = ptwX_sort( ptwX2, ptwX_sort_order_descending ) ) != nfu_Okay ) goto err;
422 }
423 }
424 }
425 return( ptwX2 );
426
427 err:
428 if( ptwX2 != NULL ) ptwX_free( ptwX2 );
429 return( NULL );
430 }
431 /*
432 ************************************************************
433 */
434 nfu_status ptwX_abs( ptwXPoints *ptwX ) {
435
436 int64_t i1;
437 double *p1;
438
439 if( ptwX->status != nfu_Okay ) return( ptwX->status );
440 for( i1 = 0, p1 = ptwX->points; i1 < ptwX->length; i1++, p1++ ) *p1 = std::fabs( *p1 );
441 return( nfu_Okay );
442 }
443 /*
444 ************************************************************
445 */
446 nfu_status ptwX_neg( ptwXPoints *ptwX ) {
447
448 return( ptwX_slopeOffset( ptwX, -1, 0 ) );
449 }
450 /*
451 ************************************************************
452 */
453 nfu_status ptwX_add_double( ptwXPoints *ptwX, double value ) {
454
455 return( ptwX_slopeOffset( ptwX, 1, value ) );
456 }
457 /*
458 ************************************************************
459 */
460 nfu_status ptwX_mul_double( ptwXPoints *ptwX, double value ) {
461
462 return( ptwX_slopeOffset( ptwX, value, 0 ) );
463 }
464 /*
465 ************************************************************
466 */
467 nfu_status ptwX_slopeOffset( ptwXPoints *ptwX, double slope, double offset ) {
468
469 int64_t i1;
470 double *p1;
471
472 if( ptwX->status != nfu_Okay ) return( ptwX->status );
473 for( i1 = 0, p1 = ptwX->points; i1 < ptwX->length; i1++, p1++ ) *p1 = slope * *p1 + offset;
474 return( nfu_Okay );
475 }
476 /*
477 ************************************************************
478 */
479 nfu_status ptwX_add_ptwX( ptwXPoints *ptwX1, ptwXPoints *ptwX2 ) {
480
481 int64_t i1;
482 double *p1 = ptwX1->points, *p2 = ptwX2->points;
483
484 if( ptwX1->status != nfu_Okay ) return( ptwX1->status );
485 if( ptwX2->status != nfu_Okay ) return( ptwX2->status );
486 if( ptwX1->length != ptwX2->length ) return( nfu_domainsNotMutual );
487
488 for( i1 = 0; i1 < ptwX1->length; i1++, p1++, p2++ ) *p1 += *p2;
489 return( nfu_Okay );
490 }
491 /*
492 ************************************************************
493 */
494 nfu_status ptwX_sub_ptwX( ptwXPoints *ptwX1, ptwXPoints *ptwX2 ) {
495
496 int64_t i1;
497 double *p1 = ptwX1->points, *p2 = ptwX2->points;
498
499 if( ptwX1->status != nfu_Okay ) return( ptwX1->status );
500 if( ptwX2->status != nfu_Okay ) return( ptwX2->status );
501 if( ptwX1->length != ptwX2->length ) return( nfu_domainsNotMutual );
502
503 for( i1 = 0; i1 < ptwX1->length; i1++, p1++, p2++ ) *p1 -= *p2;
504 return( nfu_Okay );
505 }
506 /*
507 ************************************************************
508 */
509 nfu_status ptwX_xMinMax( ptwXPoints *ptwX, double *xMin, double *xMax ) {
510
511 int64_t i1, n1 = ptwX->length;
512 *xMin = *xMax = 0;
513 double *p1 = ptwX->points;
514
515 if( ptwX->status != nfu_Okay ) return( ptwX->status );
516 if( n1 > 0 ) {
517 *xMin = *xMax = *(p1++);
518 for( i1 = 1; i1 < n1; ++i1, ++p1 ) {
519 if( *p1 < *xMin ) *xMin = *p1;
520 if( *p1 > *xMax ) *xMax = *p1;
521 }
522 }
523 return( nfu_Okay );
524 }
525 /*
526 ************************************************************
527 */
528 nfu_status ptwX_compare( ptwXPoints *ptwX1, ptwXPoints *ptwX2, int *comparison ) {
529
530 int64_t i1, n1 = ptwX1->length, n2 = ptwX2->length, nn = n1;
531 double *p1 = ptwX1->points, *p2 = ptwX2->points;
532
533 *comparison = 0;
534 if( ptwX1->status != nfu_Okay ) return( ptwX1->status );
535 if( ptwX2->status != nfu_Okay ) return( ptwX2->status );
536 if( nn > n2 ) nn = n2;
537 for( i1 = 0; i1 < nn; i1++, p1++, p2++ ) {
538 if( *p1 == *p2 ) continue;
539 *comparison = 1;
540 if( *p1 < *p2 ) *comparison = -1;
541 return( nfu_Okay );
542 }
543 if( n1 < n2 ) {
544 *comparison = -1; }
545 else if( n1 > n2 ) {
546 *comparison = 1;
547 }
548 return( nfu_Okay );
549 }
550 /*
551 ************************************************************
552 */
553 int ptwX_close( ptwXPoints *ptwX1, ptwXPoints *ptwX2, int epsilonFactor, double epsilon, nfu_status *status ) {
554
555 int64_t i1, n1 = ptwX1->length;
556 double larger;
557 double *p1 = ptwX1->points, *p2 = ptwX2->points;
558
559 epsilon = std::fabs( epsilon ) + std::abs( epsilonFactor ) * DBL_EPSILON;
560
561 *status = ptwX1->status;
562 if( ptwX1->status != nfu_Okay ) return( -1 );
563 *status = ptwX2->status;
564 if( ptwX2->status != nfu_Okay ) return( -1 );
565 *status = nfu_domainsNotMutual;
566 if( n1 != ptwX2->length ) return( -1 );
567
568 *status = nfu_Okay;
569 for( i1 = 0; i1 < n1; i1++, p1++, p2++ ) {
570 larger = std::fabs( *p1 );
571 if( std::fabs( *p2 ) > larger ) larger = std::fabs( *p2 );
572 if( std::fabs( *p2 - *p1 ) > epsilon * larger ) return( (int) ( i1 + 1 ) );
573 }
574 return( 0 );
575 }
576
577 #if defined __cplusplus
578 }
579 #endif
580