ptwXY_interpolation.cc File Reference

#include <cmath>
#include <float.h>
#include "ptwXY.h"

Go to the source code of this file.

Macros
#define	minEps   5e-16

Typedefs
typedef nfu_status(*	interpolation_func) (ptwXYPoints *desc, double x1, double y1, double x2, double y2, int depth)

Functions
nfu_status	ptwXY_interpolatePoint (ptwXY_interpolation interpolation, double x, double *y, double x1, double y1, double x2, double y2)
ptwXYPoints *	ptwXY_flatInterpolationToLinear (ptwXYPoints *ptwXY, double lowerEps, double upperEps, nfu_status *status)
ptwXYPoints *	ptwXY_toOtherInterpolation (ptwXYPoints *ptwXY, ptwXY_interpolation interpolationTo, double accuracy, nfu_status *status)
ptwXYPoints *	ptwXY_toUnitbase (ptwXYPoints *ptwXY, nfu_status *status)
ptwXYPoints *	ptwXY_fromUnitbase (ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status)
ptwXYPoints *	ptwXY_unitbaseInterpolate (double w, double w1, ptwXYPoints *ptwXY1, double w2, ptwXYPoints *ptwXY2, nfu_status *status)

Macro Definition Documentation

minEps

#define minEps   5e-16

Typedef Documentation

interpolation_func

typedef nfu_status(* interpolation_func) (ptwXYPoints *desc, double x1, double y1, double x2, double y2, int depth)

Definition at line 19 of file ptwXY_interpolation.cc.

Function Documentation

ptwXY_flatInterpolationToLinear()

ptwXYPoints * ptwXY_flatInterpolationToLinear	(	ptwXYPoints *	ptwXY,
		double	lowerEps,
		double	upperEps,
		nfu_status *	status )

Definition at line 74 of file ptwXY_interpolation.cc.

                                                                                                                         {
 
    int64_t i, length;
    double x;
    ptwXYPoints *n;
    ptwXYPoint *p1 = NULL, *p2 = NULL, *p3;
 
#define minEps 5e-16
 
    if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( NULL );
    *status = nfu_invalidInterpolation;
    if( ptwXY->interpolation != ptwXY_interpolationFlat ) return( NULL );
    *status = nfu_badInput;
    if( ( lowerEps < 0 ) || ( upperEps < 0 ) || ( ( lowerEps == 0 ) && ( upperEps == 0 ) ) ) return( NULL );
    if( ( lowerEps != 0 ) && ( lowerEps < minEps ) ) lowerEps = minEps;
    if( ( upperEps != 0 ) && ( upperEps < minEps ) ) upperEps = minEps;
 
    length = ptwXY->length * ( 1 + ( lowerEps == 0 ? 0 : 1 ) + ( lowerEps == 0 ? 0 : 1 ) );
    if( ( n = ptwXY_new( ptwXY_interpolationLinLin, NULL, ptwXY->biSectionMax, ptwXY->accuracy, length, ptwXY->overflowLength, status, ptwXY->userFlag ) ) == NULL ) return( NULL );
 
    p3 = ptwXY->points;
    if( ptwXY->length > 0 ) ptwXY_setValueAtX( n, p3->x, p3->y );
    for( i = 0; i < ptwXY->length; i++, p3++ ) {
        if( i > 1 ) {
            if( lowerEps > 0 ) {
                x = ptwXY_flatInterpolationToLinear_eps( p2->x, -lowerEps );
                if( x > p1->x ) {
                    if( ( *status = ptwXY_setValueAtX( n, x, p1->y ) ) != nfu_Okay ) goto Err;
                }
            }
            if( lowerEps == 0 ) if( ( *status = ptwXY_setValueAtX( n, p2->x, p1->y ) ) != nfu_Okay ) goto Err;
            if( upperEps == 0 ) if( ( *status = ptwXY_setValueAtX( n, p2->x, p2->y ) ) != nfu_Okay ) goto Err;
            if( upperEps > 0 ) {
                x = ptwXY_flatInterpolationToLinear_eps( p2->x, upperEps );
                if( x < p3->x ) {
                    if( ( *status = ptwXY_setValueAtX( n, x, p2->y ) ) != nfu_Okay ) goto Err;
                }
            }
        }
        p1 = p2;
        p2 = p3;
    }
    if( ptwXY->length > 1 ) {
        if( ( lowerEps != 0 ) && ( p1->y != p2->y ) ) {
            x = ptwXY_flatInterpolationToLinear_eps( p2->x, -lowerEps );
            if( x > p1->x ) {
                if( ( *status = ptwXY_setValueAtX( n, x, p1->y ) ) != nfu_Okay ) goto Err;
            }
        }
        if( ( *status = ptwXY_setValueAtX( n, p2->x, p2->y ) ) != nfu_Okay ) goto Err;
    }
 
    return( n );
 
Err:
    ptwXY_free( n );
    return( NULL );
 
#undef minEps
}

