| 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 namespace GIDI {
13 using namespace GIDI;
14 #endif
15
16 static double ptwXY_mod2( double v, double m,
17 static nfu_status ptwXY_mul2_s_ptwXY( ptwXYPoi
18 static nfu_status ptwXY_div_s_ptwXY( ptwXYPoin
19 int level, int isNAN1, int isNAN2 );
20 static ptwXYPoints *ptwXY_div_ptwXY_forFlats(
21 static nfu_status ptwXY_getValueAtX_ignore_XOu
22 /*
23 **********************************************
24 */
25 nfu_status ptwXY_slopeOffset( ptwXYPoints *ptw
26
27 int64_t i, nonOverflowLength = ptwXY_getNo
28 ptwXYPoint *p;
29 ptwXYOverflowPoint *o, *overflowHeader = &
30
31 if( ptwXY->status != nfu_Okay ) return( pt
32
33 for( i = 0, p = ptwXY->points; i < nonOver
34 for( o = overflowHeader->next; o != overfl
35 return( ptwXY->status );
36 }
37 /*
38 **********************************************
39 */
40 nfu_status ptwXY_add_double( ptwXYPoints *ptwX
41 nfu_status ptwXY_sub_doubleFrom( ptwXYPoints *
42 nfu_status ptwXY_sub_fromDouble( ptwXYPoints *
43 nfu_status ptwXY_mul_double( ptwXYPoints *ptwX
44 nfu_status ptwXY_div_doubleFrom( ptwXYPoints *
45
46 if( value == 0. ) {
47 ptwXY->status = nfu_divByZero; }
48 else {
49 ptwXY_slopeOffset( ptwXY, 1. / value,
50 }
51 return( ptwXY->status );
52 }
53 nfu_status ptwXY_div_fromDouble( ptwXYPoints *
54 /*
55 * This does not do any infilling and it shou
56 */
57
58 int64_t i, nonOverflowLength = ptwXY_getNo
59 ptwXYPoint *p;
60 ptwXYOverflowPoint *o, *overflowHeader = &
61
62 if( ptwXY->status != nfu_Okay ) return( pt
63 if( ptwXY->interpolation == ptwXY_interpol
64
65 for( i = 0, p = ptwXY->points; i < nonOver
66 for( o = overflowHeader->next; o != overfl
67 if( ptwXY->status != nfu_divByZero ) {
68 for( i = 0, p = ptwXY->points; i < non
69 for( o = overflowHeader->next; o != ov
70 }
71 return( ptwXY->status );
72 }
73 /*
74 **********************************************
75 */
76 nfu_status ptwXY_mod( ptwXYPoints *ptwXY, doub
77
78 int64_t i, nonOverflowLength = ptwXY_getNo
79 ptwXYPoint *p;
80 ptwXYOverflowPoint *o, *overflowHeader = &
81
82 if( ptwXY->status != nfu_Okay ) return( pt
83 if( m == 0 ) return( ptwXY->status = nfu_d
84
85 for( i = 0, p = ptwXY->points; i < nonOver
86 for( o = overflowHeader->next; o != overfl
87 return( ptwXY->status );
88 }
89 /*
90 **********************************************
91 */
92 static double ptwXY_mod2( double v, double m,
93
94 double r = std::fmod( std::fabs( v ), std:
95
96 if( pythonMod ) {
97 if( ( v * m ) < 0. ) r = std::fabs( m
98 if( m < 0. ) r *= -1.; }
99 else {
100 if( v < 0. ) r *= -1.;
101 }
102
103 return( r );
104 }
105 /*
106 **********************************************
107 */
108 ptwXYPoints *ptwXY_binary_ptwXY( ptwXYPoints *
109
110 int64_t i;
111 int unionOptions = ptwXY_union_fill | ptwX
112 double y;
113 ptwXYPoints *n;
114 ptwXYPoint *p;
115
116 *status = nfu_otherInterpolation;
117 if( ptwXY1->interpolation == ptwXY_interpo
118 if( ptwXY2->interpolation == ptwXY_interpo
119 if( ( *status = ptwXY_areDomainsMutual( pt
120 if( ( ptwXY1->interpolation == ptwXY_inter
121 *status = nfu_invalidInterpolation;
122 if( ( ptwXY1->interpolation != ptwXY2-
123 }
124 if( ( n = ptwXY_union( ptwXY1, ptwXY2, sta
125 for( i = 0, p = n->points; i < n->leng
126 if( ( *status = ptwXY_getValueAtX_
127 p->y = v1 * p->y + v2 * y + v1v2 *
128 }
129 }
130 return( n );
131 Err:
132 if( n ) ptwXY_free( n );
133 return( NULL );
134 }
135 /*
136 **********************************************
137 */
138 ptwXYPoints *ptwXY_add_ptwXY( ptwXYPoints *ptw
139
140 ptwXYPoints *sum;
141
142 if( ptwXY1->length == 0 ) {
143 sum = ptwXY_clone( ptwXY2, status ); }
144 else if( ptwXY2->length == 0 ) {
145 sum = ptwXY_clone( ptwXY1, status ); }
146 else {
147 sum = ptwXY_binary_ptwXY( ptwXY1, ptwX
148 }
149 return( sum );
150 }
151 /*
152 **********************************************
153 */
154 ptwXYPoints *ptwXY_sub_ptwXY( ptwXYPoints *ptw
155
156 ptwXYPoints *diff;
157
158 if( ptwXY1->length == 0 ) {
159 diff = ptwXY_clone( ptwXY2, status );
160 if( ( *status = ptwXY_neg( diff ) ) !=
161 else if( ptwXY2->length == 0 ) {
162 diff = ptwXY_clone( ptwXY1, status );
163 else {
164 diff = ptwXY_binary_ptwXY( ptwXY1, ptw
165 }
166 return( diff );
167 }
168 /*
169 **********************************************
170 */
171 ptwXYPoints *ptwXY_mul_ptwXY( ptwXYPoints *ptw
172
173 ptwXYPoints *mul;
174
175 if( ptwXY1->length == 0 ) {
176 mul = ptwXY_clone( ptwXY1, status ); }
177 else if( ptwXY2->length == 0 ) {
178 mul = ptwXY_clone( ptwXY2, status ); }
179 else {
180 mul = ptwXY_binary_ptwXY( ptwXY1, ptwX
181 }
182 return( mul );
183 }
184 /*
185 **********************************************
186 */
187 ptwXYPoints *ptwXY_mul2_ptwXY( ptwXYPoints *pt
188
189 int64_t i, length;
190 ptwXYPoints *n = NULL;
191 int found;
192 double x1, y1, x2, y2, u1, u2, v1, v2, xz1
193
194 *status = nfu_otherInterpolation;
195 if( ptwXY1->interpolation == ptwXY_interpo
196 if( ptwXY2->interpolation == ptwXY_interpo
197 if( ( n = ptwXY_mul_ptwXY( ptwXY1, ptwXY2,
198 if( ptwXY1->interpolation == ptwXY_interpo
199 if( ptwXY2->interpolation == ptwXY_interpo
200 length = n->length - 1;
201 if( length > 0 ) {
202 x2 = n->points[length].x;
203 for( i = length - 1; i >= 0; i-- ) {
204 x1 = n->points[i].x;
205 if( ( *status = ptwXY_getValueAtX_
206 if( ( *status = ptwXY_getValueAtX_
207 if( ( *status = ptwXY_getValueAtX_
208 if( ( *status = ptwXY_getValueAtX_
209 found = 0;
210 if( u1 * u2 < 0 ) {
211 xz1 = ( u1 * x2 - u2 * x1 ) /
212 if( ( *status = ptwXY_setValue
213 found = 1;
214 }
215 if( v1 * v2 < 0 ) {
216 xz2 = ( v1 * x2 - v2 * x1 ) /
217 if( ( *status = ptwXY_setValue
218 found += 1;
219 }
220 if( found > 1 ) {
221 x = 0.5 * ( xz1 + xz2 );
222 if( ( *status = ptwXY_getValue
223 if( ( *status = ptwXY_getValue
224 if( ( *status = ptwXY_setValue
225 }
226 x2 = x1;
227 }
228
229 if( ( *status = ptwXY_simpleCoalescePo
230 length = n->length;
231 x2 = n->points[n->length-1].x;
232 y2 = n->points[n->length-1].y;
233 for( i = n->length - 2; i >= 0; i-- )
234 x1 = n->points[i].x;
235 y1 = n->points[i].y;
236 if( ( *status = ptwXY_mul2_s_ptwXY
237 x2 = x1;
238 y2 = y1;
239 }
240 ptwXY_update_biSectionMax( n, (double)
241 }
242 return( n );
243
244 Err:
245 if( n ) ptwXY_free( n );
246 return( NULL );
247 }
248 /*
249 **********************************************
250 */
251 static nfu_status ptwXY_mul2_s_ptwXY( ptwXYPoi
252
253 nfu_status status;
254 double u1, u2, v1, v2, x, y, yp, dx, a1, a
255
256 if( ( x2 - x1 ) < ClosestAllowXFactor * DB
257 if( level >= n->biSectionMax ) return( nfu
258 level++;
259 if( ( status = ptwXY_getValueAtX_ignore_XO
260 if( ( status = ptwXY_getValueAtX_ignore_XO
261 if( ( status = ptwXY_getValueAtX_ignore_XO
262 if( ( status = ptwXY_getValueAtX_ignore_XO
263 if( ( u1 == u2 ) || ( v1 == v2 ) ) return(
264 a1 = u1 * v1;
265 if( y1 == 0 ) a1 = 0.;
266 a2 = u2 * v2;
267 if( y2 == 0 ) a2 = 0.;
268 if( ( a1 == 0. ) || ( a2 == 0. ) ) {
269 x = 0.5 * ( x1 + x2 ); }
270 else {
271 if( ( a1 * a2 < 0. ) ) return( nfu_Oka
272 a1 = std::sqrt( std::fabs( a1 ) );
273 a2 = std::sqrt( std::fabs( a2 ) );
274 x = ( a2 * x1 + a1 * x2 ) / ( a2 + a1
275 }
276 dx = x2 - x1;
277 yp = ( u1 * v1 * ( x2 - x ) + u2 * v2 * (
278 y = ( u1 * ( x2 - x ) + u2 * ( x - x1 ) )
279 if( std::fabs( y - yp ) < std::fabs( y * n
280 if( ( status = ptwXY_setValueAtX( n, x, y
281 if( ( status = ptwXY_mul2_s_ptwXY( n, ptwX
282 status = ptwXY_mul2_s_ptwXY( n, ptwXY1, pt
283 return( status );
284 }
285 /*
286 **********************************************
287 */
288 ptwXYPoints *ptwXY_div_ptwXY( ptwXYPoints *ptw
289
290 int isNAN1, isNAN2;
291 int64_t i, j, k, zeros = 0, length, iYs;
292 double x1, x2, y1, y2, u1, u2, v1, v2, y,
293 ptwXYPoints *n = NULL;
294 ptwXYPoint *p;
295
296 if( ( *status = ptwXY_simpleCoalescePoints
297 if( ( *status = ptwXY_simpleCoalescePoints
298 *status = nfu_otherInterpolation;
299 if( ptwXY1->interpolation == ptwXY_interpo
300 if( ptwXY2->interpolation == ptwXY_interpo
301 if( ( ptwXY1->interpolation == ptwXY_inter
302 return( ptwXY_div_ptwXY_forFlats( ptwX
303
304 if( ( *status = ptwXY_areDomainsMutual( pt
305 if( ( n = ptwXY_union( ptwXY1, ptwXY2, sta
306 for( i = 0, p = n->points; i < n->leng
307 if( ( *status = ptwXY_getValueAtX_
308 if( y == 0. ) {
309 if( p->y == 0. ) {
310 iYs = 0;
311 y1 = 0.;
312 y2 = 0.;
313 if( i > 0 ) {
314 if( ( *status = ptwXY_
315 if( *status != nfu
316 s1 = 0.;
317 }
318 if( ( *status = ptwXY_
319 if( s2 == 0. ) {
320 y1 = nan; }
321 else {
322 y1 = s1 / s2;
323 }
324 iYs++;
325 }
326 if( i < ( n->length - 1 )
327 if( ( *status = ptwXY_
328 if( *status != nfu
329 s1 = 0.;
330 }
331 if( ( *status = ptwXY_
332 if( s2 == 0. ) {
333 y2 = nan; }
334 else {
335 y2 = s1 / s2;
336 }
337 iYs++;
338 }
339 p->y = ( y1 + y2 ) / iYs;
340 if( nfu_isNAN( p->y ) ) ze
341 else {
342 if( !safeDivide ) {
343 *status = nfu_divByZer
344 goto Err;
345 }
346 zeros++;
347 p->y = nan;
348 } }
349 else {
350 p->y /= y;
351 }
352 }
353 length = n->length - 1;
354 if( length > 0 ) {
355 x2 = n->points[length].x;
356 for( i = length - 1; i >= 0; i-- )
357 x1 = n->points[i].x;
358 if( ( *status = ptwXY_getValue
359 if( ( *status = ptwXY_getValue
360 if( ( *status = ptwXY_getValue
361 if( ( *status = ptwXY_getValue
362 if( u1 * u2 < 0 ) {
363 xz = ( u1 * x2 - u2 * x1 )
364 if( ( *status = ptwXY_setV
365 }
366 if( v1 * v2 < 0 ) {
367 if( !safeDivide ) {
368 *status = nfu_divByZer
369 goto Err;
370 }
371 zeros++;
372 xz = ( v1 * x2 - v2 * x1 )
373 if( ( *status = ptwXY_setV
374 }
375 x2 = x1;
376 }
377 if( ( *status = ptwXY_simpleCoales
378 length = n->length;
379 x2 = n->points[n->length-1].x;
380 y2 = n->points[n->length-1].y;
381 isNAN2 = nfu_isNAN( y2 );
382 for( i = n->length - 2; i >= 0; i-
383 x1 = n->points[i].x;
384 y1 = n->points[i].y;
385 isNAN1 = nfu_isNAN( y1 );
386 if( !isNAN1 || !isNAN2 ) {
387 if( ( *status = ptwXY_div_
388 }
389 x2 = x1;
390 y2 = y1;
391 isNAN2 = isNAN1;
392 }
393 ptwXY_update_biSectionMax( n, (dou
394 if( zeros ) {
395 if( ( *status = ptwXY_simpleCo
396 for( i = 0; i < n->length; i++
397 if( nfu_isNAN( n->points[0].y
398 if( i == n->length ) {
399 zeros = 0;
400 for( i = 0; i < n->len
401 else {
402 n->points[0].y = 2. *
403 zeros--;
404 }
405 }
406 for( i = n->length - 1; i > 0;
407 if( nfu_isNAN( n->points[n->le
408 n->points[n->length - 1].y
409 zeros--;
410 }
411 if( zeros ) {
412 for( i = 0; i < n->length;
413 for( k = i + 1, j = i; k <
414 if( nfu_isNAN( n->points[k
415 n->points[j] = n->poin
416 j++;
417 }
418 n->length = j;
419 }
420 }
421 }
422 }
423 return( n );
424
425 Err:
426 if( n ) ptwXY_free( n );
427 return( NULL );
428 }
429 /*
430 **********************************************
431 */
432 static nfu_status ptwXY_div_s_ptwXY( ptwXYPoin
433 int level, int isNAN1, int isNAN2 ) {
434
435 nfu_status status;
436 double u1, u2, v1, v2, v, x, y, yp, dx, a1
437
438 if( ( x2 - x1 ) < ClosestAllowXFactor * DB
439 if( level >= n->biSectionMax ) return( nfu
440 level++;
441 if( ( status = ptwXY_getValueAtX_ignore_XO
442 if( ( status = ptwXY_getValueAtX_ignore_XO
443 if( ( status = ptwXY_getValueAtX( ptwXY2,
444 if( ( status = ptwXY_getValueAtX( ptwXY2,
445 if( isNAN1 ) {
446 x = 0.5 * ( x1 + x2 );
447 if( ( status = ptwXY_getValueAtX_ignor
448 if( ( status = ptwXY_getValueAtX( ptwX
449 y = u1 / v1; }
450 else if( isNAN2 ) {
451 x = 0.5 * ( x1 + x2 );
452 if( ( status = ptwXY_getValueAtX_ignor
453 if( ( status = ptwXY_getValueAtX( ptwX
454 y = u2 / v2; }
455 else {
456 if( ( u1 == u2 ) || ( v1 == v2 ) ) ret
457 if( ( y1 == 0. ) || ( y2 == 0. ) ) {
458 x = 0.5 * ( x1 + x2 ); }
459 else {
460 if( ( u1 * u2 < 0. ) ) return( nfu
461 a1 = std::sqrt( std::fabs( u1 ) );
462 a2 = std::sqrt( std::fabs( u2 ) );
463 x = ( a2 * x1 + a1 * x2 ) / ( a2 +
464 }
465 dx = x2 - x1;
466 v = v1 * ( x2 - x ) + v2 * ( x - x1 );
467 if( ( v1 == 0. ) || ( v2 == 0. ) || (
468 yp = ( u1 / v1 * ( x2 - x ) + u2 / v2
469 y = ( u1 * ( x2 - x ) + u2 * ( x - x1
470 if( std::fabs( y - yp ) < std::fabs( y
471 }
472 if( ( status = ptwXY_setValueAtX( n, x, y
473 if( ( status = ptwXY_div_s_ptwXY( n, ptwXY
474 status = ptwXY_div_s_ptwXY( n, ptwXY1, ptw
475 return( status );
476 }
477 /*
478 **********************************************
479 */
480 static ptwXYPoints *ptwXY_div_ptwXY_forFlats(
481
482 int64_t i;
483 ptwXYPoints *n = NULL;
484 ptwXYPoint *p;
485 double y;
486
487 *status = nfu_invalidInterpolation;
488 if( ptwXY1->interpolation != ptwXY_interpo
489 if( ptwXY2->interpolation != ptwXY_interpo
490 if( ( n = ptwXY_union( ptwXY1, ptwXY2, sta
491 for( i = 0, p = n->points; i < n->leng
492 if( ( *status = ptwXY_getValueAtX_
493 if( y == 0. ) {
494 if( ( safeDivide ) && ( p->y =
495 *status = nfu_divByZero;
496 goto Err;
497 } }
498 else {
499 p->y /= y;
500 }
501 }
502 }
503 return( n );
504
505 Err:
506 if( n ) ptwXY_free( n );
507 return( NULL );
508 }
509 /*
510 **********************************************
511 */
512 static nfu_status ptwXY_getValueAtX_ignore_XOu
513
514 nfu_status status = ptwXY_getValueAtX( ptw
515
516 if( status == nfu_XOutsideDomain ) status
517 return( status );
518 }
519
520 #if defined __cplusplus
521 }
522 #endif
523