| Geant4 Cross Reference |
1 /*
2 # <<BEGIN-copyright>>
3 # <<END-copyright>>
4 */
5
6 #include <cmath>
7 #include <float.h>
8
9 #include "ptwXY.h"
10
11 #if defined __cplusplus
12 #include <cmath>
13 #include "G4Exp.hh"
14 #include "G4Log.hh"
15 namespace GIDI {
16 using namespace GIDI;
17 #endif
18
19 static nfu_status ptwXY_clip2( ptwXYPoints *ptwXY1, double y, double x1, double y1, double x2, double y2 );
20 static double ptwXY_thicken_linear_dx( int sectionSubdivideMax, double dxMax, double x1, double x2 );
21 static nfu_status ptwXY_thin2( ptwXYPoints *thinned, char *thin, double accuracy, int64_t i1, int64_t i2 );
22 /*
23 ************************************************************
24 */
25 nfu_status ptwXY_clip( ptwXYPoints *ptwXY1, double yMin, double yMax ) {
26 /*
27 This function acts oddly for xy = [ [ 1, 0 ], [ 3, -2 ], [ 4, 1 ] ] and yMin = 0.2, why???????
28 This function probably only works for linear, linear interpolation (mainly because of ptwXY_clip2).
29 */
30 int64_t i, j, n;
31 double x2, y2;
32 nfu_status status;
33 ptwXYPoints *clipped;
34 ptwXYPoint *points;
35
36 if( ( status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( status );
37 if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
38 n = ptwXY1->length;
39 if( n > 0 ) {
40 i = 0;
41 if( ptwXY_getYMax( ptwXY1 ) < yMin ) i = 1;
42 if( ptwXY_getYMin( ptwXY1 ) > yMax ) i = 1;
43 if( i == 1 ) return( ptwXY_clear( ptwXY1 ) );
44 }
45 if( n == 1 ) {
46 y2 = ptwXY1->points[0].y;
47 if( y2 < yMin ) {
48 ptwXY1->points[0].y = yMin; }
49 else if( y2 > yMax ) {
50 ptwXY1->points[0].y = yMax;
51 } }
52 else if( n > 1 ) {
53 if( ( clipped = ptwXY_new( ptwXY1->interpolation, &(ptwXY1->interpolationOtherInfo),
54 ptwXY1->biSectionMax, ptwXY1->accuracy, n, 10, &status, ptwXY1->userFlag ) ) == NULL )
55 return( ptwXY1->status = status );
56 for( i = 0; i < n; i++ ) {
57 x2 = ptwXY1->points[i].x;
58 y2 = ptwXY1->points[i].y;
59 if( y2 < yMin ) {
60 if( i > 0 ) {
61 points = ptwXY_getPointAtIndex_Unsafely( clipped, clipped->length - 1 );
62 if( points->y > yMin ) {
63 if( ( status = ptwXY_clip2( clipped, yMin, points->x, points->y, x2, y2 ) ) != nfu_Okay ) goto Err;
64 }
65 }
66 if( ( status = ptwXY_setValueAtX( clipped, x2, yMin ) ) != nfu_Okay ) goto Err;
67 j = i;
68 for( i++; i < n; i++ ) if( !( ptwXY1->points[i].y < yMin ) ) break;
69 if( i < n ) {
70 x2 = ptwXY1->points[i].x;
71 y2 = ptwXY1->points[i].y;
72 if( ( status = ptwXY_clip2( clipped, yMin, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) goto Err;
73 if( y2 > yMax ) {
74 if( ( status = ptwXY_clip2( clipped, yMax, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) goto Err;
75 } }
76 else if( j != n - 1 ) {
77 if( ( status = ptwXY_setValueAtX( clipped, ptwXY1->points[n - 1].x, yMin ) ) != nfu_Okay ) goto Err;
78 }
79 i--; }
80 else if( y2 > yMax ) {
81 if( i > 0 ) {
82 points = ptwXY_getPointAtIndex_Unsafely( clipped, clipped->length - 1 );
83 if( points->y < yMax ) {
84 if( ( status = ptwXY_clip2( clipped, yMax, points->x, points->y, x2, y2 ) ) != nfu_Okay ) goto Err;
85 }
86 }
87 if( ( status = ptwXY_setValueAtX( clipped, x2, yMax ) ) != nfu_Okay ) goto Err;
88 j = i;
89 for( i++; i < n; i++ ) if( !( ptwXY1->points[i].y > yMax ) ) break;
90 if( i < n ) {
91 x2 = ptwXY1->points[i].x;
92 y2 = ptwXY1->points[i].y;
93 if( ( status = ptwXY_clip2( clipped, yMax, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) goto Err;
94 if( y2 < yMin ) {
95 if( ( status = ptwXY_clip2( clipped, yMin, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) goto Err;
96 } }
97 else if( j != n - 1 ) {
98 if( ( status = ptwXY_setValueAtX( clipped, ptwXY1->points[n - 1].x, yMax ) ) != nfu_Okay ) goto Err;
99 }
100 i--; }
101 else {
102 if( ( status = ptwXY_setValueAtX( clipped, x2, y2 ) ) != nfu_Okay ) goto Err;
103 }
104 }
105 if( ( status = ptwXY_simpleCoalescePoints( clipped ) ) != nfu_Okay ) goto Err;
106 ptwXY1->length = clipped->length; /* The squeamish may want to skip the next few lines. */
107 clipped->length = n;
108 n = ptwXY1->allocatedSize;
109 ptwXY1->allocatedSize = clipped->allocatedSize;
110 clipped->allocatedSize = n;
111 points = clipped->points;
112 clipped->points = ptwXY1->points;
113 ptwXY1->points = points;
114 ptwXY_free( clipped );
115 }
116
117 return( ptwXY1->status );
118
119 Err:
120 ptwXY_free( clipped );
121 return( ptwXY1->status = status );
122 }
123 /*
124 ************************************************************
125 */
126 static nfu_status ptwXY_clip2( ptwXYPoints *clipped, double y, double x1, double y1, double x2, double y2 ) {
127
128 double x;
129 nfu_status status = nfu_Okay;
130
131 x = ( y - y1 ) * ( x2 - x1 ) / ( y2 - y1 ) + x1;
132 if( x <= x1 ) {
133 x = x1; }
134 else if( x >= x2 ) {
135 x = x1; }
136 else {
137 status = ptwXY_setValueAtX( clipped, x, y );
138 }
139 return( status );
140 }
141 /*
142 ************************************************************
143 */
144 nfu_status ptwXY_thicken( ptwXYPoints *ptwXY1, int sectionSubdivideMax, double dxMax, double fxMax ) {
145
146 double x1, x2 = 0., y1, y2 = 0., fx = 1.1, x, dx, dxp, lfx, y; /* fx initialized so compilers want complain. */
147 int64_t i, notFirstPass = 0;
148 int nfx, nDone, doLinear;
149 nfu_status status;
150
151 if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
152 if( ( sectionSubdivideMax < 1 ) || ( dxMax < 0. ) || ( fxMax < 1. ) ) return( nfu_badInput );
153 if( sectionSubdivideMax > ptwXY_sectionSubdivideMax ) sectionSubdivideMax = ptwXY_sectionSubdivideMax;
154 if( ( status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( status );
155 for( i = ptwXY1->length - 1; i >= 0; i-- ) {
156 x1 = ptwXY1->points[i].x;
157 y1 = ptwXY1->points[i].y;
158 if( notFirstPass ) {
159 dx = ptwXY_thicken_linear_dx( sectionSubdivideMax, dxMax, x1, x2 );
160
161 if( x1 == 0. ) {
162 doLinear = 1; }
163 else {
164 fx = x2 / x1;
165 if( fx > 0. ) {
166 lfx = G4Log( fx );
167 if( fxMax == 1. ) {
168 nfx = sectionSubdivideMax; }
169 else {
170 nfx = ( (int) ( lfx / G4Log( fxMax ) ) ) + 1;
171 if( nfx > sectionSubdivideMax ) nfx = sectionSubdivideMax;
172 }
173 if( nfx > 0 ) fx = G4Exp( lfx / nfx );
174 doLinear = 0;
175 if( dx < ( fx - 1 ) * x1 ) doLinear = 1; }
176 else {
177 doLinear = 1;
178 }
179 }
180 x = x1;
181 dxp = dx;
182 nDone = 0;
183 while( 1 ) {
184 if( doLinear ) {
185 x += dx; }
186 else {
187 dx = ptwXY_thicken_linear_dx( sectionSubdivideMax - nDone, dxMax, x, x2 );
188 if( dx <= ( fx - 1 ) * x ) {
189 dxp = dx;
190 doLinear = 1;
191 continue;
192 }
193 dxp = ( fx - 1. ) * x;
194 x *= fx;
195 }
196 if( ( x2 - x ) < 0.05 * std::fabs( dxp ) ) break;
197 if( ( status = ptwXY_interpolatePoint( ptwXY1->interpolation, x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( status );
198 if( ( status = ptwXY_setValueAtX( ptwXY1, x, y ) ) != nfu_Okay ) return( status );
199 nDone++;
200 } // Loop checking, 11.06.2015, T. Koi
201 }
202 notFirstPass = 1;
203 x2 = x1;
204 y2 = y1;
205 }
206 return( status );
207 }
208 /*
209 ************************************************************
210 */
211 ptwXYPoints *ptwXY_thin( ptwXYPoints *ptwXY1, double accuracy, nfu_status *status ) {
212
213 int64_t i, j, length = ptwXY1->length;
214 ptwXYPoints *thinned = NULL;
215 double y1, y2, y3;
216 char *thin = NULL;
217
218 if( length < 3 ) return( ptwXY_clone( ptwXY1, status ) ); /* Logic below requires at least 2 points. */
219 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
220 *status = nfu_otherInterpolation;
221 if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
222 if( accuracy < ptwXY1->accuracy ) accuracy = ptwXY1->accuracy;
223 if( ( thinned = ptwXY_new( ptwXY1->interpolation, &(ptwXY1->interpolationOtherInfo),
224 ptwXY1->biSectionMax, accuracy, length, ptwXY1->overflowLength, status, ptwXY1->userFlag ) ) == NULL ) return( NULL );
225
226 thinned->points[0] = ptwXY1->points[0]; /* This sections removes middle point if surrounding points have the same y-value. */
227 y1 = ptwXY1->points[0].y;
228 y2 = ptwXY1->points[1].y;
229 for( i = 2, j = 1; i < length; i++ ) {
230 y3 = ptwXY1->points[i].y;
231 if( ( y1 != y2 ) || ( y2 != y3 ) ) {
232 thinned->points[j++] = ptwXY1->points[i - 1];
233 y1 = y2;
234 y2 = y3;
235 }
236 }
237 thinned->points[j++] = ptwXY1->points[length - 1];
238
239 if( ptwXY1->interpolation != ptwXY_interpolationFlat ) { /* Now call ptwXY_thin2 for more thinning. */
240 length = thinned->length = j;
241 if( ( thin = (char *) nfu_calloc( 1, (size_t) length ) ) == NULL ) goto Err;
242 if( ( *status = ptwXY_thin2( thinned, thin, accuracy, 0, length - 1 ) ) != nfu_Okay ) goto Err;
243 for( j = 1; j < length; j++ ) if( thin[j] != 0 ) break;
244 for( i = j + 1; i < length; i++ ) {
245 if( thin[i] == 0 ) {
246 thinned->points[j] = thinned->points[i];
247 j++;
248 }
249 }
250 nfu_free( thin );
251 }
252 thinned->length = j;
253
254 return( thinned );
255
256 Err:
257 ptwXY_free( thinned );
258 if( thin != NULL ) nfu_free( thin );
259 return( NULL );
260 }
261 /*
262 ************************************************************
263 */
264 static nfu_status ptwXY_thin2( ptwXYPoints *thinned, char *thin, double accuracy, int64_t i1, int64_t i2 ) {
265
266 int64_t i, iMax = 0;
267 double y, s, dY, dYMax = 0., dYR, dYRMax = 0;
268 double x1 = thinned->points[i1].x, y1 = thinned->points[i1].y, x2 = thinned->points[i2].x, y2 = thinned->points[i2].y;
269 nfu_status status = nfu_Okay;
270
271 if( i1 + 1 >= i2 ) return( nfu_Okay );
272 for( i = i1 + 1; i < i2; i++ ) {
273 if( ( status = ptwXY_interpolatePoint( thinned->interpolation, thinned->points[i].x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( status );
274 s = 0.5 * ( std::fabs( y ) + std::fabs( thinned->points[i].y ) );
275 dY = std::fabs( y - thinned->points[i].y );
276 dYR = 0;
277 if( s != 0 ) dYR = dY / s;
278 if( ( dYR > dYRMax ) || ( ( dYR >= 0.9999 * dYRMax ) && ( dY > dYMax ) ) ) { /* The choice of 0.9999 is not exact science. */
279 iMax = i;
280 if( dY > dYMax ) dYMax = dY;
281 if( dYR > dYRMax ) dYRMax = dYR;
282 }
283 }
284 if( dYRMax < accuracy ) {
285 for( i = i1 + 1; i < i2; i++ ) thin[i] = 1; }
286 else {
287 if( ( status = ptwXY_thin2( thinned, thin, accuracy, i1, iMax ) ) != nfu_Okay ) return( status );
288 status = ptwXY_thin2( thinned, thin, accuracy, iMax, i2 );
289 }
290 return( status );
291 }
292 /*
293 ************************************************************
294 */
295 static double ptwXY_thicken_linear_dx( int sectionSubdivideMax, double dxMax, double x1, double x2 ) {
296
297 int ndx;
298 double dx = x2 - x1, dndx;
299
300 if( dxMax == 0. ) {
301 dx = ( x2 - x1 ) / sectionSubdivideMax; }
302 else {
303 dndx = dx / dxMax;
304 ndx = (int) dndx;
305 if( ( dndx - ndx ) > 1e-6 ) ndx++;
306 if( ndx > sectionSubdivideMax ) ndx = sectionSubdivideMax;
307 if( ndx > 0 ) dx /= ndx;
308 }
309
310 return( dx );
311 }
312 /*
313 ************************************************************
314 */
315 nfu_status ptwXY_trim( ptwXYPoints *ptwXY ) {
316 /*
317 c Remove extra zeros at beginning and end.
318 */
319
320 int64_t i, i1, i2;
321 nfu_status status;
322
323 if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
324 if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status );
325 for( i1 = 0; i1 < ptwXY->length; i1++ ) {
326 if( ptwXY->points[i1].y != 0 ) break;
327 }
328 if( i1 > 0 ) i1--;
329 for( i2 = ptwXY->length - 1; i2 >= 0; i2-- ) {
330 if( ptwXY->points[i2].y != 0 ) break;
331 }
332 i2++;
333 if( i2 < ptwXY->length ) i2++;
334 if( i2 > i1 ) {
335 if( i1 > 0 ) {
336 for( i = i1; i < i2; i++ ) ptwXY->points[i - i1] = ptwXY->points[i];
337 }
338 ptwXY->length = i2 - i1; }
339 else if( i2 < i1 ) { /* Remove all zeros between endpoints. */
340 ptwXY->points[1] = ptwXY->points[ptwXY->length - 1];
341 ptwXY->length = 2;
342 }
343
344 return( nfu_Okay );
345 }
346 /*
347 ************************************************************
348 */
349 ptwXYPoints *ptwXY_union( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int unionOptions ) {
350
351 int64_t overflowSize, i, i1 = 0, i2 = 0, n1 = ptwXY1->length, n2 = ptwXY2->length, length;
352 int fillWithFirst = unionOptions & ptwXY_union_fill, trim = unionOptions & ptwXY_union_trim;
353 ptwXYPoints *n;
354 double x1 = 0., x2 = 0., y1 = 0., y2 = 0., y, biSectionMax, accuracy;
355
356 if( ( *status = ptwXY1->status ) != nfu_Okay ) return( NULL );
357 if( ( *status = ptwXY2->status ) != nfu_Okay ) return( NULL );
358 *status = nfu_otherInterpolation;
359 if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
360 /*
361 * Many other routines use the fact that ptwXY_union calls ptwXY_coalescePoints for ptwXY1 and ptwXY2 so do not change it.
362 */
363 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
364 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY2 ) ) != nfu_Okay ) return( NULL );
365
366 if( ( n1 == 1 ) || ( n2 == 1 ) ) {
367 *status = nfu_tooFewPoints;
368 return( NULL );
369 }
370 if( trim ) {
371 if( n1 > 0 ) {
372 if( n2 > 0 ) {
373 if( ptwXY1->points[0].x < ptwXY2->points[0].x ) {
374 while( i1 < n1 ) { // Loop checking, 11.05.2015, T. Koi
375 if( ptwXY1->points[i1].x >= ptwXY2->points[0].x ) break;
376 if( fillWithFirst ) {
377 if( i1 < ( ptwXY1->length - 1 ) ) {
378 x1 = ptwXY1->points[i1].x;
379 y1 = ptwXY1->points[i1].y;
380 x2 = ptwXY1->points[i1+1].x;
381 y2 = ptwXY1->points[i1+1].y;
382 }
383 }
384 i1++;
385 } }
386 else {
387 while( i2 < n2 ) { // Loop checking, 11.06.2015, T. Koi
388 if( ptwXY2->points[i2].x >= ptwXY1->points[0].x ) break;
389 i2++;
390 }
391 }
392 if( ptwXY1->points[n1-1].x > ptwXY2->points[n2-1].x ) {
393 while( i1 < n1 ) { // Loop checking, 11.06.2015, T. Koi
394 if( ptwXY1->points[n1-1].x <= ptwXY2->points[n2-1].x ) break;
395 n1--;
396 } }
397 else {
398 while( i2 < n2 ) { // Loop checking, 11.06.2015, T. Koi
399 if( ptwXY2->points[n2-1].x <= ptwXY1->points[n1-1].x ) break;
400 n2--;
401 }
402 } }
403 else {
404 n1 = 0;
405 } }
406 else {
407 n2 = 0;
408 }
409 }
410 overflowSize = ptwXY1->overflowAllocatedSize;
411 if( overflowSize < ptwXY2->overflowAllocatedSize ) overflowSize = ptwXY2->overflowAllocatedSize;
412 length = ( n1 - i1 ) + ( n2 - i2 );
413 if( length == 0 ) length = ptwXY_minimumSize;
414 biSectionMax = ptwXY1->biSectionMax;
415 if( biSectionMax < ptwXY2->biSectionMax ) biSectionMax = ptwXY2->biSectionMax;
416 accuracy = ptwXY1->accuracy;
417 if( accuracy < ptwXY2->accuracy ) accuracy = ptwXY2->accuracy;
418 n = ptwXY_new( ptwXY1->interpolation, NULL, biSectionMax, accuracy, length, overflowSize, status, ptwXY1->userFlag );
419 if( n == NULL ) return( NULL );
420
421 for( i = 0; ( i1 < n1 ) && ( i2 < n2 ); i++ ) {
422 y = 0.;
423 if( ptwXY1->points[i1].x <= ptwXY2->points[i2].x ) {
424 n->points[i].x = ptwXY1->points[i1].x;
425 if( fillWithFirst ) {
426 y = ptwXY1->points[i1].y;
427 if( i1 < ( ptwXY1->length - 1 ) ) {
428 x1 = ptwXY1->points[i1].x;
429 y1 = ptwXY1->points[i1].y;
430 x2 = ptwXY1->points[i1+1].x;
431 y2 = ptwXY1->points[i1+1].y; }
432 else {
433 y1 = 0.;
434 y2 = 0.;
435 }
436 }
437 if( ptwXY1->points[i1].x == ptwXY2->points[i2].x ) i2++;
438 i1++; }
439 else {
440 n->points[i].x = ptwXY2->points[i2].x;
441 if( fillWithFirst && ( ( y1 != 0. ) || ( y2 != 0. ) ) ) {
442 if( ( *status = ptwXY_interpolatePoint( ptwXY1->interpolation, ptwXY2->points[i2].x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) {
443 ptwXY_free( n );
444 return( NULL );
445 }
446 }
447 i2++;
448 }
449 n->points[i].y = y;
450 }
451
452 y = 0.;
453 for( ; i1 < n1; i1++, i++ ) {
454 n->points[i].x = ptwXY1->points[i1].x;
455 if( fillWithFirst ) y = ptwXY1->points[i1].y;
456 n->points[i].y = y;
457 }
458 for( ; i2 < n2; i2++, i++ ) {
459 n->points[i].x = ptwXY2->points[i2].x;
460 if( fillWithFirst && trim && ( n->points[i].x <= x2 ) ) {
461 if( ( *status = ptwXY_interpolatePoint( ptwXY1->interpolation, n->points[i].x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) {
462 ptwXY_free( n );
463 return( NULL );
464 }
465 }
466 n->points[i].y = y;
467 }
468 n->length = i;
469
470 if( unionOptions & ptwXY_union_mergeClosePoints ) {
471 if( ( *status = ptwXY_mergeClosePoints( n, 4 * DBL_EPSILON ) ) != nfu_Okay ) {
472 ptwXY_free( n );
473 return( NULL );
474 }
475 }
476 return( n );
477 }
478 /*
479 ************************************************************
480 */
481 nfu_status ptwXY_scaleOffsetXAndY( ptwXYPoints *ptwXY, double xScale, double xOffset, double yScale, double yOffset ) {
482
483 int64_t i1, length = ptwXY->length;
484 ptwXYPoint *p1;
485 nfu_status status;
486
487 if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
488 if( xScale == 0 ) return( nfu_XNotAscending );
489
490 if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status );
491
492 for( i1 = 0, p1 = ptwXY->points; i1 < length; i1++, p1++ ) {
493 p1->x = xScale * p1->x + xOffset;
494 p1->y = yScale * p1->y + yOffset;
495 }
496
497 if( xScale < 0 ) {
498 int64_t length_2 = length / 2;
499 ptwXYPoint tmp, *p2;
500
501 for( i1 = 0, p1 = ptwXY->points, p2 = &(ptwXY->points[length-1]); i1 < length_2; i1++ ) {
502 tmp = *p1;
503 *p1 = *p2;
504 *p2 = tmp;
505 }
506 }
507
508 return( ptwXY->status );
509 }
510
511 #if defined __cplusplus
512 }
513 #endif
514