src/numafunc1.c
l_int32
numaIntegrateInterval(NUMA       *nax,
                      NUMA       *nay,
                      l_float32   x0,
                      l_float32   x1,
                      l_int32     npts,
                      l_float32  *psum)
{
l_int32     i, nx, ny;
l_float32   minx, maxx, sum, del;
l_float32  *fay;
NUMA       *naiy;

    PROCNAME("numaIntegrateInterval");

    if (!psum)
        return ERROR_INT("&sum not defined", procName, 1);
    *psum = 0.0;
    if (!nay)
        return ERROR_INT("nay not defined", procName, 1);
    if (!nax)
        return ERROR_INT("nax not defined", procName, 1);
    if (x0 > x1)
        return ERROR_INT("x0 > x1", procName, 1);
    if (npts < 2)
        return ERROR_INT("npts < 2", procName, 1);
    ny = numaGetCount(nay);
    nx = numaGetCount(nax);
    if (nx != ny)
        return ERROR_INT("nax and nay not same size arrays", procName, 1);
    if (ny < 2)
        return ERROR_INT("not enough points", procName, 1);
    numaGetMin(nax, &minx, NULL);
    numaGetMax(nax, &maxx, NULL);
    if (x0 < minx || x1 > maxx)
        return ERROR_INT("xval is out of bounds", procName, 1);

        /* Generate interpolated array over specified interval */
    if (numaInterpolateArbxInterval(nax, nay, L_LINEAR_INTERP, x0, x1,
                                    npts, NULL, &naiy))
        return ERROR_INT("interpolation failed", procName, 1);

    del = (x1 - x0) / ((l_float32)npts - 1.0);
    fay = numaGetFArray(naiy, L_NOCOPY);

        /* Compute integral (simple trapezoid) */
    sum = 0.5 * (fay[0] + fay[npts - 1]);
    for (i = 1; i < npts - 1; i++)
        sum += fay[i];
    *psum = del * sum;

    numaDestroy(&naiy);
    return 0;
}
/*!
 *  numaIntegrateInterval()
 *
 *      Input:  nax (numa of abscissa values)
 *              nay (numa of ordinate values, corresponding to nax)
 *              x0 (start value of interval)
 *              x1 (end value of interval)
 *              npts (number of points to evaluate function in interval)
 *              &sum (<return> integral of function over interval)
 *      Return: 0 if OK, 1 on error (e.g., if x0 or x1 is outside range)
 *
 *  Notes:
 *      (1) The values in nax must be sorted in increasing order.
 *          If they are not sorted, it is done in the interpolation
 *          step, and a warning is issued.
 *      (2) Caller should check for valid return.
 */


(1) The values in nax must be sorted in increasing order.
If they are not sorted, it is done in the interpolation
step, and a warning is issued.
(2) Caller should check for valid return.
( 1 ) The values in nax must be sorted in increasing order . If they are not sorted it is done in the interpolation step and a warning is issued . ( 2 ) Caller should check for valid return .
l_int32
numaIntegrateInterval(NUMA       *nax,
                      NUMA       *nay,
                      l_float32   x0,
                      l_float32   x1,
                      l_int32     npts,
                      l_float32  *psum)
{
l_int32     i, nx, ny;
l_float32   minx, maxx, sum, del;
l_float32  *fay;
NUMA       *naiy;
    PROCNAME("numaIntegrateInterval");
    if (!psum)
        return ERROR_INT("&sum not defined", procName, 1);
    *psum = 0.0;
    if (!nay)
        return ERROR_INT("nay not defined", procName, 1);
    if (!nax)
        return ERROR_INT("nax not defined", procName, 1);
    if (x0 > x1)
        return ERROR_INT("x0 > x1", procName, 1);
    if (npts < 2)
        return ERROR_INT("npts < 2", procName, 1);
    ny = numaGetCount(nay);
    nx = numaGetCount(nax);
    if (nx != ny)
        return ERROR_INT("nax and nay not same size arrays", procName, 1);
    if (ny < 2)
        return ERROR_INT("not enough points", procName, 1);
    numaGetMin(nax, &minx, NULL);
    numaGetMax(nax, &maxx, NULL);
    if (x0 < minx || x1 > maxx)
        return ERROR_INT("xval is out of bounds", procName, 1);
    if (numaInterpolateArbxInterval(nax, nay, L_LINEAR_INTERP, x0, x1,
                                    npts, NULL, &naiy))
        return ERROR_INT("interpolation failed", procName, 1);
    del = (x1 - x0) / ((l_float32)npts - 1.0);
    fay = numaGetFArray(naiy, L_NOCOPY);
    sum = 0.5 * (fay[0] + fay[npts - 1]);
    for (i = 1; i < npts - 1; i++)
        sum += fay[i];
    *psum = del * sum;
    numaDestroy(&naiy);
    return 0;
}
l_int32 numaIntegrateInterval ( NUMA * nax , NUMA * nay , l_float32 x0 , l_float32 x1 , l_int32 npts , l_float32 * psum ) { l_int32 i , nx , ny ; l_float32 minx , maxx , sum , del ; l_float32 * fay ; NUMA * naiy ; PROCNAME ( " " ) ; if ( ! psum ) return ERROR_INT ( " " , procName , 1 ) ; * psum = 0.0 ; if ( ! nay ) return ERROR_INT ( " " , procName , 1 ) ; if ( ! nax ) return ERROR_INT ( " " , procName , 1 ) ; if ( x0 > x1 ) return ERROR_INT ( " " , procName , 1 ) ; if ( npts < 2 ) return ERROR_INT ( " " , procName , 1 ) ; ny = numaGetCount ( nay ) ; nx = numaGetCount ( nax ) ; if ( nx != ny ) return ERROR_INT ( " " , procName , 1 ) ; if ( ny < 2 ) return ERROR_INT ( " " , procName , 1 ) ; numaGetMin ( nax , & minx , NULL ) ; numaGetMax ( nax , & maxx , NULL ) ; if ( x0 < minx || x1 > maxx ) return ERROR_INT ( " " , procName , 1 ) ; if ( numaInterpolateArbxInterval ( nax , nay , L_LINEAR_INTERP , x0 , x1 , npts , NULL , & naiy ) ) return ERROR_INT ( " " , procName , 1 ) ; del = ( x1 - x0 ) / ( ( l_float32 ) npts - 1.0 ) ; fay = numaGetFArray ( naiy , L_NOCOPY ) ; sum = 0.5 * ( fay [ 0 ] + fay [ npts - 1 ] ) ; for ( i = 1 ; i < npts - 1 ; i ++ ) sum += fay [ i ] ; * psum = del * sum ; numaDestroy ( & naiy ) ; return 0 ; }
numaIntegrateInterval()
      Input:  nax (numa of abscissa values)
              nay (numa of ordinate values, corresponding to nax)
              x0 (start value of interval)
              x1 (end value of interval)
              npts (number of points to evaluate function in interval)
              &sum (<return> integral of function over interval)
      Return: 0 if OK, 1 on error (e.g., if x0 or x1 is outside range)
numaIntegrateInterval () Input : nax ( numa of abscissa values ) nay ( numa of ordinate values corresponding to nax ) x0 ( start value of interval ) x1 ( end value of interval ) npts ( number of points to evaluate function in interval ) &sum ( <return > integral of function over interval ) Return : 0 if OK 1 on error ( e . g . if x0 or x1 is outside range )
/* Generate interpolated array over specified interval */ /* Compute integral (simple trapezoid) */