| Geant4 Cross Reference |
1 /* 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 *pt
20 static double ptwXY_thicken_linear_dx( int sec
21 static nfu_status ptwXY_thin2( ptwXYPoints *th
22 /*
23 **********************************************
24 */
25 nfu_status ptwXY_clip( ptwXYPoints *ptwXY1, do
26 /*
27 This function acts oddly for xy = [ [ 1, 0
28 This function probably only works for line
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(
37 if( ptwXY1->interpolation == ptwXY_interpo
38 n = ptwXY1->length;
39 if( n > 0 ) {
40 i = 0;
41 if( ptwXY_getYMax( ptwXY1 ) < yMin ) i
42 if( ptwXY_getYMin( ptwXY1 ) > yMax ) i
43 if( i == 1 ) return( ptwXY_clear( ptwX
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->int
54 ptwXY1->biSectionMax, ptwXY1->
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_getPointAtI
62 if( points->y > yMin ) {
63 if( ( status = ptwXY_c
64 }
65 }
66 if( ( status = ptwXY_setValueA
67 j = i;
68 for( i++; i < n; i++ ) if( !(
69 if( i < n ) {
70 x2 = ptwXY1->points[i].x;
71 y2 = ptwXY1->points[i].y;
72 if( ( status = ptwXY_clip2
73 if( y2 > yMax ) {
74 if( ( status = ptwXY_c
75 } }
76 else if( j != n - 1 ) {
77 if( ( status = ptwXY_setVa
78 }
79 i--; }
80 else if( y2 > yMax ) {
81 if( i > 0 ) {
82 points = ptwXY_getPointAtI
83 if( points->y < yMax ) {
84 if( ( status = ptwXY_c
85 }
86 }
87 if( ( status = ptwXY_setValueA
88 j = i;
89 for( i++; i < n; i++ ) if( !(
90 if( i < n ) {
91 x2 = ptwXY1->points[i].x;
92 y2 = ptwXY1->points[i].y;
93 if( ( status = ptwXY_clip2
94 if( y2 < yMin ) {
95 if( ( status = ptwXY_c
96 } }
97 else if( j != n - 1 ) {
98 if( ( status = ptwXY_setVa
99 }
100 i--; }
101 else {
102 if( ( status = ptwXY_setValueA
103 }
104 }
105 if( ( status = ptwXY_simpleCoalescePoi
106 ptwXY1->length = clipped->length; /*
107 clipped->length = n;
108 n = ptwXY1->allocatedSize;
109 ptwXY1->allocatedSize = clipped->alloc
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 *cl
127
128 double x;
129 nfu_status status = nfu_Okay;
130
131 x = ( y - y1 ) * ( x2 - x1 ) / ( y2 - y1 )
132 if( x <= x1 ) {
133 x = x1; }
134 else if( x >= x2 ) {
135 x = x1; }
136 else {
137 status = ptwXY_setValueAtX( clipped, x
138 }
139 return( status );
140 }
141 /*
142 **********************************************
143 */
144 nfu_status ptwXY_thicken( ptwXYPoints *ptwXY1,
145
146 double x1, x2 = 0., y1, y2 = 0., fx = 1.1,
147 int64_t i, notFirstPass = 0;
148 int nfx, nDone, doLinear;
149 nfu_status status;
150
151 if( ptwXY1->interpolation == ptwXY_interpo
152 if( ( sectionSubdivideMax < 1 ) || ( dxMax
153 if( sectionSubdivideMax > ptwXY_sectionSub
154 if( ( status = ptwXY_simpleCoalescePoints(
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( sect
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 = sectionSubdivide
169 else {
170 nfx = ( (int) ( lfx /
171 if( nfx > sectionSubdi
172 }
173 if( nfx > 0 ) fx = G4Exp(
174 doLinear = 0;
175 if( dx < ( fx - 1 ) * x1 )
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_
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::f
197 if( ( status = ptwXY_interpola
198 if( ( status = ptwXY_setValueA
199 nDone++;
200 } // Loop checking, 11.06.2015, T.
201 }
202 notFirstPass = 1;
203 x2 = x1;
204 y2 = y1;
205 }
206 return( status );
207 }
208 /*
209 **********************************************
210 */
211 ptwXYPoints *ptwXY_thin( ptwXYPoints *ptwXY1,
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( ptwX
219 if( ( *status = ptwXY_simpleCoalescePoints
220 *status = nfu_otherInterpolation;
221 if( ptwXY1->interpolation == ptwXY_interpo
222 if( accuracy < ptwXY1->accuracy ) accuracy
223 if( ( thinned = ptwXY_new( ptwXY1->interpo
224 ptwXY1->biSectionMax, accuracy, length
225
226 thinned->points[0] = ptwXY1->points[0];
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->poi
233 y1 = y2;
234 y2 = y3;
235 }
236 }
237 thinned->points[j++] = ptwXY1->points[leng
238
239 if( ptwXY1->interpolation != ptwXY_interpo
240 length = thinned->length = j;
241 if( ( thin = (char *) nfu_calloc( 1, (
242 if( ( *status = ptwXY_thin2( thinned,
243 for( j = 1; j < length; j++ ) if( thin
244 for( i = j + 1; i < length; i++ ) {
245 if( thin[i] == 0 ) {
246 thinned->points[j] = thinned->
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 *th
265
266 int64_t i, iMax = 0;
267 double y, s, dY, dYMax = 0., dYR, dYRMax =
268 double x1 = thinned->points[i1].x, y1 = th
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(
274 s = 0.5 * ( std::fabs( y ) + std::fabs
275 dY = std::fabs( y - thinned->points[i]
276 dYR = 0;
277 if( s != 0 ) dYR = dY / s;
278 if( ( dYR > dYRMax ) || ( ( dYR >= 0.9
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]
286 else {
287 if( ( status = ptwXY_thin2( thinned, t
288 status = ptwXY_thin2( thinned, thin, a
289 }
290 return( status );
291 }
292 /*
293 **********************************************
294 */
295 static double ptwXY_thicken_linear_dx( int sec
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 =
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( pt
324 if( ( status = ptwXY_simpleCoalescePoints(
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->
337 }
338 ptwXY->length = i2 - i1; }
339 else if( i2 < i1 ) { /* Rem
340 ptwXY->points[1] = ptwXY->points[ptwXY
341 ptwXY->length = 2;
342 }
343
344 return( nfu_Okay );
345 }
346 /*
347 **********************************************
348 */
349 ptwXYPoints *ptwXY_union( ptwXYPoints *ptwXY1,
350
351 int64_t overflowSize, i, i1 = 0, i2 = 0, n
352 int fillWithFirst = unionOptions & ptwXY_u
353 ptwXYPoints *n;
354 double x1 = 0., x2 = 0., y1 = 0., y2 = 0.,
355
356 if( ( *status = ptwXY1->status ) != nfu_Ok
357 if( ( *status = ptwXY2->status ) != nfu_Ok
358 *status = nfu_otherInterpolation;
359 if( ptwXY1->interpolation == ptwXY_interpo
360 /*
361 * Many other routines use the fact that ptwX
362 */
363 if( ( *status = ptwXY_simpleCoalescePoints
364 if( ( *status = ptwXY_simpleCoalescePoints
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 < ptwX
374 while( i1 < n1 ) { // Loop
375 if( ptwXY1->points[i1]
376 if( fillWithFirst ) {
377 if( i1 < ( ptwXY1-
378 x1 = ptwXY1->p
379 y1 = ptwXY1->p
380 x2 = ptwXY1->p
381 y2 = ptwXY1->p
382 }
383 }
384 i1++;
385 } }
386 else {
387 while( i2 < n2 ) { // Loop
388 if( ptwXY2->points[i2]
389 i2++;
390 }
391 }
392 if( ptwXY1->points[n1-1].x > p
393 while( i1 < n1 ) { // Loop
394 if( ptwXY1->points[n1-
395 n1--;
396 } }
397 else {
398 while( i2 < n2 ) { // Loop
399 if( ptwXY2->points[n2-
400 n2--;
401 }
402 } }
403 else {
404 n1 = 0;
405 } }
406 else {
407 n2 = 0;
408 }
409 }
410 overflowSize = ptwXY1->overflowAllocatedSi
411 if( overflowSize < ptwXY2->overflowAllocat
412 length = ( n1 - i1 ) + ( n2 - i2 );
413 if( length == 0 ) length = ptwXY_minimumSi
414 biSectionMax = ptwXY1->biSectionMax;
415 if( biSectionMax < ptwXY2->biSectionMax )
416 accuracy = ptwXY1->accuracy;
417 if( accuracy < ptwXY2->accuracy ) accuracy
418 n = ptwXY_new( ptwXY1->interpolation, NULL
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->po
424 n->points[i].x = ptwXY1->points[i1
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].
431 y2 = ptwXY1->points[i1+1].
432 else {
433 y1 = 0.;
434 y2 = 0.;
435 }
436 }
437 if( ptwXY1->points[i1].x == ptwXY2
438 i1++; }
439 else {
440 n->points[i].x = ptwXY2->points[i2
441 if( fillWithFirst && ( ( y1 != 0.
442 if( ( *status = ptwXY_interpol
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
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->poin
461 if( ( *status = ptwXY_interpolateP
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_mergeCloseP
471 if( ( *status = ptwXY_mergeClosePoints
472 ptwXY_free( n );
473 return( NULL );
474 }
475 }
476 return( n );
477 }
478 /*
479 **********************************************
480 */
481 nfu_status ptwXY_scaleOffsetXAndY( ptwXYPoints
482
483 int64_t i1, length = ptwXY->length;
484 ptwXYPoint *p1;
485 nfu_status status;
486
487 if( ptwXY->status != nfu_Okay ) return( pt
488 if( xScale == 0 ) return( nfu_XNotAscendin
489
490 if( ( status = ptwXY_simpleCoalescePoints(
491
492 for( i1 = 0, p1 = ptwXY->points; i1 < leng
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 =
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