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