| 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 #include <nf_Legendre.h>
11 #include <nf_integration.h>
12
13 #if defined __cplusplus
14 #include "G4Log.hh"
15 #include "G4Pow.hh"
16 namespace GIDI {
17 using namespace GIDI;
18 #endif
19
20 typedef struct ptwXY_integrateWithFunctionInfo_s {
21 int degree;
22 ptwXY_createFromFunction_callback func;
23 void *argList;
24 ptwXY_interpolation interpolation;
25 double x1, x2, y1, y2;
26 } ptwXY_integrateWithFunctionInfo;
27
28 static nfu_status ptwXY_integrateWithFunction2( nf_Legendre_GaussianQuadrature_callback integrandFunction, void *argList, double x1,
29 double x2, double *integral );
30 static nfu_status ptwXY_integrateWithFunction3( double x, double *y, void *argList );
31 /*
32 ************************************************************
33 */
34 nfu_status ptwXY_f_integrate( ptwXY_interpolation interpolation, double x1, double y1, double x2, double y2, double *value ) {
35
36 nfu_status status = nfu_Okay;
37 double r;
38
39 *value = 0.;
40 switch( interpolation ) {
41 case ptwXY_interpolationLinLin : /* x linear, y linear */
42 *value = 0.5 * ( y1 + y2 ) * ( x2 - x1 );
43 break;
44 case ptwXY_interpolationLinLog : /* x linear, y log */
45 if( ( y1 <= 0. ) || ( y2 <= 0. ) ) {
46 status = nfu_badIntegrationInput; }
47 else {
48 r = y2 / y1;
49 if( std::fabs( r - 1. ) < 1e-4 ) {
50 r = r - 1.;
51 *value = y1 * ( x2 - x1 ) / ( 1. + r * ( -0.5 + r * ( 1. / 3. + r * ( -0.25 + .2 * r ) ) ) ); }
52 else {
53 *value = ( y2 - y1 ) * ( x2 - x1 ) / G4Log( r );
54 }
55 }
56 break;
57 case ptwXY_interpolationLogLin : /* x log, y linear */
58 if( ( x1 <= 0. ) || ( x2 <= 0. ) ) {
59 status = nfu_badIntegrationInput; }
60 else {
61 r = x2 / x1;
62 if( std::fabs( r - 1. ) < 1e-4 ) {
63 r = r - 1.;
64 r = r * ( -0.5 + r * ( 1. / 3. + r * ( -0.25 + .2 * r ) ) );
65 *value = x1 * ( y2 - y1 ) * r / ( 1. + r ) + y2 * ( x2 - x1 ); }
66 else {
67 *value = ( y1 - y2 ) * ( x2 - x1 ) / G4Log( r ) + x2 * y2 - x1 * y1;
68 }
69 }
70 break;
71 case ptwXY_interpolationLogLog : /* x log, y log */
72 if( ( x1 <= 0. ) || ( x2 <= 0. ) || ( y1 <= 0. ) || ( y2 <= 0. ) ) {
73 status = nfu_badIntegrationInput; }
74 else {
75 int i, n;
76 double a, z, lx, ly, s, f;
77
78 r = y2 / y1;
79 if( std::fabs( r - 1. ) < 1e-4 ) {
80 ly = ( y2 - y1 ) / y1;
81 ly = ly * ( 1. + ly * ( -0.5 + ly * ( 1. / 3. - 0.25 * ly ) ) ); }
82 else {
83 ly = G4Log( r );
84 }
85 r = x2 / x1;
86 if( std::fabs( r - 1. ) < 1e-4 ) {
87 lx = ( x2 - x1 ) / x1;
88 lx = lx * ( 1 + lx * ( -0.5 + lx * ( 1. / 3. - 0.25 * lx ) ) ); }
89 else {
90 lx = G4Log( r );
91 }
92 a = ly / lx;
93 if( std::fabs( r - 1. ) < 1e-3 ) {
94 z = ( x2 - x1 ) / x1;
95 n = (int) a;
96 if( n > 10 ) n = 12;
97 if( n < 4 ) n = 6;
98 a = a - n + 1;
99 f = n + 1.;
100 for( i = 0, s = 0.; i < n; i++, a++, f-- ) s = ( 1. + s ) * a * z / f;
101 *value = y1 * ( x2 - x1 ) * ( 1. + s ); }
102 else {
103 *value = y1 * x1 * ( G4Pow::GetInstance()->powA( r, a + 1. ) - 1. ) / ( a + 1. );
104 }
105 }
106 break;
107 case ptwXY_interpolationFlat : /* x ?, y flat */
108 *value = y1 * ( x2 - x1 );
109 break;
110 case ptwXY_interpolationOther :
111 status = nfu_otherInterpolation;
112 }
113 return( status );
114 }
115 /*
116 ************************************************************
117 */
118 double ptwXY_integrate( ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status ) {
119
120 int64_t i, n = ptwXY->length;
121 double sum = 0., dSum, x, y, x1, x2, y1, y2, _sign = 1.;
122 ptwXYPoint *point;
123
124 if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
125 *status = nfu_otherInterpolation;
126 if( ptwXY->interpolation == ptwXY_interpolationOther ) return( 0. );
127
128 if( xMax < xMin ) {
129 x = xMin;
130 xMin = xMax;
131 xMax = x;
132 _sign = -1.;
133 }
134 if( n < 2 ) return( 0. );
135
136 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( 0. );
137 for( i = 0, point = ptwXY->points; i < n; i++, point++ ) {
138 if( point->x >= xMin ) break;
139 }
140 if( i == n ) return( 0. );
141 x2 = point->x;
142 y2 = point->y;
143 if( i > 0 ) {
144 if( x2 > xMin ) {
145 x1 = point[-1].x;
146 y1 = point[-1].y;
147 if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMin, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( 0. );
148 if( x2 > xMax ) {
149 double yMax;
150
151 if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMax, &yMax, x1, y1, x2, y2 ) ) != nfu_Okay ) return( 0. );
152 if( ( *status = ptwXY_f_integrate( ptwXY->interpolation, xMin, y, xMax, yMax, &sum ) ) != nfu_Okay ) return( 0. );
153 return( sum ); }
154 else {
155 if( ( *status = ptwXY_f_integrate( ptwXY->interpolation, xMin, y, x2, y2, &sum ) ) != nfu_Okay ) return( 0. );
156 }
157 }
158 }
159 i++;
160 point++;
161 for( ; i < n; i++, point++ ) {
162 x1 = x2;
163 y1 = y2;
164 x2 = point->x;
165 y2 = point->y;
166 if( x2 > xMax ) {
167 if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMax, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( 0. );
168 if( ( *status = ptwXY_f_integrate( ptwXY->interpolation, x1, y1, xMax, y, &dSum ) ) != nfu_Okay ) return( 0. );
169 sum += dSum;
170 break;
171 }
172 if( ( *status = ptwXY_f_integrate( ptwXY->interpolation, x1, y1, x2, y2, &dSum ) ) != nfu_Okay ) return( 0. );
173 sum += dSum;
174 }
175
176 return( _sign * sum );
177 }
178 /*
179 ************************************************************
180 */
181 double ptwXY_integrateDomain( ptwXYPoints *ptwXY, nfu_status *status ) {
182
183 if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
184 if( ptwXY->length > 0 ) return( ptwXY_integrate( ptwXY, ptwXY_getXMin( ptwXY ), ptwXY_getXMax( ptwXY ), status ) );
185 return( 0. );
186 }
187 /*
188 ************************************************************
189 */
190 nfu_status ptwXY_normalize( ptwXYPoints *ptwXY ) {
191 /*
192 * This function assumes ptwXY_integrateDomain checks status and coalesces the points.
193 */
194
195 int64_t i;
196 nfu_status status;
197 double sum = ptwXY_integrateDomain( ptwXY, &status );
198
199 if( status != nfu_Okay ) return( status );
200 if( sum == 0. ) {
201 status = nfu_badNorm; }
202 else {
203 for( i = 0; i < ptwXY->length; i++ ) ptwXY->points[i].y /= sum;
204 }
205 return( status );
206 }
207 /*
208 ************************************************************
209 */
210 double ptwXY_integrateDomainWithWeight_x( ptwXYPoints *ptwXY, nfu_status *status ) {
211
212 if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
213 if( ptwXY->length < 2 ) return( 0. );
214 return( ptwXY_integrateWithWeight_x( ptwXY, ptwXY_getXMin( ptwXY ), ptwXY_getXMax( ptwXY ), status ) );
215 }
216 /*
217 ************************************************************
218 */
219 double ptwXY_integrateWithWeight_x( ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status ) {
220
221 int64_t i, n = ptwXY->length;
222 double sum = 0., x, y, x1, x2, y1, y2, _sign = 1.;
223 ptwXYPoint *point;
224
225 if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
226 *status = nfu_unsupportedInterpolation;
227 if( ( ptwXY->interpolation != ptwXY_interpolationLinLin ) &&
228 ( ptwXY->interpolation != ptwXY_interpolationFlat ) ) return( 0. );
229
230 if( n < 2 ) return( 0. );
231 if( xMax < xMin ) {
232 x = xMin;
233 xMin = xMax;
234 xMax = x;
235 _sign = -1.;
236 }
237
238 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( 0. );
239 for( i = 0, point = ptwXY->points; i < n; ++i, ++point ) {
240 if( point->x >= xMin ) break;
241 }
242 if( i == n ) return( 0. );
243 x2 = point->x;
244 y2 = point->y;
245 if( i > 0 ) {
246 if( x2 > xMin ) {
247 if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMin, &y, point[-1].x, point[-1].y, x2, y2 ) ) != nfu_Okay ) return( 0. );
248 x2 = xMin;
249 y2 = y;
250 --i;
251 --point;
252 }
253 }
254 ++i;
255 ++point;
256 for( ; i < n; ++i, ++point ) {
257 x1 = x2;
258 y1 = y2;
259 x2 = point->x;
260 y2 = point->y;
261 if( x2 > xMax ) {
262 if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMax, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( 0. );
263 x2 = xMax;
264 y2 = y;
265 }
266 switch( ptwXY->interpolation ) {
267 case ptwXY_interpolationFlat :
268 sum += ( x2 - x1 ) * y1 * 3 * ( x1 + x2 );
269 break;
270 case ptwXY_interpolationLinLin :
271 sum += ( x2 - x1 ) * ( y1 * ( 2 * x1 + x2 ) + y2 * ( x1 + 2 * x2 ) );
272 break;
273 default : /* Only to stop compilers from complaining. */
274 break;
275 }
276 if( x2 == xMax ) break;
277 }
278
279 return( _sign * sum / 6 );
280 }
281 /*
282 ************************************************************
283 */
284 double ptwXY_integrateDomainWithWeight_sqrt_x( ptwXYPoints *ptwXY, nfu_status *status ) {
285
286 if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
287 if( ptwXY->length < 2 ) return( 0. );
288 return( ptwXY_integrateWithWeight_sqrt_x( ptwXY, ptwXY_getXMin( ptwXY ), ptwXY_getXMax( ptwXY ), status ) );
289 }
290 /*
291 ************************************************************
292 */
293 double ptwXY_integrateWithWeight_sqrt_x( ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status ) {
294
295 int64_t i, n = ptwXY->length;
296 double sum = 0., x, y, x1, x2, y1, y2, _sign = 1., sqrt_x1, sqrt_x2, inv_apb, c;
297 ptwXYPoint *point;
298
299 if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
300 *status = nfu_unsupportedInterpolation;
301 if( ( ptwXY->interpolation != ptwXY_interpolationLinLin ) &&
302 ( ptwXY->interpolation != ptwXY_interpolationFlat ) ) return( 0. );
303
304 if( n < 2 ) return( 0. );
305 if( xMax < xMin ) {
306 x = xMin;
307 xMin = xMax;
308 xMax = x;
309 _sign = -1.;
310 }
311
312 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( 0. );
313 for( i = 0, point = ptwXY->points; i < n; ++i, ++point ) {
314 if( point->x >= xMin ) break;
315 }
316 if( i == n ) return( 0. );
317 x2 = point->x;
318 y2 = point->y;
319 if( i > 0 ) {
320 if( x2 > xMin ) {
321 if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMin, &y, point[-1].x, point[-1].y, x2, y2 ) ) != nfu_Okay ) return( 0. );
322 x2 = xMin;
323 y2 = y;
324 --i;
325 --point;
326 }
327 }
328 ++i;
329 ++point;
330 sqrt_x2 = std::sqrt( x2 );
331 for( ; i < n; ++i, ++point ) {
332 x1 = x2;
333 y1 = y2;
334 sqrt_x1 = sqrt_x2;
335 x2 = point->x;
336 y2 = point->y;
337 if( x2 > xMax ) {
338 if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMax, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( 0. );
339 x2 = xMax;
340 y2 = y;
341 }
342 sqrt_x2 = std::sqrt( x2 );
343 inv_apb = sqrt_x1 + sqrt_x2;
344 c = 2. * ( sqrt_x1 * sqrt_x2 + x1 + x2 );
345 switch( ptwXY->interpolation ) {
346 case ptwXY_interpolationFlat :
347 sum += ( sqrt_x2 - sqrt_x1 ) * y1 * 2.5 * c;
348 break;
349 case ptwXY_interpolationLinLin :
350 sum += ( sqrt_x2 - sqrt_x1 ) * ( y1 * ( c + x1 * ( 1. + sqrt_x2 / inv_apb ) ) + y2 * ( c + x2 * ( 1. + sqrt_x1 / inv_apb ) ) );
351 break;
352 default : /* Only to stop compilers from complaining. */
353 break;
354 }
355 if( x2 == xMax ) break;
356 }
357
358 return( 2. / 15. * _sign * sum );
359 }
360 /*
361 ************************************************************
362 */
363 ptwXPoints *ptwXY_groupOneFunction( ptwXYPoints *ptwXY, ptwXPoints *groupBoundaries, ptwXY_group_normType normType, ptwXPoints *ptwX_norm, nfu_status *status ) {
364
365 int64_t i, igs, ngs;
366 double x1, y1, x2, y2, y2p, xg1, xg2, sum;
367 ptwXYPoints *f;
368 ptwXPoints *groupedData = NULL;
369
370 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( NULL );
371 if( ( *status = groupBoundaries->status ) != nfu_Okay ) return( NULL );
372 *status = nfu_otherInterpolation;
373 if( ptwXY->interpolation == ptwXY_interpolationOther ) return( NULL );
374
375 ngs = ptwX_length( groupBoundaries ) - 1;
376 if( normType == ptwXY_group_normType_norm ) {
377 if( ptwX_norm == NULL ) {
378 *status = nfu_badNorm;
379 return( NULL );
380 }
381 *status = ptwX_norm->status;
382 if( ptwX_norm->status != nfu_Okay ) return( NULL );
383 if( ptwX_length( ptwX_norm ) != ngs ) {
384 *status = nfu_badNorm;
385 return( NULL );
386 }
387 }
388
389 if( ( f = ptwXY_intersectionWith_ptwX( ptwXY, groupBoundaries, status ) ) == NULL ) return( NULL );
390 if( f->length == 0 ) return( ptwX_createLine( ngs, ngs, 0, 0, status ) );
391
392 if( ( groupedData = ptwX_new( ngs, status ) ) == NULL ) goto err;
393 xg1 = groupBoundaries->points[0];
394 x1 = f->points[0].x;
395 y1 = f->points[0].y;
396 for( igs = 0, i = 1; igs < ngs; igs++ ) {
397 xg2 = groupBoundaries->points[igs+1];
398 sum = 0;
399 if( xg2 > x1 ) {
400 for( ; i < f->length; i++, x1 = x2, y1 = y2 ) {
401 x2 = f->points[i].x;
402 if( x2 > xg2 ) break;
403 y2p = y2 = f->points[i].y;
404 if( f->interpolation == ptwXY_interpolationFlat ) y2p = y1;
405 sum += ( y1 + y2p ) * ( x2 - x1 );
406 }
407 }
408 if( sum != 0. ) {
409 if( normType == ptwXY_group_normType_dx ) {
410 sum /= ( xg2 - xg1 ); }
411 else if( normType == ptwXY_group_normType_norm ) {
412 if( ptwX_norm->points[igs] == 0. ) {
413 *status = nfu_divByZero;
414 goto err;
415 }
416 sum /= ptwX_norm->points[igs];
417 }
418 }
419 groupedData->points[igs] = 0.5 * sum;
420 groupedData->length++;
421 xg1 = xg2;
422 }
423
424 ptwXY_free( f );
425 return( groupedData );
426
427 err:
428 ptwXY_free( f );
429 if( groupedData != NULL ) ptwX_free( groupedData );
430 return( NULL );
431 }
432 /*
433 ************************************************************
434 */
435 ptwXPoints *ptwXY_groupTwoFunctions( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, ptwXPoints *groupBoundaries, ptwXY_group_normType normType,
436 ptwXPoints *ptwX_norm, nfu_status *status ) {
437
438 int64_t i, igs, ngs;
439 double x1, fy1, gy1, x2, fy2, gy2, fy2p, gy2p, xg1, xg2, sum;
440 ptwXYPoints *f = NULL, *ff, *g = NULL, *gg = NULL;
441 ptwXPoints *groupedData = NULL;
442
443 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
444 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY2 ) ) != nfu_Okay ) return( NULL );
445 if( ( *status = groupBoundaries->status ) != nfu_Okay ) return( NULL );
446 *status = nfu_otherInterpolation;
447 if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
448 if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( NULL );
449
450 ngs = ptwX_length( groupBoundaries ) - 1;
451 if( normType == ptwXY_group_normType_norm ) {
452 if( ptwX_norm == NULL ) {
453 *status = nfu_badNorm;
454 return( NULL );
455 }
456 if( ( *status = ptwX_norm->status ) != nfu_Okay ) return( NULL );
457 if( ptwX_length( ptwX_norm ) != ngs ) {
458 *status = nfu_badNorm;
459 return( NULL );
460 }
461 }
462
463 if( ( ff = ptwXY_intersectionWith_ptwX( ptwXY1, groupBoundaries, status ) ) == NULL ) return( NULL );
464 if( ( gg = ptwXY_intersectionWith_ptwX( ptwXY2, groupBoundaries, status ) ) == NULL ) goto err;
465 if( ( ff->length == 0 ) || ( gg->length == 0 ) ) {
466 ptwXY_free( ff );
467 ptwXY_free( gg );
468 return( ptwX_createLine( ngs, ngs, 0, 0, status ) );
469 }
470
471 if( ( *status = ptwXY_tweakDomainsToMutualify( ff, gg, 4, 0 ) ) != nfu_Okay ) goto err;
472 if( ( f = ptwXY_union( ff, gg, status, ptwXY_union_fill ) ) == NULL ) goto err;
473 if( ( g = ptwXY_union( gg, f, status, ptwXY_union_fill ) ) == NULL ) goto err;
474
475 if( ( groupedData = ptwX_new( ngs, status ) ) == NULL ) goto err;
476 xg1 = groupBoundaries->points[0];
477 x1 = f->points[0].x;
478 fy1 = f->points[0].y;
479 gy1 = g->points[0].y;
480 for( igs = 0, i = 1; igs < ngs; igs++ ) {
481 xg2 = groupBoundaries->points[igs+1];
482 sum = 0;
483 if( xg2 > x1 ) {
484 for( ; i < f->length; i++, x1 = x2, fy1 = fy2, gy1 = gy2 ) {
485 x2 = f->points[i].x;
486 if( x2 > xg2 ) break;
487 fy2p = fy2 = f->points[i].y;
488 if( f->interpolation == ptwXY_interpolationFlat ) fy2p = fy1;
489 gy2p = gy2 = g->points[i].y;
490 if( g->interpolation == ptwXY_interpolationFlat ) gy2p = gy1;
491 sum += ( ( fy1 + fy2p ) * ( gy1 + gy2p ) + fy1 * gy1 + fy2p * gy2p ) * ( x2 - x1 );
492 }
493 }
494 if( sum != 0. ) {
495 if( normType == ptwXY_group_normType_dx ) {
496 sum /= ( xg2 - xg1 ); }
497 else if( normType == ptwXY_group_normType_norm ) {
498 if( ptwX_norm->points[igs] == 0. ) {
499 *status = nfu_divByZero;
500 goto err;
501 }
502 sum /= ptwX_norm->points[igs];
503 }
504 }
505 groupedData->points[igs] = sum / 6.;
506 groupedData->length++;
507 xg1 = xg2;
508 }
509
510 ptwXY_free( f );
511 ptwXY_free( g );
512 ptwXY_free( ff );
513 ptwXY_free( gg );
514 return( groupedData );
515
516 err:
517 ptwXY_free( ff );
518 if( gg != NULL ) ptwXY_free( gg );
519 // Coverity #63063
520 if( f != NULL ) ptwXY_free( f );
521 if( g != NULL ) ptwXY_free( g );
522 if( groupedData != NULL ) ptwX_free( groupedData );
523 return( NULL );
524 }
525 /*
526 ************************************************************
527 */
528 ptwXPoints *ptwXY_groupThreeFunctions( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, ptwXYPoints *ptwXY3, ptwXPoints *groupBoundaries,
529 ptwXY_group_normType normType, ptwXPoints *ptwX_norm, nfu_status *status ) {
530
531 int64_t i, igs, ngs;
532 double x1, fy1, gy1, hy1, x2, fy2, gy2, hy2, fy2p, gy2p, hy2p, xg1, xg2, sum;
533 ptwXYPoints *f = NULL, *ff, *fff = NULL, *g = NULL, *gg = NULL, *h = NULL, *hh = NULL;
534 ptwXPoints *groupedData = NULL;
535
536 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
537 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY2 ) ) != nfu_Okay ) return( NULL );
538 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY3 ) ) != nfu_Okay ) return( NULL );
539 if( ( *status = groupBoundaries->status ) != nfu_Okay ) return( NULL );
540 *status = nfu_otherInterpolation;
541 if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
542 if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( NULL );
543 if( ptwXY3->interpolation == ptwXY_interpolationOther ) return( NULL );
544
545 ngs = ptwX_length( groupBoundaries ) - 1;
546 if( normType == ptwXY_group_normType_norm ) {
547 if( ptwX_norm == NULL ) {
548 *status = nfu_badNorm;
549 return( NULL );
550 }
551 if( ( *status = ptwX_norm->status ) != nfu_Okay ) return( NULL );
552 if( ptwX_length( ptwX_norm ) != ngs ) {
553 *status = nfu_badNorm;
554 return( NULL );
555 }
556 }
557
558 if( ( ff = ptwXY_intersectionWith_ptwX( ptwXY1, groupBoundaries, status ) ) == NULL ) return( NULL );
559 if( ( gg = ptwXY_intersectionWith_ptwX( ptwXY2, groupBoundaries, status ) ) == NULL ) goto err;
560 if( ( hh = ptwXY_intersectionWith_ptwX( ptwXY3, groupBoundaries, status ) ) == NULL ) goto err;
561 if( ( ff->length == 0 ) || ( gg->length == 0 ) || ( hh->length == 0 ) ) return( ptwX_createLine( ngs, ngs, 0, 0, status ) );
562
563 if( ( *status = ptwXY_tweakDomainsToMutualify( ff, gg, 4, 0 ) ) != nfu_Okay ) goto err;
564 if( ( *status = ptwXY_tweakDomainsToMutualify( ff, hh, 4, 0 ) ) != nfu_Okay ) goto err;
565 if( ( *status = ptwXY_tweakDomainsToMutualify( gg, hh, 4, 0 ) ) != nfu_Okay ) goto err;
566 if( ( fff = ptwXY_union( ff, gg, status, ptwXY_union_fill ) ) == NULL ) goto err;
567 if( ( h = ptwXY_union( hh, fff, status, ptwXY_union_fill ) ) == NULL ) goto err;
568 if( ( f = ptwXY_union( fff, h, status, ptwXY_union_fill ) ) == NULL ) goto err;
569 if( ( g = ptwXY_union( gg, h, status, ptwXY_union_fill ) ) == NULL ) goto err;
570
571 if( ( groupedData = ptwX_new( ngs, status ) ) == NULL ) goto err;
572 xg1 = groupBoundaries->points[0];
573 x1 = f->points[0].x;
574 fy1 = f->points[0].y;
575 gy1 = g->points[0].y;
576 hy1 = h->points[0].y;
577 for( igs = 0, i = 1; igs < ngs; igs++ ) {
578 xg2 = groupBoundaries->points[igs+1];
579 sum = 0;
580 if( xg2 > x1 ) {
581 for( ; i < f->length; i++, x1 = x2, fy1 = fy2, gy1 = gy2, hy1 = hy2 ) {
582 x2 = f->points[i].x;
583 if( x2 > xg2 ) break;
584 fy2p = fy2 = f->points[i].y;
585 if( f->interpolation == ptwXY_interpolationFlat ) fy2p = fy1;
586 gy2p = gy2 = g->points[i].y;
587 if( g->interpolation == ptwXY_interpolationFlat ) gy2p = gy1;
588 hy2p = hy2 = h->points[i].y;
589 if( h->interpolation == ptwXY_interpolationFlat ) hy2p = hy1;
590 sum += ( ( fy1 + fy2p ) * ( gy1 + gy2p ) * ( hy1 + hy2p ) + 2 * fy1 * gy1 * hy1 + 2 * fy2p * gy2p * hy2p ) * ( x2 - x1 );
591 }
592 }
593 if( sum != 0. ) {
594 if( normType == ptwXY_group_normType_dx ) {
595 sum /= ( xg2 - xg1 ); }
596 else if( normType == ptwXY_group_normType_norm ) {
597 if( ptwX_norm->points[igs] == 0. ) {
598 *status = nfu_divByZero;
599 goto err;
600 }
601 sum /= ptwX_norm->points[igs];
602 }
603 }
604 groupedData->points[igs] = sum / 12.;
605 groupedData->length++;
606 xg1 = xg2;
607 }
608
609 ptwXY_free( f );
610 ptwXY_free( g );
611 ptwXY_free( h );
612 ptwXY_free( ff );
613 ptwXY_free( gg );
614 ptwXY_free( hh );
615 ptwXY_free( fff );
616 return( groupedData );
617
618 err:
619 ptwXY_free( ff );
620 if( fff != NULL ) ptwXY_free( fff );
621 if( gg != NULL ) ptwXY_free( gg );
622 if( hh != NULL ) ptwXY_free( hh );
623 if( f != NULL ) ptwXY_free( f );
624 if( g != NULL ) ptwXY_free( g );
625 if( h != NULL ) ptwXY_free( h );
626 if( groupedData != NULL ) ptwX_free( groupedData );
627 return( NULL );
628 }
629 /*
630 ************************************************************
631 */
632 ptwXPoints *ptwXY_runningIntegral( ptwXYPoints *ptwXY, nfu_status *status ) {
633
634 int i;
635 ptwXPoints *runningIntegral = NULL;
636 double integral = 0., sum;
637
638 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( NULL );
639 if( ( runningIntegral = ptwX_new( ptwXY->length, status ) ) == NULL ) goto err;
640
641 if( ( *status = ptwX_setPointAtIndex( runningIntegral, 0, 0. ) ) != nfu_Okay ) goto err;
642 for( i = 1; i < ptwXY->length; i++ ) {
643 if( ( *status = ptwXY_f_integrate( ptwXY->interpolation, ptwXY->points[i-1].x, ptwXY->points[i-1].y,
644 ptwXY->points[i].x, ptwXY->points[i].y, &sum ) ) != nfu_Okay ) goto err;
645 integral += sum;
646 if( ( *status = ptwX_setPointAtIndex( runningIntegral, i, integral ) ) != nfu_Okay ) goto err;
647 }
648 return( runningIntegral );
649
650 err:
651 if( runningIntegral != NULL ) ptwX_free( runningIntegral );
652 return( NULL );
653 }
654 /*
655 ************************************************************
656 */
657 double ptwXY_integrateWithFunction( ptwXYPoints *ptwXY, ptwXY_createFromFunction_callback func, void *argList,
658 double xMin, double xMax, int degree, int recursionLimit, double tolerance, nfu_status *status ) {
659
660 int64_t i1, i2, n1 = ptwXY->length;
661 long evaluations;
662 double integral = 0., integral_, sign = -1., xa, xb;
663 ptwXY_integrateWithFunctionInfo integrateWithFunctionInfo;
664 ptwXYPoint *point;
665
666 if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
667
668 if( xMin == xMax ) return( 0. );
669 if( n1 < 2 ) return( 0. );
670
671 ptwXY_simpleCoalescePoints( ptwXY );
672
673 if( xMin > xMax ) {
674 sign = xMin;
675 xMin = xMax;
676 xMax = sign;
677 sign = -1.;
678 }
679 if( xMin >= ptwXY->points[n1-1].x ) return( 0. );
680 if( xMax <= ptwXY->points[0].x ) return( 0. );
681
682 for( i1 = 0; i1 < ( n1 - 1 ); i1++ ) {
683 if( ptwXY->points[i1+1].x > xMin ) break;
684 }
685 for( i2 = n1 - 1; i2 > i1; i2-- ) {
686 if( ptwXY->points[i2-1].x < xMax ) break;
687 }
688 point = &(ptwXY->points[i1]);
689
690 integrateWithFunctionInfo.degree = degree;
691 integrateWithFunctionInfo.func = func;
692 integrateWithFunctionInfo.argList = argList;
693 integrateWithFunctionInfo.interpolation = ptwXY->interpolation;
694 integrateWithFunctionInfo.x2 = point->x;
695 integrateWithFunctionInfo.y2 = point->y;
696
697 xa = xMin;
698 for( ; i1 < i2; i1++ ) {
699 integrateWithFunctionInfo.x1 = integrateWithFunctionInfo.x2;
700 integrateWithFunctionInfo.y1 = integrateWithFunctionInfo.y2;
701 ++point;
702 integrateWithFunctionInfo.x2 = point->x;
703 integrateWithFunctionInfo.y2 = point->y;
704 xb = point->x;
705 if( xb > xMax ) xb = xMax;
706 *status = nf_GnG_adaptiveQuadrature( ptwXY_integrateWithFunction2, ptwXY_integrateWithFunction3, &integrateWithFunctionInfo,
707 xa, xb, recursionLimit, tolerance, &integral_, &evaluations );
708 if( *status != nfu_Okay ) return( 0. );
709 integral += integral_;
710 xa = xb;
711 }
712
713 return( integral );
714 }
715 /*
716 ************************************************************
717 */
718 static nfu_status ptwXY_integrateWithFunction2( nf_Legendre_GaussianQuadrature_callback integrandFunction, void *argList, double x1,
719 double x2, double *integral ) {
720
721 ptwXY_integrateWithFunctionInfo *integrateWithFunctionInfo = (ptwXY_integrateWithFunctionInfo *) argList;
722 nfu_status status;
723
724 status = nf_Legendre_GaussianQuadrature( integrateWithFunctionInfo->degree, x1, x2, integrandFunction, argList, integral );
725 return( status );
726 }
727 /*
728 ************************************************************
729 */
730 static nfu_status ptwXY_integrateWithFunction3( double x, double *y, void *argList ) {
731
732 double yf;
733 ptwXY_integrateWithFunctionInfo *integrateWithFunctionInfo = (ptwXY_integrateWithFunctionInfo *) argList;
734 nfu_status status;
735
736 if( ( status = ptwXY_interpolatePoint( integrateWithFunctionInfo->interpolation, x, &yf,
737 integrateWithFunctionInfo->x1, integrateWithFunctionInfo->y1,
738 integrateWithFunctionInfo->x2, integrateWithFunctionInfo->y2 ) ) == nfu_Okay ) {
739 status = integrateWithFunctionInfo->func( x, y, integrateWithFunctionInfo->argList );
740 *y *= yf;
741 }
742 return( status );
743 }
744
745 #if defined __cplusplus
746 }
747 #endif
748