| Geant4 Cross Reference |
1 /* 1
2 # <<BEGIN-copyright>>
3 # <<END-copyright>>
4 */
5
6 #include <stdlib.h>
7 #include <cmath>
8 #include <float.h>
9
10 #include "ptwXY.h"
11
12 #if defined __cplusplus
13 #include <cmath>
14 #include "G4Exp.hh"
15 #include "G4Log.hh"
16 namespace GIDI {
17 using namespace GIDI;
18 #endif
19
20 static nfu_status ptwXY_createGaussianCentered
21 /*
22 **********************************************
23 */
24 ptwXPoints *ptwXY_getXArray( ptwXYPoints *ptwX
25
26 int64_t i, n;
27 ptwXPoints *xArray;
28
29 if( ( *status = ptwXY->status ) != nfu_Oka
30 n = ptwXY->length;
31
32 if( ( *status = ptwXY_simpleCoalescePoints
33 if( ( xArray = ptwX_new( n, status ) ) ==
34 for( i = 0; i < n; i++ ) xArray->points[i]
35 xArray->length = n;
36
37 return( xArray );
38 }
39 /*
40 **********************************************
41 */
42 nfu_status ptwXY_dullEdges( ptwXYPoints *ptwXY
43
44 #define minEps 5e-16
45
46 nfu_status status;
47 double xm, xp, dx, y, x1, y1, x2, y2, sign
48 ptwXYPoint *p;
49
50 /* This routine can only be used for linear in
51 This needs to be fixed and documented. */
52 if( ( status = ptwXY->status ) != nfu_Okay
53 if( ptwXY->interpolation == ptwXY_interpol
54 if( ptwXY->interpolation == ptwXY_interpol
55
56 if( ptwXY->length < 2 ) return( nfu_Okay )
57
58 if( lowerEps != 0. ) {
59 if( std::fabs( lowerEps ) < minEps ) {
60 sign = 1;
61 if( lowerEps < 0. ) sign = -1;
62 lowerEps = sign * minEps;
63 }
64
65 p = ptwXY_getPointAtIndex_Unsafely( pt
66 x1 = p->x;
67 y1 = p->y;
68 p = ptwXY_getPointAtIndex_Unsafely( pt
69 x2 = p->x;
70 y2 = p->y;
71
72 if( y1 != 0. ) {
73 dx = std::fabs( x1 * lowerEps );
74 if( x1 == 0 ) dx = std::fabs( lowe
75 xm = x1 - dx;
76 xp = x1 + dx;
77 if( ( xp + dx ) < x2 ) {
78 if( ( status = ptwXY_getValueA
79 if( ( status = ptwXY_setValueA
80 else {
81 xp = x2;
82 y = y2;
83 }
84 if( lowerEps > 0 ) {
85 if( ( status = ptwXY_setValueA
86 else {
87 if( ( xm < 0. ) && ( x1 >= 0.
88 if( ( status = ptwXY_setVa
89 else {
90 if( ( status = ptwXY_setVa
91 if( ( status = ptwXY_inter
92 if( ( status = ptwXY_setVa
93 }
94 }
95 }
96 }
97
98 if( upperEps != 0. ) {
99 if( std::fabs( upperEps ) < minEps ) {
100 sign = 1;
101 if( upperEps < 0. ) sign = -1;
102 upperEps = sign * minEps;
103 }
104
105 p = ptwXY_getPointAtIndex_Unsafely( pt
106 x1 = p->x;
107 y1 = p->y;
108 p = ptwXY_getPointAtIndex_Unsafely( pt
109 x2 = p->x;
110 y2 = p->y;
111
112 if( y2 != 0. ) {
113 dx = std::fabs( x2 * upperEps );
114 if( x2 == 0 ) dx = std::fabs( uppe
115 xm = x2 - dx;
116 xp = x2 + dx;
117 if( ( xm - dx ) > x1 ) {
118 if( ( status = ptwXY_getValueA
119 if( ( status = ptwXY_setValueA
120 else {
121 xm = x1;
122 y = y1;
123 }
124 if( upperEps < 0 ) {
125 if( ( status = ptwXY_setValueA
126 else {
127 if( ( status = ptwXY_setValueA
128 if( ( status = ptwXY_interpola
129 if( ( status = ptwXY_setValueA
130 }
131 }
132 }
133
134 return( ptwXY->status );
135
136 #undef minEps
137 }
138 /*
139 **********************************************
140 */
141 nfu_status ptwXY_mergeClosePoints( ptwXYPoints
142
143 int64_t i, i1, j, k, n = ptwXY->length;
144 double x, y;
145 ptwXYPoint *p1, *p2;
146
147 if( n < 2 ) return( ptwXY->status );
148 if( epsilon < 4 * DBL_EPSILON ) epsilon =
149 if( ptwXY_simpleCoalescePoints( ptwXY ) !=
150
151 p2 = ptwXY->points;
152 x = p2->x;
153 for( i1 = 1, p2++; i1 < ( n - 1 ); i1++, p
154 if( ( p2->x - x ) > 0.5 * epsilon * (
155 }
156 if( i1 != 1 ) {
157 for( i = i1, p1 = &(ptwXY->points[1]);
158 n = ptwXY->length = ptwXY->length - i1
159 }
160
161 p1 = &(ptwXY->points[n-1]);
162 x = p1->x;
163 for( i1 = n - 2, p1--; i1 > 0; i1--, p1--
164 if( x - p1->x > 0.5 * epsilon * ( std:
165 }
166 if( i1 != ( n - 2 ) ) {
167 ptwXY->points[i1 + 1] = ptwXY->points[
168 n = ptwXY->length = i1 + 2;
169 }
170
171 for( i = 1; i < n - 1; i++ ) {
172 p1 = &(ptwXY->points[i]);
173 x = p1->x;
174 y = p1->y;
175 for( j = i + 1, p2 = &(ptwXY->points[i
176 if( ( p2->x - p1->x ) > 0.5 * epsi
177 x += p2->x;
178 y += p2->y;
179 }
180 if( ( k = ( j - i ) ) > 1 ) {
181 p1->x = x / k;
182 p1->y = y / k;
183 for( p1 = &(ptwXY->points[i+1]); j
184 n -= ( k - 1 );
185 }
186 }
187 ptwXY->length = n;
188
189 return( ptwXY->status );
190 }
191 /*
192 **********************************************
193 */
194 ptwXYPoints *ptwXY_intersectionWith_ptwX( ptwX
195
196 int64_t i, i1, i2, lengthX = ptwX_length(
197 double x, y, xMin, xMax;
198 ptwXYPoints *n = NULL;
199
200 if( ( *status = ptwXY->status ) != nfu_Oka
201 if( ( *status = ptwX->status ) != nfu_Okay
202 if( ( *status = ptwXY_simpleCoalescePoints
203 *status = nfu_otherInterpolation;
204 if( ptwXY->interpolation == ptwXY_interpol
205 if( ( n = ptwXY_clone( ptwXY, status ) ) =
206 if( ptwXY->length == 0 ) return( n );
207 xMin = ptwXY->points[0].x;
208 xMax = ptwXY->points[ptwXY->length - 1].x;
209
210 if( ( xMin >= ptwX->points[lengthX-1] ) ||
211 n->length = 0;
212 return( n );
213 }
214
215 for( i = 0; i < lengthX; i++ ) { /*
216 x = ptwX->points[i];
217 if( x <= xMin ) continue;
218 if( x >= xMax ) break;
219 if( ( *status = ptwXY_getValueAtX( ptw
220 if( ( *status = ptwXY_setValueAtX( n,
221 }
222 if( ( *status = ptwXY_simpleCoalescePoints
223
224 i1 = 0;
225 i2 = n->length - 1;
226 if( lengthX > 0 ) {
227 x = ptwX->points[0];
228 if( x > n->points[i1].x ) {
229 for( ; i1 < n->length; i1++ ) {
230 if( n->points[i1].x == x ) bre
231 }
232 }
233
234 x = ptwX->points[lengthX - 1];
235 if( x < n->points[i2].x ) {
236 for( ; i2 > i1; i2-- ) {
237 if( n->points[i2].x == x ) bre
238 }
239 }
240 }
241 i2++;
242
243 if( i1 != 0 ) {
244 for( i = i1; i < i2; i++ ) n->points[i
245 }
246 n->length = i2 - i1;
247
248 return( n );
249
250 Err:
251 ptwXY_free( n );
252 return( NULL );
253 }
254 /*
255 **********************************************
256 */
257 nfu_status ptwXY_areDomainsMutual( ptwXYPoints
258
259 nfu_status status;
260 int64_t n1 = ptwXY1->length, n2 = ptwXY2->
261 ptwXYPoint *xy1, *xy2;
262
263 if( ( status = ptwXY1->status ) != nfu_Oka
264 if( ( status = ptwXY2->status ) != nfu_Oka
265 if( n1 == 0 ) return( nfu_empty );
266 if( n2 == 0 ) return( nfu_empty );
267 if( n1 < 2 ) {
268 status = nfu_tooFewPoints; }
269 else if( n2 < 2 ) {
270 status = nfu_tooFewPoints; }
271 else {
272 xy1 = ptwXY_getPointAtIndex_Unsafely(
273 xy2 = ptwXY_getPointAtIndex_Unsafely(
274 if( xy1->x < xy2->x ) {
275 if( xy2->y != 0. ) status = nfu_do
276 else if( xy1->x > xy2->x ) {
277 if( xy1->y != 0. ) status = nfu_do
278 }
279
280 if( status == nfu_Okay ) {
281 xy1 = ptwXY_getPointAtIndex_Unsafe
282 xy2 = ptwXY_getPointAtIndex_Unsafe
283 if( xy1->x < xy2->x ) {
284 if( xy1->y != 0. ) status = nf
285 else if( xy1->x > xy2->x ) {
286 if( xy2->y != 0. ) status = nf
287 }
288 }
289 }
290 return( status );
291 }
292 /*
293 **********************************************
294 */
295 nfu_status ptwXY_tweakDomainsToMutualify( ptwX
296
297 nfu_status status;
298 int64_t n1 = ptwXY1->length, n2 = ptwXY2->
299 double sum, diff;
300 ptwXYPoint *xy1, *xy2;
301
302 epsilon = std::fabs( epsilon ) + std::fabs
303
304 if( ( status = ptwXY1->status ) != nfu_Oka
305 if( ( status = ptwXY2->status ) != nfu_Oka
306 if( n1 == 0 ) return( nfu_empty );
307 if( n2 == 0 ) return( nfu_empty );
308 if( n1 < 2 ) {
309 status = nfu_tooFewPoints; }
310 else if( n2 < 2 ) {
311 status = nfu_tooFewPoints; }
312 else {
313 xy1 = ptwXY_getPointAtIndex_Unsafely(
314 xy2 = ptwXY_getPointAtIndex_Unsafely(
315 if( xy1->x < xy2->x ) {
316 if( xy2->y != 0. ) {
317 sum = std::fabs( xy1->x ) + st
318 diff = std::fabs( xy2->x - xy1
319 if( diff > epsilon * sum ) {
320 status = nfu_domainsNotMut
321 else {
322 xy1->x = xy2->x;
323 }
324 } }
325 else if( xy1->x > xy2->x ) {
326 if( xy1->y != 0. ) {
327 sum = std::fabs( xy1->x ) + st
328 diff = std::fabs( xy2->x - xy1
329 if( diff > epsilon * sum ) {
330 status = nfu_domainsNotMut
331 else {
332 xy2->x = xy1->x;
333 }
334 }
335 }
336
337 if( status == nfu_Okay ) {
338 xy1 = ptwXY_getPointAtIndex_Unsafe
339 xy2 = ptwXY_getPointAtIndex_Unsafe
340 if( xy1->x < xy2->x ) {
341 if( xy1->y != 0. ) {
342 sum = std::fabs( xy1->x )
343 diff = std::fabs( xy2->x -
344 if( diff > epsilon * sum )
345 status = nfu_domainsNo
346 else {
347 xy2->x = xy1->x;
348 }
349 } }
350 else if( xy1->x > xy2->x ) {
351 if( xy2->y != 0. ) {
352 sum = std::fabs( xy1->x )
353 diff = std::fabs( xy2->x -
354 if( diff > epsilon * sum )
355 status = nfu_domainsNo
356 else {
357 xy1->x = xy2->x;
358 }
359 }
360 }
361 }
362 }
363 return( status );
364 }
365 /*
366 **********************************************
367 */
368 nfu_status ptwXY_mutualifyDomains( ptwXYPoints
369 ptwX
370
371 nfu_status status;
372 int64_t n1 = ptwXY1->length, n2 = ptwXY2->
373 ptwXYPoint *xy1, *xy2;
374
375 switch( status = ptwXY_areDomainsMutual( p
376 case nfu_Okay :
377 case nfu_empty :
378 return( nfu_Okay );
379 case nfu_domainsNotMutual :
380 break;
381 default :
382 return( status );
383 }
384
385 if( ptwXY1->interpolation == ptwXY_interpo
386 if( ptwXY2->interpolation == ptwXY_interpo
387 if( ptwXY1->interpolation == ptwXY_interpo
388 if( ptwXY2->interpolation == ptwXY_interpo
389
390 xy1 = ptwXY_getPointAtIndex_Unsafely( ptwX
391 xy2 = ptwXY_getPointAtIndex_Unsafely( ptwX
392 if( xy1->x < xy2->x ) {
393 lowerEps1 = 0.;
394 if( xy2->y == 0. ) lowerEps2 = 0.; }
395 else if( xy1->x > xy2->x ) {
396 lowerEps2 = 0.;
397 if( xy1->y == 0. ) lowerEps1 = 0.; }
398 else {
399 lowerEps1 = lowerEps2 = 0.;
400 }
401
402 xy1 = ptwXY_getPointAtIndex_Unsafely( ptwX
403 xy2 = ptwXY_getPointAtIndex_Unsafely( ptwX
404 if( xy1->x < xy2->x ) {
405 upperEps2 = 0.;
406 if( xy1->y == 0. ) upperEps1 = 0.; }
407 else if( xy1->x > xy2->x ) {
408 upperEps1 = 0.;
409 if( xy2->y == 0. ) upperEps2 = 0.; }
410 else {
411 upperEps1 = upperEps2 = 0.;
412 }
413
414 if( ( lowerEps1 != 0. ) || ( upperEps1 !=
415 if( ( status = ptwXY_dullEdges( ptwXY1
416 if( ( lowerEps2 != 0. ) || ( upperEps2 !=
417 if( ( status = ptwXY_dullEdges( ptwXY2
418
419 return( status );
420 }
421 /*
422 **********************************************
423 */
424 nfu_status ptwXY_copyToC_XY( ptwXYPoints *ptwX
425
426 int64_t i;
427 double *d = xys;
428 nfu_status status;
429 ptwXYPoint *pointFrom;
430
431 if( ptwXY->status != nfu_Okay ) return( pt
432 if( ( status = ptwXY_simpleCoalescePoints(
433
434 if( index1 < 0 ) index1 = 0;
435 if( index2 > ptwXY->length ) index2 = ptwX
436 if( index2 < index1 ) index2 = index1;
437 *numberOfPoints = index2 - index1;
438 if( allocatedSize < ( index2 - index1 ) )
439
440 for( i = index1, pointFrom = ptwXY->points
441 *(d++) = pointFrom->x;
442 *(d++) = pointFrom->y;
443 }
444
445 return( nfu_Okay );
446 }
447 /*
448 **********************************************
449 */
450 nfu_status ptwXY_valueTo_ptwXAndY( ptwXYPoints
451
452 int64_t i1, length = ptwXY_length( ptwXY )
453 double *xps, *yps;
454 ptwXYPoint *pointFrom;
455 nfu_status status;
456
457 if( ptwXY->status != nfu_Okay ) return( pt
458 if( ( status = ptwXY_simpleCoalescePoints(
459
460 if( ( *xs = (double *) malloc( length * si
461 if( ( *ys = (double *) malloc( length * si
462 free( *xs );
463 *xs = NULL;
464 return( nfu_mallocError );
465 }
466
467 for( i1 = 0, pointFrom = ptwXY->points, xp
468 *xps = pointFrom->x;
469 *yps = pointFrom->y;
470 }
471
472 return( nfu_Okay );
473 }
474 /*
475 **********************************************
476 */
477 ptwXYPoints *ptwXY_valueTo_ptwXY( double x1, d
478
479 ptwXYPoints *n;
480
481 *status = nfu_XNotAscending;
482 if( x1 >= x2 ) return( NULL );
483 *status = nfu_Okay;
484 if( ( n = ptwXY_new( ptwXY_interpolationLi
485 ptwXY_setValueAtX( n, x1, y );
486 ptwXY_setValueAtX( n, x2, y );
487 return( n );
488 }
489 /*
490 **********************************************
491 */
492 ptwXYPoints *ptwXY_createGaussianCenteredSigma
493
494 int64_t i, n;
495 ptwXYPoint *pm, *pp;
496 double x1, y1, x2, y2, accuracy2, yMin = 1
497 ptwXYPoints *gaussian;
498
499 if( accuracy < 1e-5 ) accuracy = 1e-5;
500 if( accuracy > 1e-1 ) accuracy = 1e-1;
501 if( ( gaussian = ptwXY_new( ptwXY_interpol
502 accuracy2 = accuracy = gaussian->accuracy;
503 if( accuracy2 > 5e-3 ) accuracy2 = 5e-3;
504
505 x1 = -std::sqrt( -2. * G4Log( yMin ) );
506 y1 = yMin;
507 x2 = -5.2;
508 y2 = G4Exp( -0.5 * x2 * x2 );
509 if( ( *status = ptwXY_setValueAtX( gaussia
510 gaussian->accuracy = 20 * accuracy2;
511 if( ( *status = ptwXY_createGaussianCenter
512 x1 = x2;
513 y1 = y2;
514 x2 = -4.;
515 y2 = G4Exp( -0.5 * x2 * x2 );
516 gaussian->accuracy = 5 * accuracy2;
517 if( ( *status = ptwXY_createGaussianCenter
518 x1 = x2;
519 y1 = y2;
520 x2 = -1;
521 y2 = G4Exp( -0.5 * x2 * x2 );
522 gaussian->accuracy = accuracy;
523 if( ( *status = ptwXY_createGaussianCenter
524 x1 = x2;
525 y1 = y2;
526 x2 = 0;
527 y2 = G4Exp( -0.5 * x2 * x2 );
528 if( ( *status = ptwXY_createGaussianCenter
529
530 n = gaussian->length;
531 if( ( *status = ptwXY_coalescePoints( gaus
532 if( ( *status = ptwXY_setValueAtX( gaussia
533 pp = &(gaussian->points[gaussian->length])
534 for( i = 0, pm = pp - 2; i < n; i++, pp++,
535 *pp = *pm;
536 pp->x *= -1;
537 }
538 gaussian->length = 2 * n + 1;
539
540 return( gaussian );
541
542 Err:
543 ptwXY_free( gaussian );
544 return( NULL );
545 }
546 /*
547 **********************************************
548 */
549 static nfu_status ptwXY_createGaussianCentered
550
551 nfu_status status = nfu_Okay;
552 int morePoints = 0;
553 double x = 0.5 * ( x1 + x2 );
554 double y = G4Exp( -x * x / 2 ), yMin = ( y
555
556 if( std::fabs( y - yMin ) > y * ptwXY->acc
557 if( morePoints && ( status = ptwXY_createG
558 if( ( status = ptwXY_setValueAtX( ptwXY, x
559 if( morePoints && ( status = ptwXY_createG
560 if( addX1Point ) status = ptwXY_setValueAt
561 return( status );
562 }
563 /*
564 **********************************************
565 */
566 ptwXYPoints *ptwXY_createGaussian( double accu
567 double /*dullEps*/, nfu_status *status
568
569 int64_t i;
570 ptwXYPoints *gaussian, *sliced;
571 ptwXYPoint *point;
572
573 if( ( gaussian = ptwXY_createGaussianCente
574 for( i = 0, point = gaussian->points; i <
575 point->x = point->x * sigma + xCenter;
576 point->y *= amplitude;
577 }
578 if( ( gaussian->points[0].x < xMin ) || (
579 if( ( sliced = ptwXY_xSlice( gaussian,
580 ptwXY_free( gaussian );
581 gaussian = sliced;
582 }
583
584 return( gaussian );
585
586 Err:
587 ptwXY_free( gaussian );
588 return( NULL );
589 }
590
591 #if defined __cplusplus
592 }
593 #endif
594