ptwXY_fromUnitbase()

ptwXYPoints * ptwXY_fromUnitbase	(	ptwXYPoints *	ptwXY,
		double	xMin,
		double	xMax,
		nfu_status *	status )

Definition at line 331 of file ptwXY_interpolation.cc.

                                                                                                    {
 
    int64_t i, length;
    ptwXYPoints *n;
    ptwXYPoint *p, *p2;
    double dx, inverseDx, xLast = 0.;
 
    *status = nfu_tooFewPoints;
    if( ptwXY->length < 2 ) return( NULL );
    if( ( n = ptwXY_clone( ptwXY, status ) ) == NULL ) return( NULL );
 
    dx = xMax - xMin;
    inverseDx = 1. / dx;
    length = n->length;
    for( i = 0, p2 = p = n->points; i < length; ++i, ++p ) {
        p2->x = p->x * dx + xMin;
        if( i > 0 ) {
            if( std::fabs( p2->x - xLast ) <= 10. * DBL_EPSILON * ( std::fabs( p2->x ) + std::fabs( xLast ) ) ) {
                --(n->length);
                continue;
            }
        }
        p2->y = p->y * inverseDx;
        xLast = p2->x;
        ++p2;
    }
    n->points[n->length-1].x = xMax;                          /* Make sure last point is realy xMax. */
    return( n );
}

Referenced by ptwXY_unitbaseInterpolate().

ptwXY_interpolatePoint()

nfu_status ptwXY_interpolatePoint	(	ptwXY_interpolation	interpolation,
		double	x,
		double *	y,
		double	x1,
		double	y1,
		double	x2,
		double	y2 )

