| 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 #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
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
29 double x2, double *integral );
30 static nfu_status ptwXY_integrateWithFunction3
31 /*
32 **********************************************
33 */
34 nfu_status ptwXY_f_integrate( ptwXY_interpolat
35
36 nfu_status status = nfu_Okay;
37 double r;
38
39 *value = 0.;
40 switch( interpolation ) {
41 case ptwXY_interpolationLinLin :
42 *value = 0.5 * ( y1 + y2 ) * ( x2 - x1
43 break;
44 case ptwXY_interpolationLinLog :
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 ) / (
52 else {
53 *value = ( y2 - y1 ) * ( x2 -
54 }
55 }
56 break;
57 case ptwXY_interpolationLogLin :
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.
65 *value = x1 * ( y2 - y1 ) * r
66 else {
67 *value = ( y1 - y2 ) * ( x2 -
68 }
69 }
70 break;
71 case ptwXY_interpolationLogLog :
72 if( ( x1 <= 0. ) || ( x2 <= 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 +
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 +
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++
101 *value = y1 * ( x2 - x1 ) * (
102 else {
103 *value = y1 * x1 * ( G4Pow::Ge
104 }
105 }
106 break;
107 case ptwXY_interpolationFlat :
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, do
119
120 int64_t i, n = ptwXY->length;
121 double sum = 0., dSum, x, y, x1, x2, y1, y
122 ptwXYPoint *point;
123
124 if( ( *status = ptwXY->status ) != nfu_Oka
125 *status = nfu_otherInterpolation;
126 if( ptwXY->interpolation == ptwXY_interpol
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
137 for( i = 0, point = ptwXY->points; i < n;
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_interpolateP
148 if( x2 > xMax ) {
149 double yMax;
150
151 if( ( *status = ptwXY_interpol
152 if( ( *status = ptwXY_f_integr
153 return( sum ); }
154 else {
155 if( ( *status = ptwXY_f_integr
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_interpolateP
168 if( ( *status = ptwXY_f_integrate(
169 sum += dSum;
170 break;
171 }
172 if( ( *status = ptwXY_f_integrate( ptw
173 sum += dSum;
174 }
175
176 return( _sign * sum );
177 }
178 /*
179 **********************************************
180 */
181 double ptwXY_integrateDomain( ptwXYPoints *ptw
182
183 if( ( *status = ptwXY->status ) != nfu_Oka
184 if( ptwXY->length > 0 ) return( ptwXY_inte
185 return( 0. );
186 }
187 /*
188 **********************************************
189 */
190 nfu_status ptwXY_normalize( ptwXYPoints *ptwXY
191 /*
192 * This function assumes ptwXY_integrateDomai
193 */
194
195 int64_t i;
196 nfu_status status;
197 double sum = ptwXY_integrateDomain( ptwXY,
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++ ) p
204 }
205 return( status );
206 }
207 /*
208 **********************************************
209 */
210 double ptwXY_integrateDomainWithWeight_x( ptwX
211
212 if( ( *status = ptwXY->status ) != nfu_Oka
213 if( ptwXY->length < 2 ) return( 0. );
214 return( ptwXY_integrateWithWeight_x( ptwXY
215 }
216 /*
217 **********************************************
218 */
219 double ptwXY_integrateWithWeight_x( ptwXYPoint
220
221 int64_t i, n = ptwXY->length;
222 double sum = 0., x, y, x1, x2, y1, y2, _si
223 ptwXYPoint *point;
224
225 if( ( *status = ptwXY->status ) != nfu_Oka
226 *status = nfu_unsupportedInterpolation;
227 if( ( ptwXY->interpolation != ptwXY_interp
228 ( ptwXY->interpolation != ptwXY_interp
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
239 for( i = 0, point = ptwXY->points; i < n;
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_interpolateP
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_interpolateP
263 x2 = xMax;
264 y2 = y;
265 }
266 switch( ptwXY->interpolation ) {
267 case ptwXY_interpolationFlat :
268 sum += ( x2 - x1 ) * y1 * 3 * ( x1
269 break;
270 case ptwXY_interpolationLinLin :
271 sum += ( x2 - x1 ) * ( y1 * ( 2 *
272 break;
273 default : /* Only to stop compil
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(
285
286 if( ( *status = ptwXY->status ) != nfu_Oka
287 if( ptwXY->length < 2 ) return( 0. );
288 return( ptwXY_integrateWithWeight_sqrt_x(
289 }
290 /*
291 **********************************************
292 */
293 double ptwXY_integrateWithWeight_sqrt_x( ptwXY
294
295 int64_t i, n = ptwXY->length;
296 double sum = 0., x, y, x1, x2, y1, y2, _si
297 ptwXYPoint *point;
298
299 if( ( *status = ptwXY->status ) != nfu_Oka
300 *status = nfu_unsupportedInterpolation;
301 if( ( ptwXY->interpolation != ptwXY_interp
302 ( ptwXY->interpolation != ptwXY_interp
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
313 for( i = 0, point = ptwXY->points; i < n;
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_interpolateP
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_interpolateP
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
348 break;
349 case ptwXY_interpolationLinLin :
350 sum += ( sqrt_x2 - sqrt_x1 ) * ( y
351 break;
352 default : /* Only to stop compil
353 break;
354 }
355 if( x2 == xMax ) break;
356 }
357
358 return( 2. / 15. * _sign * sum );
359 }
360 /*
361 **********************************************
362 */
363 ptwXPoints *ptwXY_groupOneFunction( ptwXYPoint
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
371 if( ( *status = groupBoundaries->status )
372 *status = nfu_otherInterpolation;
373 if( ptwXY->interpolation == ptwXY_interpol
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 ) re
383 if( ptwX_length( ptwX_norm ) != ngs )
384 *status = nfu_badNorm;
385 return( NULL );
386 }
387 }
388
389 if( ( f = ptwXY_intersectionWith_ptwX( ptw
390 if( f->length == 0 ) return( ptwX_createLi
391
392 if( ( groupedData = ptwX_new( ngs, status
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
401 x2 = f->points[i].x;
402 if( x2 > xg2 ) break;
403 y2p = y2 = f->points[i].y;
404 if( f->interpolation == ptwXY_
405 sum += ( y1 + y2p ) * ( x2 - x
406 }
407 }
408 if( sum != 0. ) {
409 if( normType == ptwXY_group_normTy
410 sum /= ( xg2 - xg1 ); }
411 else if( normType == ptwXY_group_n
412 if( ptwX_norm->points[igs] ==
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( group
430 return( NULL );
431 }
432 /*
433 **********************************************
434 */
435 ptwXPoints *ptwXY_groupTwoFunctions( ptwXYPoin
436 ptwXPoints *ptwX_norm, nfu_status *sta
437
438 int64_t i, igs, ngs;
439 double x1, fy1, gy1, x2, fy2, gy2, fy2p, g
440 ptwXYPoints *f = NULL, *ff, *g = NULL, *gg
441 ptwXPoints *groupedData = NULL;
442
443 if( ( *status = ptwXY_simpleCoalescePoints
444 if( ( *status = ptwXY_simpleCoalescePoints
445 if( ( *status = groupBoundaries->status )
446 *status = nfu_otherInterpolation;
447 if( ptwXY1->interpolation == ptwXY_interpo
448 if( ptwXY2->interpolation == ptwXY_interpo
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 ) !=
457 if( ptwX_length( ptwX_norm ) != ngs )
458 *status = nfu_badNorm;
459 return( NULL );
460 }
461 }
462
463 if( ( ff = ptwXY_intersectionWith_ptwX( pt
464 if( ( gg = ptwXY_intersectionWith_ptwX( pt
465 if( ( ff->length == 0 ) || ( gg->length ==
466 ptwXY_free( ff );
467 ptwXY_free( gg );
468 return( ptwX_createLine( ngs, ngs, 0,
469 }
470
471 if( ( *status = ptwXY_tweakDomainsToMutual
472 if( ( f = ptwXY_union( ff, gg, status, ptw
473 if( ( g = ptwXY_union( gg, f, status, ptwX
474
475 if( ( groupedData = ptwX_new( ngs, status
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
485 x2 = f->points[i].x;
486 if( x2 > xg2 ) break;
487 fy2p = fy2 = f->points[i].y;
488 if( f->interpolation == ptwXY_
489 gy2p = gy2 = g->points[i].y;
490 if( g->interpolation == ptwXY_
491 sum += ( ( fy1 + fy2p ) * ( gy
492 }
493 }
494 if( sum != 0. ) {
495 if( normType == ptwXY_group_normTy
496 sum /= ( xg2 - xg1 ); }
497 else if( normType == ptwXY_group_n
498 if( ptwX_norm->points[igs] ==
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( group
523 return( NULL );
524 }
525 /*
526 **********************************************
527 */
528 ptwXPoints *ptwXY_groupThreeFunctions( ptwXYPo
529 ptwXY_group_normType normType, ptwXPoi
530
531 int64_t i, igs, ngs;
532 double x1, fy1, gy1, hy1, x2, fy2, gy2, hy
533 ptwXYPoints *f = NULL, *ff, *fff = NULL, *
534 ptwXPoints *groupedData = NULL;
535
536 if( ( *status = ptwXY_simpleCoalescePoints
537 if( ( *status = ptwXY_simpleCoalescePoints
538 if( ( *status = ptwXY_simpleCoalescePoints
539 if( ( *status = groupBoundaries->status )
540 *status = nfu_otherInterpolation;
541 if( ptwXY1->interpolation == ptwXY_interpo
542 if( ptwXY2->interpolation == ptwXY_interpo
543 if( ptwXY3->interpolation == ptwXY_interpo
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 ) !=
552 if( ptwX_length( ptwX_norm ) != ngs )
553 *status = nfu_badNorm;
554 return( NULL );
555 }
556 }
557
558 if( ( ff = ptwXY_intersectionWith_ptwX( pt
559 if( ( gg = ptwXY_intersectionWith_ptwX( pt
560 if( ( hh = ptwXY_intersectionWith_ptwX( pt
561 if( ( ff->length == 0 ) || ( gg->length ==
562
563 if( ( *status = ptwXY_tweakDomainsToMutual
564 if( ( *status = ptwXY_tweakDomainsToMutual
565 if( ( *status = ptwXY_tweakDomainsToMutual
566 if( ( fff = ptwXY_union( ff, gg, status,
567 if( ( h = ptwXY_union( hh, fff, status,
568 if( ( f = ptwXY_union( fff, h, status,
569 if( ( g = ptwXY_union( gg, h, status,
570
571 if( ( groupedData = ptwX_new( ngs, status
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
582 x2 = f->points[i].x;
583 if( x2 > xg2 ) break;
584 fy2p = fy2 = f->points[i].y;
585 if( f->interpolation == ptwXY_
586 gy2p = gy2 = g->points[i].y;
587 if( g->interpolation == ptwXY_
588 hy2p = hy2 = h->points[i].y;
589 if( h->interpolation == ptwXY_
590 sum += ( ( fy1 + fy2p ) * ( gy
591 }
592 }
593 if( sum != 0. ) {
594 if( normType == ptwXY_group_normTy
595 sum /= ( xg2 - xg1 ); }
596 else if( normType == ptwXY_group_n
597 if( ptwX_norm->points[igs] ==
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( group
627 return( NULL );
628 }
629 /*
630 **********************************************
631 */
632 ptwXPoints *ptwXY_runningIntegral( ptwXYPoints
633
634 int i;
635 ptwXPoints *runningIntegral = NULL;
636 double integral = 0., sum;
637
638 if( ( *status = ptwXY_simpleCoalescePoints
639 if( ( runningIntegral = ptwX_new( ptwXY->l
640
641 if( ( *status = ptwX_setPointAtIndex( runn
642 for( i = 1; i < ptwXY->length; i++ ) {
643 if( ( *status = ptwXY_f_integrate( ptw
644 ptwXY->points[i].x, ptwXY->points[
645 integral += sum;
646 if( ( *status = ptwX_setPointAtIndex(
647 }
648 return( runningIntegral );
649
650 err:
651 if( runningIntegral != NULL ) ptwX_free( r
652 return( NULL );
653 }
654 /*
655 **********************************************
656 */
657 double ptwXY_integrateWithFunction( ptwXYPoint
658 double xMin, double xMax, int degree,
659
660 int64_t i1, i2, n1 = ptwXY->length;
661 long evaluations;
662 double integral = 0., integral_, sign = -1
663 ptwXY_integrateWithFunctionInfo integrateW
664 ptwXYPoint *point;
665
666 if( ( *status = ptwXY->status ) != nfu_Oka
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
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 ) bre
684 }
685 for( i2 = n1 - 1; i2 > i1; i2-- ) {
686 if( ptwXY->points[i2-1].x < xMax ) bre
687 }
688 point = &(ptwXY->points[i1]);
689
690 integrateWithFunctionInfo.degree = degree;
691 integrateWithFunctionInfo.func = func;
692 integrateWithFunctionInfo.argList = argLis
693 integrateWithFunctionInfo.interpolation =
694 integrateWithFunctionInfo.x2 = point->x;
695 integrateWithFunctionInfo.y2 = point->y;
696
697 xa = xMin;
698 for( ; i1 < i2; i1++ ) {
699 integrateWithFunctionInfo.x1 = integra
700 integrateWithFunctionInfo.y1 = integra
701 ++point;
702 integrateWithFunctionInfo.x2 = point->
703 integrateWithFunctionInfo.y2 = point->
704 xb = point->x;
705 if( xb > xMax ) xb = xMax;
706 *status = nf_GnG_adaptiveQuadrature( p
707 xa, xb, recursionLimit, tolerance,
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
719 double x2, double *integral ) {
720
721 ptwXY_integrateWithFunctionInfo *integrate
722 nfu_status status;
723
724 status = nf_Legendre_GaussianQuadrature( i
725 return( status );
726 }
727 /*
728 **********************************************
729 */
730 static nfu_status ptwXY_integrateWithFunction3
731
732 double yf;
733 ptwXY_integrateWithFunctionInfo *integrate
734 nfu_status status;
735
736 if( ( status = ptwXY_interpolatePoint( int
737 integrateWithFunctionInfo->x1, int
738 integrateWithFunctionInfo->x2, int
739 status = integrateWithFunctionInfo->fu
740 *y *= yf;
741 }
742 return( status );
743 }
744
745 #if defined __cplusplus
746 }
747 #endif
748