ptwXY_functions.cc File Reference

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

Go to the source code of this file.

Functions
nfu_status	ptwXY_pow (ptwXYPoints *ptwXY, double v)
nfu_status	ptwXY_exp (ptwXYPoints *ptwXY, double a)
ptwXYPoints *	ptwXY_convolution (ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int mode)

Function Documentation

ptwXY_convolution()

ptwXYPoints * ptwXY_convolution	(	ptwXYPoints *	ptwXY1,
		ptwXYPoints *	ptwXY2,
		nfu_status *	status,
		int	mode )

Definition at line 96 of file ptwXY_functions.cc.

                                                                                                         {
/*
*   Currently, only supports linear-linear interpolation.
*
*   This function calculates        c(y) = integral dx f1(x) * f2(y-x)
*
*/
    int64_t i1, i2, n1, n2, n;
    ptwXYPoints *f1 = ptwXY1, *f2 = ptwXY2, *convolute;
    double accuracy = ptwXY1->accuracy, yMin, yMax, c, y, dy;
 
    if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
    if( ( *status = ptwXY_simpleCoalescePoints( ptwXY2 ) ) != nfu_Okay ) return( NULL );
 
    *status = nfu_unsupportedInterpolation;
    if( ( ptwXY1->interpolation != ptwXY_interpolationLinLin ) || ( ptwXY2->interpolation != ptwXY_interpolationLinLin ) ) return( NULL );
    *status = nfu_Okay;
 
    n1 = f1->length;
    n2 = f2->length;
 
    if( ( n1 == 0 ) || ( n2 == 0 ) ) {
        convolute = ptwXY_new( ptwXY_interpolationLinLin, NULL, 1., accuracy, 0, 0, status, 0 );
        return( convolute );
    }
 
    if( ( n1 == 1 ) || ( n2 == 1 ) ) {
        *status = nfu_tooFewPoints;
        return( NULL );
    }
 
    if( accuracy < ptwXY2->accuracy ) accuracy = ptwXY2->accuracy;
    n = n1 * n2;
    if( mode == 0 ) {
        mode = 1;
        if( n > 1000 ) mode = -1;
    }
    if( n > 100000 ) mode = -1;
    if( ( convolute = ptwXY_new( ptwXY_interpolationLinLin, NULL, 1., accuracy, 400, 40, status, 0 ) ) == NULL ) return( NULL );
 
    yMin = f1->points[0].x + f2->points[0].x;
    yMax = f1->points[n1 - 1].x + f2->points[n2 - 1].x;
 
    if( ( *status = ptwXY_setValueAtX( convolute, yMin, 0. ) ) != nfu_Okay ) goto Err;
 
    if( mode < 0 ) {
        dy = ( yMax - yMin ) / 400;
        for( y = yMin + dy; y < yMax; y += dy ) {
            if( ( *status = ptwXY_convolution2( f1, f2, y, yMin, &c ) ) != nfu_Okay ) goto Err;
            if( ( *status = ptwXY_setValueAtX( convolute, y, c ) ) != nfu_Okay ) goto Err;
        } }
    else {
        for( i1 = 0; i1 < n1; i1++ ) {
            for( i2 = 0; i2 < n2; i2++ ) {
                y = yMin + ( f1->points[i1].x - f1->points[0].x ) + ( f2->points[i2].x - f2->points[0].x );
                if( y <= yMin ) continue;
                if( y >= yMax ) continue;
                if( ( *status = ptwXY_convolution2( f1, f2, y, yMin, &c ) ) != nfu_Okay ) goto Err;
                if( ( *status = ptwXY_setValueAtX( convolute, y, c ) ) != nfu_Okay ) goto Err;
            }
        }
    }
    if( ( *status = ptwXY_setValueAtX( convolute, yMax, 0. ) ) != nfu_Okay ) goto Err;
    if( ( *status = ptwXY_simpleCoalescePoints( convolute ) ) != nfu_Okay ) goto Err;
    for( i1 = convolute->length - 1; i1 > 0; i1-- ) {
        if( ( *status = ptwXY_convolution3( convolute, f1, f2, convolute->points[i1 - 1].x, convolute->points[i1 - 1].y, 
            convolute->points[i1].x, convolute->points[i1].y, yMin ) ) != nfu_Okay ) goto Err;
    }
 
    return( convolute );
 
Err:
    ptwXY_free( convolute );
    return( NULL );
}

ptwXY_exp()

nfu_status ptwXY_exp	(	ptwXYPoints *	ptwXY,
		double	a )

Definition at line 43 of file ptwXY_functions.cc.

                                                     { 
 
    int64_t i, length;
    nfu_status status;
    double x1, y1, z1, x2, y2, z2;
 
    length = ptwXY->length;
    if( length < 1 ) return( ptwXY->status );
    if( ptwXY->interpolation == ptwXY_interpolationFlat ) return( nfu_invalidInterpolation );
    if( ptwXY->interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
 
    if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status );
    x2 = ptwXY->points[length-1].x;
    y2 = a * ptwXY->points[length-1].y;
    z2 = ptwXY->points[length-1].y = G4Exp( y2 );
    for( i = length - 2; i >= 0; i-- ) {
        x1 = ptwXY->points[i].x;
        y1 = a * ptwXY->points[i].y;
        z1 = ptwXY->points[i].y = G4Exp( y1 );
        if( ( status = ptwXY_exp_s( ptwXY, x1, y1, z1, x2, y2, z2, 0 ) ) != nfu_Okay ) return( status );
        x2 = x1;
        y2 = y1;
    }
    return( status );
}

ptwXY_pow()

nfu_status ptwXY_pow	(	ptwXYPoints *	ptwXY,
		double	v )

Definition at line 24 of file ptwXY_functions.cc.

                                                     { 
 
    return( ptwXY_applyFunction( ptwXY, ptwXY_pow_callback, (void *) &v, 0 ) );
}