Definition at line 30 of file ptwXY_interpolation.cc.

                                                                                                                                        {
 
    nfu_status status = nfu_Okay;
 
    if( interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
    if( ( x1 > x2 ) || ( x < x1 ) || ( x > x2 ) ) return( nfu_invalidInterpolation );
    if( y1 == y2 ) {
        *y = y1; }
    else if( x1 == x2 ) {
        *y = 0.5 * ( y1  + y2 ); }
    else if( x == x1 ) {
        *y = y1; }
    else if( x == x2 ) {
        *y = y2; }
    else {
        switch( interpolation ) {
        case ptwXY_interpolationLinLin :
            *y = ( y1 * ( x2 - x ) + y2 * ( x - x1 ) ) / ( x2 - x1 );
            break;
        case ptwXY_interpolationLogLin :
            if( ( x <= 0. ) || ( x1 <= 0. ) || ( x2 <= 0. ) ) return( nfu_invalidInterpolation );
            *y = ( y1 * G4Log( x2 / x ) + y2 * G4Log( x / x1 ) ) / G4Log( x2 / x1 );
            break;
        case ptwXY_interpolationLinLog :
            if( ( y1 <= 0. ) || ( y2 <= 0. ) ) return( nfu_invalidInterpolation );
            *y = G4Exp( ( G4Log( y1 ) * ( x2 - x ) + G4Log( y2 ) * ( x - x1 ) ) / ( x2 - x1 ) );
            break;
        case ptwXY_interpolationLogLog :
            if( ( x <= 0. ) || ( x1 <= 0. ) || ( x2 <= 0. ) ) return( nfu_invalidInterpolation );
            if( ( y1 <= 0. ) || ( y2 <= 0. ) ) return( nfu_invalidInterpolation );
            *y = G4Exp( ( G4Log( y1 ) * G4Log( x2 / x ) + G4Log( y2 ) * G4Log( x / x1 ) ) / G4Log( x2 / x1 ) );
            break;
        case ptwXY_interpolationFlat :
            *y = y1;
            break;
        default :
            status = nfu_invalidInterpolation;
        }
    }
    return( status );
}

Referenced by ptwXY_dullEdges(), ptwXY_getValueAtX(), ptwXY_integrate(), ptwXY_integrateWithWeight_sqrt_x(), ptwXY_integrateWithWeight_x(), ptwXY_thicken(), and ptwXY_union().

ptwXY_toOtherInterpolation()

ptwXYPoints * ptwXY_toOtherInterpolation	(	ptwXYPoints *	ptwXY,
		ptwXY_interpolation	interpolationTo,
		double	accuracy,
		nfu_status *	status )

Definition at line 153 of file ptwXY_interpolation.cc.

                                                                                                                                        {
/*
*   This function only works when 'ptwXY->interpolation == interpolationTo' or when interpolationTo is ptwXY_interpolationLinLin.
*/
    ptwXYPoints *n1;
    interpolation_func func = NULL;
 
    if( ( *status = ptwXY->status ) != nfu_Okay ) return( NULL );
    if( ptwXY->interpolation == interpolationTo ) {
        *status = nfu_Okay;
        return( ptwXY_clone( ptwXY, status ) ); }
    else {
        if( interpolationTo == ptwXY_interpolationLinLin ) {
            switch( ptwXY->interpolation ) {
            case ptwXY_interpolationLogLog :
                func = ptwXY_LogLogToLinLin; break;
            case ptwXY_interpolationLinLog :
                func = ptwXY_LinLogToLinLin; break;
            case ptwXY_interpolationLogLin :
                func = ptwXY_LogLinToLinLin; break;
            case ptwXY_interpolationOther :
                if( ptwXY->interpolationOtherInfo.getValueFunc != NULL ) func = ptwXY_otherToLinLin;
                break;
            case ptwXY_interpolationLinLin :        /* Stops compilers from complaining. */
            case ptwXY_interpolationFlat :
                break;
            }
        }
    }
    *status = nfu_unsupportedInterpolationConversion;
    if( func == NULL ) return( NULL );
 
    *status = nfu_Okay;
    if( ( n1 = ptwXY_cloneToInterpolation( ptwXY, interpolationTo, status ) ) == NULL ) return( NULL );
    if( accuracy < ptwXY->accuracy ) accuracy = ptwXY->accuracy;
    n1->accuracy = accuracy;
 
    n1->interpolationOtherInfo.getValueFunc = ptwXY->interpolationOtherInfo.getValueFunc;
    n1->interpolationOtherInfo.argList = ptwXY->interpolationOtherInfo.argList;
    *status = ptwXY_toOtherInterpolation2( n1, ptwXY, func );
    n1->interpolationOtherInfo.getValueFunc = NULL;
    n1->interpolationOtherInfo.argList = NULL;
    if( *status != nfu_Okay ) n1 = ptwXY_free( n1 );
    return( n1 );
}

ptwXY_toUnitbase()

ptwXYPoints * ptwXY_toUnitbase	(	ptwXYPoints *	ptwXY,
		nfu_status *	status )

Definition at line 306 of file ptwXY_interpolation.cc.

                                                                        {
 
    int64_t i;
    ptwXYPoints *n;
    ptwXYPoint *p;
    double xMin, xMax, dx, inverseDx;
 
    *status = nfu_tooFewPoints;
    if( ptwXY->length < 2 ) return( NULL );
    if( ( n = ptwXY_clone( ptwXY, status ) ) == NULL ) return( NULL );
 
    xMin = n->points[0].x;
    xMax = n->points[n->length-1].x;
    dx = xMax - xMin;
    inverseDx = 1. / dx;
    for( i = 0, p = n->points; i < n->length; i++, p++ ) {
        p->x = ( p->x - xMin ) * inverseDx;
        p->y = p->y * dx;
    }
    n->points[n->length-1].x = 1.;                          /* Make sure last point is realy 1. */
    return( n );
}

Referenced by ptwXY_unitbaseInterpolate().

ptwXY_unitbaseInterpolate()

ptwXYPoints * ptwXY_unitbaseInterpolate	(	double	w,
		double	w1,
		ptwXYPoints *	ptwXY1,
		double	w2,
		ptwXYPoints *	ptwXY2,
		nfu_status *	status )

Definition at line 363 of file ptwXY_interpolation.cc.

                                                                                                                                       {
/* 
*   Should we not be checking the interpolation members???????
*/
    int64_t i;
    ptwXYPoints *n1 = NULL, *n2 = NULL, *a = NULL, *r = NULL;
    ptwXYPoint *p;
    double f, g, xMin, xMax;
 
    *status = nfu_XOutsideDomain;
    if( w <= w1 ) {
        if( w < w1 ) return( NULL );
        return( ptwXY_clone( ptwXY1, status ) );
    }
    if( w >= w2 ) {
        if( w > w2 ) return( NULL );
        return( ptwXY_clone( ptwXY2, status ) );
    }
    if( ( n1 = ptwXY_toUnitbase( ptwXY1, status ) ) == NULL ) return( NULL );
    if( ( n2 = ptwXY_toUnitbase( ptwXY2, status ) ) == NULL ) goto Err;
    f = ( w - w1 ) / ( w2 - w1 );
    g = 1. - f;
    for( i = 0, p = n1->points; i < n1->length; i++, p++ ) p->y *= g;
    for( i = 0, p = n2->points; i < n2->length; i++, p++ ) p->y *= f;
    if( ( a = ptwXY_add_ptwXY( n1, n2, status ) ) == NULL ) goto Err;
 
    xMin = g * ptwXY1->points[0].x + f * ptwXY2->points[0].x;
    xMax = g * ptwXY1->points[ptwXY1->length-1].x + f * ptwXY2->points[ptwXY2->length-1].x;
    if( ( r = ptwXY_fromUnitbase( a, xMin, xMax, status ) ) == NULL ) goto Err;
    ptwXY_free( n1 );
    ptwXY_free( n2 );
    ptwXY_free( a );
    return( r );
 
Err:
    if( n1 != NULL ) ptwXY_free( n1 );
    if( n2 != NULL ) ptwXY_free( n2 );
    if( a  != NULL ) ptwXY_free( a );
    return( NULL );
}