Geant4  10.02.p02
ptwXY_integration.cc
Go to the documentation of this file.
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  if( f != NULL ) ptwXY_free( ff );
520  if( g != NULL ) ptwXY_free( g );
521  if( groupedData != NULL ) ptwX_free( groupedData );
522  return( NULL );
523 }
524 /*
525 ************************************************************
526 */
527 ptwXPoints *ptwXY_groupThreeFunctions( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, ptwXYPoints *ptwXY3, ptwXPoints *groupBoundaries,
528  ptwXY_group_normType normType, ptwXPoints *ptwX_norm, nfu_status *status ) {
529 
530  int64_t i, igs, ngs;
531  double x1, fy1, gy1, hy1, x2, fy2, gy2, hy2, fy2p, gy2p, hy2p, xg1, xg2, sum;
532  ptwXYPoints *f = NULL, *ff, *fff = NULL, *g = NULL, *gg = NULL, *h = NULL, *hh = NULL;
533  ptwXPoints *groupedData = NULL;
534 
535  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
536  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY2 ) ) != nfu_Okay ) return( NULL );
537  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY3 ) ) != nfu_Okay ) return( NULL );
538  if( ( *status = groupBoundaries->status ) != nfu_Okay ) return( NULL );
539  *status = nfu_otherInterpolation;
540  if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
541  if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( NULL );
542  if( ptwXY3->interpolation == ptwXY_interpolationOther ) return( NULL );
543 
544  ngs = ptwX_length( groupBoundaries ) - 1;
545  if( normType == ptwXY_group_normType_norm ) {
546  if( ptwX_norm == NULL ) {
547  *status = nfu_badNorm;
548  return( NULL );
549  }
550  if( ( *status = ptwX_norm->status ) != nfu_Okay ) return( NULL );
551  if( ptwX_length( ptwX_norm ) != ngs ) {
552  *status = nfu_badNorm;
553  return( NULL );
554  }
555  }
556 
557  if( ( ff = ptwXY_intersectionWith_ptwX( ptwXY1, groupBoundaries, status ) ) == NULL ) return( NULL );
558  if( ( gg = ptwXY_intersectionWith_ptwX( ptwXY2, groupBoundaries, status ) ) == NULL ) goto err;
559  if( ( hh = ptwXY_intersectionWith_ptwX( ptwXY3, groupBoundaries, status ) ) == NULL ) goto err;
560  if( ( ff->length == 0 ) || ( gg->length == 0 ) || ( hh->length == 0 ) ) return( ptwX_createLine( ngs, ngs, 0, 0, status ) );
561 
562  if( ( *status = ptwXY_tweakDomainsToMutualify( ff, gg, 4, 0 ) ) != nfu_Okay ) goto err;
563  if( ( *status = ptwXY_tweakDomainsToMutualify( ff, hh, 4, 0 ) ) != nfu_Okay ) goto err;
564  if( ( *status = ptwXY_tweakDomainsToMutualify( gg, hh, 4, 0 ) ) != nfu_Okay ) goto err;
565  if( ( fff = ptwXY_union( ff, gg, status, ptwXY_union_fill ) ) == NULL ) goto err;
566  if( ( h = ptwXY_union( hh, fff, status, ptwXY_union_fill ) ) == NULL ) goto err;
567  if( ( f = ptwXY_union( fff, h, status, ptwXY_union_fill ) ) == NULL ) goto err;
568  if( ( g = ptwXY_union( gg, h, status, ptwXY_union_fill ) ) == NULL ) goto err;
569 
570  if( ( groupedData = ptwX_new( ngs, status ) ) == NULL ) goto err;
571  xg1 = groupBoundaries->points[0];
572  x1 = f->points[0].x;
573  fy1 = f->points[0].y;
574  gy1 = g->points[0].y;
575  hy1 = h->points[0].y;
576  for( igs = 0, i = 1; igs < ngs; igs++ ) {
577  xg2 = groupBoundaries->points[igs+1];
578  sum = 0;
579  if( xg2 > x1 ) {
580  for( ; i < f->length; i++, x1 = x2, fy1 = fy2, gy1 = gy2, hy1 = hy2 ) {
581  x2 = f->points[i].x;
582  if( x2 > xg2 ) break;
583  fy2p = fy2 = f->points[i].y;
584  if( f->interpolation == ptwXY_interpolationFlat ) fy2p = fy1;
585  gy2p = gy2 = g->points[i].y;
586  if( g->interpolation == ptwXY_interpolationFlat ) gy2p = gy1;
587  hy2p = hy2 = h->points[i].y;
588  if( h->interpolation == ptwXY_interpolationFlat ) hy2p = hy1;
589  sum += ( ( fy1 + fy2p ) * ( gy1 + gy2p ) * ( hy1 + hy2p ) + 2 * fy1 * gy1 * hy1 + 2 * fy2p * gy2p * hy2p ) * ( x2 - x1 );
590  }
591  }
592  if( sum != 0. ) {
593  if( normType == ptwXY_group_normType_dx ) {
594  sum /= ( xg2 - xg1 ); }
595  else if( normType == ptwXY_group_normType_norm ) {
596  if( ptwX_norm->points[igs] == 0. ) {
597  *status = nfu_divByZero;
598  goto err;
599  }
600  sum /= ptwX_norm->points[igs];
601  }
602  }
603  groupedData->points[igs] = sum / 12.;
604  groupedData->length++;
605  xg1 = xg2;
606  }
607 
608  ptwXY_free( f );
609  ptwXY_free( g );
610  ptwXY_free( h );
611  ptwXY_free( ff );
612  ptwXY_free( gg );
613  ptwXY_free( hh );
614  ptwXY_free( fff );
615  return( groupedData );
616 
617 err:
618  ptwXY_free( ff );
619  if( fff != NULL ) ptwXY_free( fff );
620  if( gg != NULL ) ptwXY_free( gg );
621  if( hh != NULL ) ptwXY_free( hh );
622  if( f != NULL ) ptwXY_free( f );
623  if( g != NULL ) ptwXY_free( g );
624  if( h != NULL ) ptwXY_free( h );
625  if( groupedData != NULL ) ptwX_free( groupedData );
626  return( NULL );
627 }
628 /*
629 ************************************************************
630 */
631 ptwXPoints *ptwXY_runningIntegral( ptwXYPoints *ptwXY, nfu_status *status ) {
632 
633  int i;
634  ptwXPoints *runningIntegral = NULL;
635  double integral = 0., sum;
636 
637  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( NULL );
638  if( ( runningIntegral = ptwX_new( ptwXY->length, status ) ) == NULL ) goto err;
639 
640  if( ( *status = ptwX_setPointAtIndex( runningIntegral, 0, 0. ) ) != nfu_Okay ) goto err;
641  for( i = 1; i < ptwXY->length; i++ ) {
642  if( ( *status = ptwXY_f_integrate( ptwXY->interpolation, ptwXY->points[i-1].x, ptwXY->points[i-1].y,
643  ptwXY->points[i].x, ptwXY->points[i].y, &sum ) ) != nfu_Okay ) goto err;
644  integral += sum;
645  if( ( *status = ptwX_setPointAtIndex( runningIntegral, i, integral ) ) != nfu_Okay ) goto err;
646  }
647  return( runningIntegral );
648 
649 err:
650  if( runningIntegral != NULL ) ptwX_free( runningIntegral );
651  return( NULL );
652 }
653 /*
654 ************************************************************
655 */
656 double ptwXY_integrateWithFunction( ptwXYPoints *ptwXY, ptwXY_createFromFunction_callback func, void *argList,
657  double xMin, double xMax, int degree, int recursionLimit, double tolerance, nfu_status *status ) {
658 
659  int64_t i1, i2, n1 = ptwXY->length;
660  long evaluations;
661  double integral = 0., integral_, sign = -1., xa, xb;
662  ptwXY_integrateWithFunctionInfo integrateWithFunctionInfo;
663  ptwXYPoint *point;
664 
665  if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
666 
667  if( xMin == xMax ) return( 0. );
668  if( n1 < 2 ) return( 0. );
669 
670  ptwXY_simpleCoalescePoints( ptwXY );
671 
672  if( xMin > xMax ) {
673  sign = xMin;
674  xMin = xMax;
675  xMax = sign;
676  sign = -1.;
677  }
678  if( xMin >= ptwXY->points[n1-1].x ) return( 0. );
679  if( xMax <= ptwXY->points[0].x ) return( 0. );
680 
681  for( i1 = 0; i1 < ( n1 - 1 ); i1++ ) {
682  if( ptwXY->points[i1+1].x > xMin ) break;
683  }
684  for( i2 = n1 - 1; i2 > i1; i2-- ) {
685  if( ptwXY->points[i2-1].x < xMax ) break;
686  }
687  point = &(ptwXY->points[i1]);
688 
689  integrateWithFunctionInfo.degree = degree;
690  integrateWithFunctionInfo.func = func;
691  integrateWithFunctionInfo.argList = argList;
692  integrateWithFunctionInfo.interpolation = ptwXY->interpolation;
693  integrateWithFunctionInfo.x2 = point->x;
694  integrateWithFunctionInfo.y2 = point->y;
695 
696  xa = xMin;
697  for( ; i1 < i2; i1++ ) {
698  integrateWithFunctionInfo.x1 = integrateWithFunctionInfo.x2;
699  integrateWithFunctionInfo.y1 = integrateWithFunctionInfo.y2;
700  ++point;
701  integrateWithFunctionInfo.x2 = point->x;
702  integrateWithFunctionInfo.y2 = point->y;
703  xb = point->x;
704  if( xb > xMax ) xb = xMax;
705  *status = nf_GnG_adaptiveQuadrature( ptwXY_integrateWithFunction2, ptwXY_integrateWithFunction3, &integrateWithFunctionInfo,
706  xa, xb, recursionLimit, tolerance, &integral_, &evaluations );
707  if( *status != nfu_Okay ) return( 0. );
708  integral += integral_;
709  xa = xb;
710  }
711 
712  return( integral );
713 }
714 /*
715 ************************************************************
716 */
717 static nfu_status ptwXY_integrateWithFunction2( nf_Legendre_GaussianQuadrature_callback integrandFunction, void *argList, double x1,
718  double x2, double *integral ) {
719 
720  ptwXY_integrateWithFunctionInfo *integrateWithFunctionInfo = (ptwXY_integrateWithFunctionInfo *) argList;
721  nfu_status status;
722 
723  status = nf_Legendre_GaussianQuadrature( integrateWithFunctionInfo->degree, x1, x2, integrandFunction, argList, integral );
724  return( status );
725 }
726 /*
727 ************************************************************
728 */
729 static nfu_status ptwXY_integrateWithFunction3( double x, double *y, void *argList ) {
730 
731  double yf;
732  ptwXY_integrateWithFunctionInfo *integrateWithFunctionInfo = (ptwXY_integrateWithFunctionInfo *) argList;
733  nfu_status status;
734 
735  if( ( status = ptwXY_interpolatePoint( integrateWithFunctionInfo->interpolation, x, &yf,
736  integrateWithFunctionInfo->x1, integrateWithFunctionInfo->y1,
737  integrateWithFunctionInfo->x2, integrateWithFunctionInfo->y2 ) ) == nfu_Okay ) {
738  status = integrateWithFunctionInfo->func( x, y, integrateWithFunctionInfo->argList );
739  *y *= yf;
740  }
741  return( status );
742 }
743 
744 #if defined __cplusplus
745 }
746 #endif
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4Pow::powA
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
ptwXY_integrateDomainWithWeight_sqrt_x
double ptwXY_integrateDomainWithWeight_sqrt_x(ptwXYPoints *ptwXY, nfu_status *status)
Definition: ptwXY_integration.cc:284
ptwXY_intersectionWith_ptwX
ptwXYPoints * ptwXY_intersectionWith_ptwX(ptwXYPoints *ptwXY, ptwXPoints *ptwX, nfu_status *status)
Definition: ptwXY_convenient.cc:194
ptwX_setPointAtIndex
nfu_status ptwX_setPointAtIndex(ptwXPoints *ptwX, int64_t index, double x)
Definition: ptwX_core.cc:222
ptwXY_integrateWithFunctionInfo_s::x2
double x2
Definition: ptwXY_integration.cc:25
ptwXY_normalize
nfu_status ptwXY_normalize(ptwXYPoints *ptwXY)
Definition: ptwXY_integration.cc:190
ptwXY_integrateWithFunction
double ptwXY_integrateWithFunction(ptwXYPoints *ptwXY, ptwXY_createFromFunction_callback func, void *argList, double xMin, double xMax, int degree, int recursionLimit, double tolerance, nfu_status *status)
Definition: ptwXY_integration.cc:656
ptwXY_groupThreeFunctions
ptwXPoints * ptwXY_groupThreeFunctions(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, ptwXYPoints *ptwXY3, ptwXPoints *groupBoundaries, ptwXY_group_normType normType, ptwXPoints *ptwX_norm, nfu_status *status)
Definition: ptwXY_integration.cc:527
z
G4double z
Definition: TRTMaterials.hh:39
tolerance
static const G4float tolerance
Definition: G4LevelManager.hh:54
ptwXY_simpleCoalescePoints
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
ptwXY_integrateWithFunctionInfo_s::y1
double y1
Definition: ptwXY_integration.cc:25
ptwXY_free
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
a
G4double a
Definition: TRTMaterials.hh:39
ptwXY_getXMax
double ptwXY_getXMax(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:1239
ptwXY_integrateWithFunctionInfo
struct ptwXY_integrateWithFunctionInfo_s ptwXY_integrateWithFunctionInfo
s
static const double s
Definition: G4SIunits.hh:168
ptwX_new
ptwXPoints * ptwX_new(int64_t size, nfu_status *status)
Definition: ptwX_core.cc:24
ptwXY_f_integrate
nfu_status ptwXY_f_integrate(ptwXY_interpolation interpolation, double x1, double y1, double x2, double y2, double *value)
Definition: ptwXY_integration.cc:34
ptwX_createLine
ptwXPoints * ptwX_createLine(int64_t size, int64_t length, double slope, double offset, nfu_status *status)
Definition: ptwX_core.cc:62
ptwXY_integrate
double ptwXY_integrate(ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status)
Definition: ptwXY_integration.cc:118
nf_GnG_adaptiveQuadrature
nfu_status nf_GnG_adaptiveQuadrature(nf_GnG_adaptiveQuadrature_callback quadratureFunction, nf_Legendre_GaussianQuadrature_callback integrandFunction, void *argList, double x1, double x2, int maxDepth, double tolerance, double *integral, long *evaluations)
Definition: nf_GnG_adaptiveQuadrature.cc:31
ptwXY_integrateDomainWithWeight_x
double ptwXY_integrateDomainWithWeight_x(ptwXYPoints *ptwXY, nfu_status *status)
Definition: ptwXY_integration.cc:210
ptwXY_integrateWithFunction3
static nfu_status ptwXY_integrateWithFunction3(double x, double *y, void *argList)
Definition: ptwXY_integration.cc:729
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
ptwXY_integrateWithFunction2
static nfu_status ptwXY_integrateWithFunction2(nf_Legendre_GaussianQuadrature_callback integrandFunction, void *argList, double x1, double x2, double *integral)
Definition: ptwXY_integration.cc:717
G4Log.hh
ptwX_free
ptwXPoints * ptwX_free(ptwXPoints *ptwX)
Definition: ptwX_core.cc:158
ptwXY_integrateWithFunctionInfo_s::interpolation
ptwXY_interpolation interpolation
Definition: ptwXY_integration.cc:24
ptwXY_groupOneFunction
ptwXPoints * ptwXY_groupOneFunction(ptwXYPoints *ptwXY, ptwXPoints *groupBoundaries, ptwXY_group_normType normType, ptwXPoints *ptwX_norm, nfu_status *status)
Definition: ptwXY_integration.cc:363
ptwXY_groupTwoFunctions
ptwXPoints * ptwXY_groupTwoFunctions(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, ptwXPoints *groupBoundaries, ptwXY_group_normType normType, ptwXPoints *ptwX_norm, nfu_status *status)
Definition: ptwXY_integration.cc:435
ptwXY_integrateWithFunctionInfo_s::degree
int degree
Definition: ptwXY_integration.cc:21
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:230
ptwXY_integrateWithFunctionInfo_s::y2
double y2
Definition: ptwXY_integration.cc:25
ptwX_length
int64_t ptwX_length(ptwXPoints *ptwX)
Definition: ptwX_core.cc:166
ptwXY_integrateWithWeight_sqrt_x
double ptwXY_integrateWithWeight_sqrt_x(ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status)
Definition: ptwXY_integration.cc:293
g
static const double g
Definition: G4SIunits.hh:180
x
const G4double x[NPOINTSGL]
Definition: G4PreCompoundFragment.cc:56
ptwXY_integrateWithWeight_x
double ptwXY_integrateWithWeight_x(ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status)
Definition: ptwXY_integration.cc:219
degree
static const double degree
Definition: G4SIunits.hh:143
ptwXY_union
ptwXYPoints * ptwXY_union(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int unionOptions)
Definition: ptwXY_methods.cc:349
ptwXY_runningIntegral
ptwXPoints * ptwXY_runningIntegral(ptwXYPoints *ptwXY, nfu_status *status)
Definition: ptwXY_integration.cc:631
nf_Legendre_GaussianQuadrature
nfu_status nf_Legendre_GaussianQuadrature(int degree, double x1, double x2, nf_Legendre_GaussianQuadrature_callback func, void *argList, double *integral)
Definition: nf_Legendre_GaussianQuadrature.cc:63
ptwXY_integrateWithFunctionInfo_s::x1
double x1
Definition: ptwXY_integration.cc:25
ptwXY_getXMin
double ptwXY_getXMin(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:1206
ptwXY_integrateWithFunctionInfo_s::func
ptwXY_createFromFunction_callback func
Definition: ptwXY_integration.cc:22
ptwXY_tweakDomainsToMutualify
nfu_status ptwXY_tweakDomainsToMutualify(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, int epsilonFactor, double epsilon)
Definition: ptwXY_convenient.cc:295
G4Pow.hh
ptwXY_interpolatePoint
nfu_status ptwXY_interpolatePoint(ptwXY_interpolation interpolation, double x, double *y, double x1, double y1, double x2, double y2)
Definition: ptwXY_interpolation.cc:30
ptwXY_integrateDomain
double ptwXY_integrateDomain(ptwXYPoints *ptwXY, nfu_status *status)
Definition: ptwXY_integration.cc:181
G4INCL::Math::sign
G4int sign(const T t)
A simple sign function that allows us to port fortran code to c++ more easily.
Definition: G4INCLGlobals.hh:107
ptwXY_integrateWithFunctionInfo_s
Definition: ptwXY_integration.cc:20
ptwXY_integrateWithFunctionInfo_s::argList
void * argList
Definition: ptwXY_integration.cc:23