ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ptwXY_integration.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ptwXY_integration.cc
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