ptwXY.h File Reference

#include <stdio.h>
#include <stdint.h>
#include <nf_utilities.h>
#include <ptwX.h>

Go to the source code of this file.

Classes
struct	ptwXYPoint_s
struct	ptwXY_interpolationOtherInfo
struct	ptwXYOverflowPoint_s
struct	ptwXYPoints_s

Macros
#define	ptwXY_minimumSize   10 /* This must be > 0 otherwise some logic will fail. */
#define	ptwXY_minimumOverflowSize   4 /* This must be > 0 otherwise some logic will fail. */
#define	ptwXY_maxBiSectionMax   20
#define	ptwXY_minAccuracy   1e-14
#define	ptwXY_sectionSubdivideMax   1 << 16
#define	ClosestAllowXFactor   10
#define	ptwXY_union_fill   1 /* If filling, union is filled with y value of first ptw. */
#define	ptwXY_union_trim   2 /* If trimming, union in only over common domain of ptw1 and ptw2. */
#define	ptwXY_union_mergeClosePoints   4 /* If true, union calls ptwXY_mergeClosePoints with eps = 4 * DBL_EPSILON. */

Typedefs
typedef enum ptwXY_dataFrom_e	ptwXY_dataFrom
typedef enum ptwXY_group_normType_e	ptwXY_group_normType
typedef enum ptwXY_sigma_e	ptwXY_sigma
typedef enum ptwXY_interpolation_e	ptwXY_interpolation
typedef enum ptwXY_lessEqualGreaterX_e	ptwXY_lessEqualGreaterX
typedef struct ptwXYPoint_s	ptwXYPoint
typedef nfu_status(*	ptwXY_createFromFunction_callback) (double x, double *y, void *argList)
typedef nfu_status(*	ptwXY_applyFunction_callback) (ptwXYPoint *point, void *argList)
typedef nfu_status(*	ptwXY_getValue_callback) (void *argList, double x, double *y, double x1, double y1, double x2, double y2)
typedef struct ptwXYOverflowPoint_s	ptwXYOverflowPoint
typedef struct ptwXYPoints_s	ptwXYPoints

Enumerations
enum	ptwXY_dataFrom_e { ptwXY_dataFrom_Unknown , ptwXY_dataFrom_Points , ptwXY_dataFrom_Overflow }
enum	ptwXY_group_normType_e { ptwXY_group_normType_none , ptwXY_group_normType_dx , ptwXY_group_normType_norm }
enum	ptwXY_sigma_e { ptwXY_sigma_none , ptwXY_sigma_plusMinus , ptwXY_sigma_Minus , ptwXY_sigma_plus }
enum	ptwXY_interpolation_e {   ptwXY_interpolationLinLin , ptwXY_interpolationLinLog , ptwXY_interpolationLogLin , ptwXY_interpolationLogLog ,   ptwXY_interpolationFlat , ptwXY_interpolationOther }
enum	ptwXY_lessEqualGreaterX_e {   ptwXY_lessEqualGreaterX_empty , ptwXY_lessEqualGreaterX_lessThan , ptwXY_lessEqualGreaterX_equal , ptwXY_lessEqualGreaterX_between ,   ptwXY_lessEqualGreaterX_greater }

Functions
ptwXYPoints *	ptwXY_new (ptwXY_interpolation interpolation, ptwXY_interpolationOtherInfo const *interpolationOtherInfo, double biSectionMax, double accuracy, int64_t primarySize, int64_t secondarySize, nfu_status *status, int userFlag)
nfu_status	ptwXY_setup (ptwXYPoints *ptwXY, ptwXY_interpolation interpolation, ptwXY_interpolationOtherInfo const *interpolationOtherInfo, double biSectionMax, double accuracy, int64_t primarySize, int64_t secondarySize, int userFlag)
ptwXYPoints *	ptwXY_create (ptwXY_interpolation interpolation, ptwXY_interpolationOtherInfo const *interpolationOtherInfo, double biSectionMax, double accuracy, int64_t primarySize, int64_t secondarySize, int64_t length, double const *xy, nfu_status *status, int userFlag)
ptwXYPoints *	ptwXY_createFrom_Xs_Ys (ptwXY_interpolation interpolation, ptwXY_interpolationOtherInfo const *interpolationOtherInfo, double biSectionMax, double accuracy, int64_t primarySize, int64_t secondarySize, int64_t length, double const *Xs, double const *Ys, nfu_status *status, int userFlag)
nfu_status	ptwXY_copy (ptwXYPoints *dest, ptwXYPoints *src)
ptwXYPoints *	ptwXY_clone (ptwXYPoints *ptwXY, nfu_status *status)
ptwXYPoints *	ptwXY_cloneToInterpolation (ptwXYPoints *ptwXY, ptwXY_interpolation interpolationTo, nfu_status *status)
ptwXYPoints *	ptwXY_slice (ptwXYPoints *ptwXY, int64_t index1, int64_t index2, int64_t secondarySize, nfu_status *status)
ptwXYPoints *	ptwXY_xSlice (ptwXYPoints *ptwXY, double xMin, double xMax, int64_t secondarySize, int fill, nfu_status *status)
ptwXYPoints *	ptwXY_xMinSlice (ptwXYPoints *ptwXY, double xMin, int64_t secondarySize, int fill, nfu_status *status)
ptwXYPoints *	ptwXY_xMaxSlice (ptwXYPoints *ptwXY, double xMax, int64_t secondarySize, int fill, nfu_status *status)
ptwXY_interpolation	ptwXY_getInterpolation (ptwXYPoints *ptwXY)
char const *	ptwXY_getInterpolationString (ptwXYPoints *ptwXY)
nfu_status	ptwXY_getStatus (ptwXYPoints *ptwXY)
int	ptwXY_getUserFlag (ptwXYPoints *ptwXY)
void	ptwXY_setUserFlag (ptwXYPoints *ptwXY, int userFlag)
double	ptwXY_getAccuracy (ptwXYPoints *ptwXY)
double	ptwXY_setAccuracy (ptwXYPoints *ptwXY, double accuracy)
double	ptwXY_getBiSectionMax (ptwXYPoints *ptwXY)
double	ptwXY_setBiSectionMax (ptwXYPoints *ptwXY, double biSectionMax)
nfu_status	ptwXY_reallocatePoints (ptwXYPoints *ptwXY, int64_t size, int forceSmallerResize)
nfu_status	ptwXY_reallocateOverflowPoints (ptwXYPoints *ptwXY, int64_t size)
nfu_status	ptwXY_coalescePoints (ptwXYPoints *ptwXY, int64_t size, ptwXYPoint *newPoint, int forceSmallerResize)
nfu_status	ptwXY_simpleCoalescePoints (ptwXYPoints *ptwXY)
nfu_status	ptwXY_clear (ptwXYPoints *ptwXY)
nfu_status	ptwXY_release (ptwXYPoints *ptwXY)
ptwXYPoints *	ptwXY_free (ptwXYPoints *ptwXY)
int64_t	ptwXY_length (ptwXYPoints *ptwXY)
int64_t	ptwXY_getNonOverflowLength (ptwXYPoints const *ptwXY)
nfu_status	ptwXY_setXYData (ptwXYPoints *ptwXY, int64_t length, double const *xy)
nfu_status	ptwXY_setXYDataFromXsAndYs (ptwXYPoints *ptwXY, int64_t length, double const *x, double const *y)
nfu_status	ptwXY_deletePoints (ptwXYPoints *ptwXY, int64_t i1, int64_t i2)
ptwXYPoint *	ptwXY_getPointAtIndex (ptwXYPoints *ptwXY, int64_t index)
ptwXYPoint *	ptwXY_getPointAtIndex_Unsafely (ptwXYPoints *ptwXY, int64_t index)
nfu_status	ptwXY_getXYPairAtIndex (ptwXYPoints *ptwXY, int64_t index, double *x, double *y)
ptwXY_lessEqualGreaterX	ptwXY_getPointsAroundX (ptwXYPoints *ptwXY, double x, ptwXYOverflowPoint *lessThanEqualXPoint, ptwXYOverflowPoint *greaterThanXPoint)
ptwXY_lessEqualGreaterX	ptwXY_getPointsAroundX_closeIsEqual (ptwXYPoints *ptwXY, double x, ptwXYOverflowPoint *lessThanEqualXPoint, ptwXYOverflowPoint *greaterThanXPoint, double eps, int *closeIsEqual, ptwXYPoint **closePoint)
nfu_status	ptwXY_getValueAtX (ptwXYPoints *ptwXY, double x, double *y)
nfu_status	ptwXY_setValueAtX (ptwXYPoints *ptwXY, double x, double y)
nfu_status	ptwXY_setValueAtX_overrideIfClose (ptwXYPoints *ptwXY, double x, double y, double eps, int override)
nfu_status	ptwXY_mergeFromXsAndYs (ptwXYPoints *ptwXY, int length, double *xs, double *ys)
nfu_status	ptwXY_mergeFromXYs (ptwXYPoints *ptwXY, int length, double *xys)
nfu_status	ptwXY_appendXY (ptwXYPoints *ptwXY, double x, double y)
nfu_status	ptwXY_setXYPairAtIndex (ptwXYPoints *ptwXY, int64_t index, double x, double y)
nfu_status	ptwXY_getSlopeAtX (ptwXYPoints *ptwXY, double x, const char side, double *slope)
double	ptwXY_getXMinAndFrom (ptwXYPoints *ptwXY, ptwXY_dataFrom *dataFrom)
double	ptwXY_getXMin (ptwXYPoints *ptwXY)
double	ptwXY_getXMaxAndFrom (ptwXYPoints *ptwXY, ptwXY_dataFrom *dataFrom)
double	ptwXY_getXMax (ptwXYPoints *ptwXY)
double	ptwXY_getYMin (ptwXYPoints *ptwXY)
double	ptwXY_getYMax (ptwXYPoints *ptwXY)
nfu_status	ptwXY_clip (ptwXYPoints *ptwXY1, double yMin, double yMax)
nfu_status	ptwXY_thicken (ptwXYPoints *ptwXY1, int sectionSubdivideMax, double dxMax, double fxMax)
ptwXYPoints *	ptwXY_thin (ptwXYPoints *ptwXY1, double accuracy, nfu_status *status)
nfu_status	ptwXY_trim (ptwXYPoints *ptwXY)
ptwXYPoints *	ptwXY_union (ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int unionOptions)
nfu_status	ptwXY_scaleOffsetXAndY (ptwXYPoints *ptwXY, double xScale, double xOffset, double yScale, double yOffset)
nfu_status	ptwXY_abs (ptwXYPoints *ptwXY)
nfu_status	ptwXY_neg (ptwXYPoints *ptwXY)
nfu_status	ptwXY_slopeOffset (ptwXYPoints *ptwXY, double slope, double offset)
nfu_status	ptwXY_add_double (ptwXYPoints *ptwXY, double value)
nfu_status	ptwXY_sub_doubleFrom (ptwXYPoints *ptwXY, double value)
nfu_status	ptwXY_sub_fromDouble (ptwXYPoints *ptwXY, double value)
nfu_status	ptwXY_mul_double (ptwXYPoints *ptwXY, double value)
nfu_status	ptwXY_div_doubleFrom (ptwXYPoints *ptwXY, double value)
nfu_status	ptwXY_div_fromDouble (ptwXYPoints *ptwXY, double value)
nfu_status	ptwXY_mod (ptwXYPoints *ptwXY, double m, int pythonMod)
ptwXYPoints *	ptwXY_binary_ptwXY (ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, double v1, double v2, double v1v2, nfu_status *status)
ptwXYPoints *	ptwXY_add_ptwXY (ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status)
ptwXYPoints *	ptwXY_sub_ptwXY (ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status)
ptwXYPoints *	ptwXY_mul_ptwXY (ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status)
ptwXYPoints *	ptwXY_mul2_ptwXY (ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status)
ptwXYPoints *	ptwXY_div_ptwXY (ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int safeDivide)
nfu_status	ptwXY_pow (ptwXYPoints *ptwXY, double p)
nfu_status	ptwXY_exp (ptwXYPoints *ptwXY, double a)
ptwXYPoints *	ptwXY_convolution (ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int mode)
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 interpolation, double accuracy, nfu_status *status)
ptwXYPoints *	ptwXY_unitbaseInterpolate (double w, double w1, ptwXYPoints *ptwXY1, double w2, ptwXYPoints *ptwXY2, nfu_status *status)
ptwXYPoints *	ptwXY_toUnitbase (ptwXYPoints *ptwXY, nfu_status *status)
ptwXYPoints *	ptwXY_fromUnitbase (ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status)
ptwXPoints *	ptwXY_getXArray (ptwXYPoints *ptwXY, nfu_status *status)
nfu_status	ptwXY_dullEdges (ptwXYPoints *ptwXY, double lowerEps, double upperEps, int positiveXOnly)
nfu_status	ptwXY_mergeClosePoints (ptwXYPoints *ptwXY, double epsilon)
ptwXYPoints *	ptwXY_intersectionWith_ptwX (ptwXYPoints *ptwXY, ptwXPoints *ptwX, nfu_status *status)
nfu_status	ptwXY_areDomainsMutual (ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2)
nfu_status	ptwXY_tweakDomainsToMutualify (ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, int epsilonFactor, double epsilon)
nfu_status	ptwXY_mutualifyDomains (ptwXYPoints *ptwXY1, double lowerEps1, double upperEps1, int positiveXOnly1, ptwXYPoints *ptwXY2, double lowerEps2, double upperEps2, int positiveXOnly2)
nfu_status	ptwXY_copyToC_XY (ptwXYPoints *ptwXY, int64_t index1, int64_t index2, int64_t allocatedSize, int64_t *numberOfPoints, double *xy)
nfu_status	ptwXY_valueTo_ptwXAndY (ptwXYPoints *ptwXY, double **xs, double **ys)
ptwXYPoints *	ptwXY_valueTo_ptwXY (double x1, double x2, double y, nfu_status *status)
ptwXYPoints *	ptwXY_createGaussianCenteredSigma1 (double accuracy, nfu_status *status)
ptwXYPoints *	ptwXY_createGaussian (double accuracy, double xCenter, double sigma, double amplitude, double xMin, double xMax, double dullEps, nfu_status *status)
void	ptwXY_update_biSectionMax (ptwXYPoints *ptwXY1, double oldLength)
ptwXYPoints *	ptwXY_createFromFunction (int n, double *xs, ptwXY_createFromFunction_callback func, void *argList, double accuracy, int checkForRoots, int biSectionMax, nfu_status *status)
ptwXYPoints *	ptwXY_createFromFunction2 (ptwXPoints *xs, ptwXY_createFromFunction_callback func, void *argList, double accuracy, int checkForRoots, int biSectionMax, nfu_status *status)
nfu_status	ptwXY_applyFunction (ptwXYPoints *ptwXY1, ptwXY_applyFunction_callback func, void *argList, int checkForRoots)
ptwXYPoints *	ptwXY_fromString (char const *str, ptwXY_interpolation interpolation, ptwXY_interpolationOtherInfo const *interpolationOtherInfo, double biSectionMax, double accuracy, char **endCharacter, nfu_status *status)
void	ptwXY_showInteralStructure (ptwXYPoints *ptwXY, FILE *f, int printPointersAsNull)
void	ptwXY_simpleWrite (ptwXYPoints *ptwXY, FILE *f, char *format)
void	ptwXY_simplePrint (ptwXYPoints *ptwXY, char *format)
nfu_status	ptwXY_f_integrate (ptwXY_interpolation interpolation, double x1, double y1, double x2, double y2, double *value)
double	ptwXY_integrate (ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status)
double	ptwXY_integrateDomain (ptwXYPoints *ptwXY, nfu_status *status)
nfu_status	ptwXY_normalize (ptwXYPoints *ptwXY1)
double	ptwXY_integrateDomainWithWeight_x (ptwXYPoints *ptwXY, nfu_status *status)
double	ptwXY_integrateWithWeight_x (ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status)
double	ptwXY_integrateDomainWithWeight_sqrt_x (ptwXYPoints *ptwXY, nfu_status *status)
double	ptwXY_integrateWithWeight_sqrt_x (ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status)
ptwXPoints *	ptwXY_groupOneFunction (ptwXYPoints *ptwXY, ptwXPoints *groupBoundaries, ptwXY_group_normType normType, ptwXPoints *ptwX_norm, nfu_status *status)
ptwXPoints *	ptwXY_groupTwoFunctions (ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, ptwXPoints *groupBoundaries, ptwXY_group_normType normType, ptwXPoints *ptwX_norm, nfu_status *status)
ptwXPoints *	ptwXY_groupThreeFunctions (ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, ptwXYPoints *ptwXY3, ptwXPoints *groupBoundaries, ptwXY_group_normType normType, ptwXPoints *ptwX_norm, nfu_status *status)
ptwXPoints *	ptwXY_runningIntegral (ptwXYPoints *ptwXY, nfu_status *status)
double	ptwXY_integrateWithFunction (ptwXYPoints *ptwXY, ptwXY_createFromFunction_callback func, void *argList, double xMin, double xMax, int degree, int recursionLimit, double tolerance, nfu_status *status)

Macro Definition Documentation

ClosestAllowXFactor

#define ClosestAllowXFactor   10

Definition at line 25 of file ptwXY.h.

Referenced by ptwXY_createFromFunction().

ptwXY_maxBiSectionMax

#define ptwXY_maxBiSectionMax   20

Definition at line 22 of file ptwXY.h.

Referenced by ptwXY_setBiSectionMax(), ptwXY_setup(), ptwXY_update_biSectionMax(), and ptwXY_valueTo_ptwXY().

ptwXY_minAccuracy

#define ptwXY_minAccuracy   1e-14

Definition at line 23 of file ptwXY.h.

Referenced by ptwXY_setAccuracy(), ptwXY_setup(), and ptwXY_valueTo_ptwXY().

ptwXY_minimumOverflowSize

#define ptwXY_minimumOverflowSize   4 /* This must be > 0 otherwise some logic will fail. */

Definition at line 21 of file ptwXY.h.

Referenced by ptwXY_reallocateOverflowPoints().

ptwXY_minimumSize

#define ptwXY_minimumSize   10 /* This must be > 0 otherwise some logic will fail. */

Definition at line 20 of file ptwXY.h.

Referenced by ptwXY_reallocatePoints(), and ptwXY_union().

ptwXY_sectionSubdivideMax

#define ptwXY_sectionSubdivideMax   1 << 16

Definition at line 24 of file ptwXY.h.

Referenced by ptwXY_thicken().

ptwXY_union_fill

#define ptwXY_union_fill   1 /* If filling, union is filled with y value of first ptw. */

Definition at line 31 of file ptwXY.h.

Referenced by ptwXY_binary_ptwXY(), ptwXY_div_ptwXY(), ptwXY_groupThreeFunctions(), ptwXY_groupTwoFunctions(), and ptwXY_union().

ptwXY_union_mergeClosePoints

#define ptwXY_union_mergeClosePoints   4 /* If true, union calls ptwXY_mergeClosePoints with eps = 4 * DBL_EPSILON. */

Definition at line 33 of file ptwXY.h.

Referenced by ptwXY_binary_ptwXY(), ptwXY_div_ptwXY(), and ptwXY_union().

ptwXY_union_trim

#define ptwXY_union_trim   2 /* If trimming, union in only over common domain of ptw1 and ptw2. */

Definition at line 32 of file ptwXY.h.

Referenced by ptwXY_union().

Typedef Documentation

ptwXY_applyFunction_callback

typedef nfu_status(* ptwXY_applyFunction_callback) (ptwXYPoint *point, void *argList)

Definition at line 66 of file ptwXY.h.

ptwXY_createFromFunction_callback

typedef nfu_status(* ptwXY_createFromFunction_callback) (double x, double *y, void *argList)

Definition at line 65 of file ptwXY.h.

ptwXY_dataFrom

typedef enum ptwXY_dataFrom_e ptwXY_dataFrom

ptwXY_getValue_callback

typedef nfu_status(* ptwXY_getValue_callback) (void *argList, double x, double *y, double x1, double y1, double x2, double y2)

Definition at line 67 of file ptwXY.h.

ptwXY_group_normType

typedef enum ptwXY_group_normType_e ptwXY_group_normType

ptwXY_interpolation

typedef enum ptwXY_interpolation_e ptwXY_interpolation

ptwXY_lessEqualGreaterX

typedef enum ptwXY_lessEqualGreaterX_e ptwXY_lessEqualGreaterX

ptwXY_sigma

typedef enum ptwXY_sigma_e ptwXY_sigma

ptwXYOverflowPoint

typedef struct ptwXYOverflowPoint_s ptwXYOverflowPoint

ptwXYPoint

typedef struct ptwXYPoint_s ptwXYPoint

ptwXYPoints

typedef struct ptwXYPoints_s ptwXYPoints

Enumeration Type Documentation

ptwXY_dataFrom_e

enum ptwXY_dataFrom_e

Enumerator
ptwXY_dataFrom_Unknown
ptwXY_dataFrom_Points
ptwXY_dataFrom_Overflow

Definition at line 27 of file ptwXY.h.

{ ptwXY_dataFrom_Unknown, ptwXY_dataFrom_Points, ptwXY_dataFrom_Overflow } ptwXY_dataFrom;

ptwXY_group_normType_e

enum ptwXY_group_normType_e

Enumerator
ptwXY_group_normType_none
ptwXY_group_normType_dx
ptwXY_group_normType_norm

Definition at line 28 of file ptwXY.h.

{ ptwXY_group_normType_none, ptwXY_group_normType_dx, ptwXY_group_normType_norm } ptwXY_group_normType;

ptwXY_interpolation_e

enum ptwXY_interpolation_e

Enumerator
ptwXY_interpolationLinLin
ptwXY_interpolationLinLog
ptwXY_interpolationLogLin
ptwXY_interpolationLogLog
ptwXY_interpolationFlat
ptwXY_interpolationOther

Definition at line 35 of file ptwXY.h.

                                   { ptwXY_interpolationLinLin, ptwXY_interpolationLinLog, ptwXY_interpolationLogLin, ptwXY_interpolationLogLog,
    ptwXY_interpolationFlat, ptwXY_interpolationOther } ptwXY_interpolation;

ptwXY_lessEqualGreaterX_e

enum ptwXY_lessEqualGreaterX_e

Enumerator
ptwXY_lessEqualGreaterX_empty
ptwXY_lessEqualGreaterX_lessThan
ptwXY_lessEqualGreaterX_equal
ptwXY_lessEqualGreaterX_between
ptwXY_lessEqualGreaterX_greater

Definition at line 57 of file ptwXY.h.

                                       { ptwXY_lessEqualGreaterX_empty, ptwXY_lessEqualGreaterX_lessThan, ptwXY_lessEqualGreaterX_equal,
    ptwXY_lessEqualGreaterX_between, ptwXY_lessEqualGreaterX_greater } ptwXY_lessEqualGreaterX;

ptwXY_sigma_e

enum ptwXY_sigma_e

Enumerator
ptwXY_sigma_none
ptwXY_sigma_plusMinus
ptwXY_sigma_Minus
ptwXY_sigma_plus

Definition at line 34 of file ptwXY.h.

{ ptwXY_sigma_none, ptwXY_sigma_plusMinus, ptwXY_sigma_Minus, ptwXY_sigma_plus } ptwXY_sigma;

Function Documentation

ptwXY_abs()

nfu_status ptwXY_abs

(

ptwXYPoints *

ptwXY

)

Definition at line 19 of file ptwXY_unitaryOperators.cc.

                                           {
 
    int64_t i, nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
    ptwXYPoint *p;
    ptwXYOverflowPoint *o, *overflowHeader = &(ptwXY->overflowHeader);
 
    if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
 
    for( i = 0, p = ptwXY->points; i < nonOverflowLength; i++, p++ ) p->y = std::fabs( p->y );
    for( o = overflowHeader->next; o != overflowHeader; o = o->next ) o->point.y = std::fabs( o->point.y );
    return( ptwXY->status );
}

ptwXY_add_double()

nfu_status ptwXY_add_double	(	ptwXYPoints *	ptwXY,
		double	value )

Definition at line 40 of file ptwXY_binaryOperators.cc.

{ return( ptwXY_slopeOffset( ptwXY, 1., value ) ); }

ptwXY_add_ptwXY()

ptwXYPoints * ptwXY_add_ptwXY	(	ptwXYPoints *	ptwXY1,
		ptwXYPoints *	ptwXY2,
		nfu_status *	status )

Definition at line 138 of file ptwXY_binaryOperators.cc.

                                                                                             {
 
    ptwXYPoints *sum;
 
    if( ptwXY1->length == 0 ) {
        sum = ptwXY_clone( ptwXY2, status ); }
    else if( ptwXY2->length == 0 ) {
        sum = ptwXY_clone( ptwXY1, status ); }
    else {
        sum = ptwXY_binary_ptwXY( ptwXY1, ptwXY2, 1., 1., 0., status );
    }
    return( sum );
}

Referenced by MCGIDI_target_heated_read(), and ptwXY_unitbaseInterpolate().

ptwXY_appendXY()

nfu_status ptwXY_appendXY	(	ptwXYPoints *	ptwXY,
		double	x,
		double	y )

Definition at line 1062 of file ptwXY_core.cc.

                                                                    {
 
    int64_t nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
    ptwXY_dataFrom dataFrom;
 
    if( ptwXY->length != 0 ) {
        double xMax = ptwXY_getXMaxAndFrom( ptwXY, &dataFrom );
        if( xMax >= x ) return( nfu_XNotAscending );
    }
 
    if( nonOverflowLength < ptwXY->allocatedSize ) {      /* Room at end of points. Also handles the case when length = 0. */
        ptwXY->points[nonOverflowLength].x = x;
        ptwXY->points[nonOverflowLength].y = y; }
    else {
        if( ptwXY->overflowLength == ptwXY->overflowAllocatedSize ) {
            ptwXYPoint newPoint = { x, y };
            return( ptwXY_coalescePoints( ptwXY, ptwXY->length + ptwXY->overflowAllocatedSize, &newPoint, 0 ) ); }
        else {                                              /* Add to end of overflow. */
            ptwXYOverflowPoint *overflowPoint = &(ptwXY->overflowPoints[ptwXY->overflowLength]);
 
            overflowPoint->prior = ptwXY->overflowHeader.prior;
            overflowPoint->next = overflowPoint->prior->next;
            overflowPoint->index = ptwXY->length;
            overflowPoint->prior->next = overflowPoint;
            overflowPoint->next->prior = overflowPoint;
            overflowPoint->point.x = x;
            overflowPoint->point.y = y;
            ptwXY->overflowLength++;
        }
    }
    ptwXY->length++;
    return( nfu_Okay );
}

ptwXY_applyFunction()

nfu_status ptwXY_applyFunction	(	ptwXYPoints *	ptwXY1,
		ptwXY_applyFunction_callback	func,
		void *	argList,
		int	checkForRoots )

Definition at line 146 of file ptwXY_misc.cc.

                                                                                                                           {
 
    int64_t i, originalLength = ptwXY1->length, notFirstPass = 0;
    double y1, y2 = 0;
    nfu_status status;
    ptwXYPoint p1, p2;
 
    checkForRoots = checkForRoots && ptwXY1->biSectionMax;
    if( ptwXY1->status != nfu_Okay ) return( ptwXY1->status );
    if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
    if( ptwXY1->interpolation == ptwXY_interpolationFlat ) return( nfu_invalidInterpolation );
    if( ( status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( status );
    for( i = originalLength - 1; i >= 0; i-- ) {
        y1 = ptwXY1->points[i].y;
        if( ( status = func( &(ptwXY1->points[i]), argList ) ) != nfu_Okay ) return( status );
        p1 = ptwXY1->points[i];
        if( notFirstPass ) {
            if( ( status = ptwXY_applyFunction2( ptwXY1, y1, y2, &p1, &p2, func, argList, 0, checkForRoots ) ) != nfu_Okay ) return( status );
        }
        notFirstPass = 1;
        p2 = p1;
        y2 = y1;
    }
    ptwXY_update_biSectionMax( ptwXY1, (double) originalLength );
    return( status );
}

Referenced by ptwXY_pow().

ptwXY_areDomainsMutual()

nfu_status ptwXY_areDomainsMutual	(	ptwXYPoints *	ptwXY1,
		ptwXYPoints *	ptwXY2 )

Definition at line 257 of file ptwXY_convenient.cc.

                                                                              {
 
    nfu_status status;
    int64_t n1 = ptwXY1->length, n2 = ptwXY2->length;
    ptwXYPoint *xy1, *xy2;
 
    if( ( status = ptwXY1->status ) != nfu_Okay ) return( status );
    if( ( status = ptwXY2->status ) != nfu_Okay ) return( status );
    if( n1 == 0 ) return( nfu_empty );
    if( n2 == 0 ) return( nfu_empty );
    if( n1 < 2 ) { 
        status = nfu_tooFewPoints; }
    else if( n2 < 2 ) { 
        status = nfu_tooFewPoints; }
    else {
        xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, 0 );
        xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, 0 );
        if( xy1->x < xy2->x ) {
            if( xy2->y != 0. ) status = nfu_domainsNotMutual; }
        else if( xy1->x > xy2->x ) {
            if( xy1->y != 0. ) status = nfu_domainsNotMutual;
        }
 
        if( status == nfu_Okay ) {
            xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, n1 - 1 );
            xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, n2 - 1 );
            if( xy1->x < xy2->x ) {
                if( xy1->y != 0. ) status = nfu_domainsNotMutual; }
            else if( xy1->x > xy2->x ) {
                if( xy2->y != 0. ) status = nfu_domainsNotMutual;
            }
        }
    }
    return( status );
}

Referenced by ptwXY_binary_ptwXY(), ptwXY_div_ptwXY(), and ptwXY_mutualifyDomains().

ptwXY_binary_ptwXY()

ptwXYPoints * ptwXY_binary_ptwXY	(	ptwXYPoints *	ptwXY1,
		ptwXYPoints *	ptwXY2,
		double	v1,
		double	v2,
		double	v1v2,
		nfu_status *	status )

Definition at line 108 of file ptwXY_binaryOperators.cc.

                                                                                                                                   {
 
    int64_t i;
    int unionOptions = ptwXY_union_fill | ptwXY_union_mergeClosePoints;
    double y;
    ptwXYPoints *n;
    ptwXYPoint *p;
 
    *status = nfu_otherInterpolation;
    if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
    if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( NULL );
    if( ( *status = ptwXY_areDomainsMutual( ptwXY1, ptwXY2 ) ) != nfu_Okay ) return( NULL );
    if( ( ptwXY1->interpolation == ptwXY_interpolationFlat ) || ( ptwXY2->interpolation == ptwXY_interpolationFlat ) ) {
        *status = nfu_invalidInterpolation;
        if( ( ptwXY1->interpolation != ptwXY2->interpolation ) ) return( NULL );
    }
    if( ( n = ptwXY_union( ptwXY1, ptwXY2, status, unionOptions ) ) != NULL ) {
        for( i = 0, p = n->points; i < n->length; i++, p++ ) {
            if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, p->x, &y ) ) != nfu_Okay ) goto Err;
            p->y = v1 * p->y + v2 * y + v1v2 * y * p->y;
        }
    }
    return( n );
Err:
    if( n ) ptwXY_free( n );
    return( NULL );
}

Referenced by ptwXY_add_ptwXY(), ptwXY_mul_ptwXY(), and ptwXY_sub_ptwXY().

ptwXY_clear()

nfu_status ptwXY_clear

(

ptwXYPoints *

ptwXY

)

Definition at line 536 of file ptwXY_core.cc.

                                             {
 
    if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
 
    ptwXY->length = 0;
    ptwXY->overflowLength = 0;
    ptwXY->overflowHeader.prior = &(ptwXY->overflowHeader);
    ptwXY->overflowHeader.next = &(ptwXY->overflowHeader);
    return( nfu_Okay );
}

Referenced by ptwXY_clip(), ptwXY_copy(), and ptwXY_setXYDataFromXsAndYs().

ptwXY_clip()

nfu_status ptwXY_clip	(	ptwXYPoints *	ptwXY1,
		double	yMin,
		double	yMax )

Definition at line 25 of file ptwXY_methods.cc.

                                                                       {
/*
    This function acts oddly for xy = [ [ 1, 0 ], [ 3, -2 ], [ 4, 1 ] ] and yMin = 0.2, why???????
    This function probably only works for linear, linear interpolation (mainly because of ptwXY_clip2).
*/
    int64_t i, j, n;
    double x2, y2;
    nfu_status status;
    ptwXYPoints *clipped;
    ptwXYPoint *points;
 
    if( ( status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( status );
    if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
    n = ptwXY1->length;
    if( n > 0 ) {
        i = 0;
        if( ptwXY_getYMax( ptwXY1 ) < yMin ) i = 1;
        if( ptwXY_getYMin( ptwXY1 ) > yMax ) i = 1;
        if( i == 1 ) return( ptwXY_clear( ptwXY1 ) );
    }
    if( n == 1 ) {
        y2 = ptwXY1->points[0].y;
        if( y2 < yMin ) {
            ptwXY1->points[0].y = yMin; }
        else if( y2 > yMax ) {
            ptwXY1->points[0].y = yMax;
        } }
    else if( n > 1 ) {
        if( ( clipped = ptwXY_new( ptwXY1->interpolation, &(ptwXY1->interpolationOtherInfo), 
                ptwXY1->biSectionMax, ptwXY1->accuracy, n, 10, &status, ptwXY1->userFlag ) ) == NULL )
            return( ptwXY1->status = status );
        for( i = 0; i < n; i++ ) {
            x2 = ptwXY1->points[i].x;
            y2 = ptwXY1->points[i].y;
            if( y2 < yMin ) {
                if( i > 0 ) {
                    points = ptwXY_getPointAtIndex_Unsafely( clipped, clipped->length - 1 );
                    if( points->y > yMin ) {
                        if( ( status = ptwXY_clip2( clipped, yMin, points->x, points->y, x2, y2 ) ) != nfu_Okay ) goto Err;
                    }
                }
                if( ( status = ptwXY_setValueAtX( clipped, x2, yMin ) ) != nfu_Okay ) goto Err;
                j = i;
                for( i++; i < n; i++ ) if( !( ptwXY1->points[i].y < yMin ) ) break;
                if( i < n ) {
                    x2 = ptwXY1->points[i].x;
                    y2 = ptwXY1->points[i].y;
                    if( ( status = ptwXY_clip2( clipped, yMin, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) goto Err;
                    if( y2 > yMax ) {
                        if( ( status = ptwXY_clip2( clipped, yMax, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) goto Err;
                    } }
                else if( j != n - 1 ) {
                    if( ( status = ptwXY_setValueAtX( clipped, ptwXY1->points[n - 1].x, yMin ) ) != nfu_Okay ) goto Err;
                }
                i--; }
            else if( y2 > yMax ) {
                if( i > 0 ) {
                    points = ptwXY_getPointAtIndex_Unsafely( clipped, clipped->length - 1 );
                    if( points->y < yMax ) {
                        if( ( status = ptwXY_clip2( clipped, yMax, points->x, points->y, x2, y2 ) ) != nfu_Okay ) goto Err;
                    }
                }
                if( ( status = ptwXY_setValueAtX( clipped, x2, yMax ) ) != nfu_Okay ) goto Err;
                j = i;
                for( i++; i < n; i++ ) if( !( ptwXY1->points[i].y > yMax ) ) break;
                if( i < n ) {
                    x2 = ptwXY1->points[i].x;
                    y2 = ptwXY1->points[i].y;
                    if( ( status = ptwXY_clip2( clipped, yMax, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) goto Err;
                    if( y2 < yMin ) {
                        if( ( status = ptwXY_clip2( clipped, yMin, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) goto Err;
                    } }
                else if( j != n - 1 ) {
                    if( ( status = ptwXY_setValueAtX( clipped, ptwXY1->points[n - 1].x, yMax ) ) != nfu_Okay ) goto Err;
                }
                i--; }
            else {
                if( ( status = ptwXY_setValueAtX( clipped, x2, y2 ) ) != nfu_Okay ) goto Err;
            }
        }
        if( ( status = ptwXY_simpleCoalescePoints( clipped ) ) != nfu_Okay ) goto Err;
        ptwXY1->length = clipped->length;   /* The squeamish may want to skip the next few lines. */
        clipped->length = n;
        n = ptwXY1->allocatedSize;
        ptwXY1->allocatedSize = clipped->allocatedSize;
        clipped->allocatedSize = n;
        points = clipped->points;
        clipped->points = ptwXY1->points;
        ptwXY1->points = points;
        ptwXY_free( clipped );
    }
 
    return( ptwXY1->status );
 
Err:
    ptwXY_free( clipped );
    return( ptwXY1->status = status );
}

ptwXY_clone()

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

Definition at line 208 of file ptwXY_core.cc.

                                                                   {
 
    return( ptwXY_slice( ptwXY, 0, ptwXY->length, ptwXY->overflowAllocatedSize, status ) );
}

Referenced by GIDI_settings_processedFlux::GIDI_settings_processedFlux(), GIDI_settings_processedFlux::operator=(), ptwXY_add_ptwXY(), ptwXY_cloneToInterpolation(), ptwXY_fromUnitbase(), ptwXY_intersectionWith_ptwX(), ptwXY_mul_ptwXY(), ptwXY_sub_ptwXY(), ptwXY_thin(), ptwXY_toOtherInterpolation(), ptwXY_toUnitbase(), ptwXY_unitbaseInterpolate(), and ptwXY_xSlice().

ptwXY_cloneToInterpolation()

ptwXYPoints * ptwXY_cloneToInterpolation	(	ptwXYPoints *	ptwXY,
		ptwXY_interpolation	interpolationTo,
		nfu_status *	status )

Definition at line 215 of file ptwXY_core.cc.

                                                                                                                       {
 
    ptwXYPoints *n1;
 
    if( interpolationTo == ptwXY_interpolationOther ) {
        *status = nfu_otherInterpolation;
        return( NULL );
    }
    if( ( n1 = ptwXY_clone( ptwXY, status ) ) != NULL ) {
        if( n1->interpolation == ptwXY_interpolationOther ) nfu_free( (void *) n1->interpolationOtherInfo.interpolationString );
        n1->interpolation = interpolationTo;
        switch( interpolationTo ) {
            case ptwXY_interpolationLinLin :
                n1->interpolationOtherInfo.interpolationString = linLinInterpolationString; break;
            case ptwXY_interpolationLinLog :
                n1->interpolationOtherInfo.interpolationString = linLogInterpolationString; break;
            case ptwXY_interpolationLogLin :
                n1->interpolationOtherInfo.interpolationString = logLinInterpolationString; break;
            case ptwXY_interpolationLogLog :
                n1->interpolationOtherInfo.interpolationString = logLogInterpolationString; break;
            case ptwXY_interpolationFlat :
                n1->interpolationOtherInfo.interpolationString = flatInterpolationString; break;
            case ptwXY_interpolationOther :     /* Does not happen, but needed to stop compilers from complaining. */
                break;
        }
        n1->interpolationOtherInfo.getValueFunc = NULL;
        n1->interpolationOtherInfo.argList = NULL;
    }
    return( n1 );
}

Referenced by ptwXY_toOtherInterpolation().

ptwXY_coalescePoints()

nfu_status ptwXY_coalescePoints	(	ptwXYPoints *	ptwXY,
		int64_t	size,
		ptwXYPoint *	newPoint,
		int	forceSmallerResize )

Definition at line 469 of file ptwXY_core.cc.

                                                                                                                  {
 
    int addNewPoint;
    int64_t length = ptwXY->length + ( ( newPoint != NULL ) ? 1 : 0 );
    ptwXYOverflowPoint *last = ptwXY->overflowHeader.prior;
    ptwXYPoint *pointsFrom, *pointsTo;
 
    if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
    if( ptwXY->overflowLength == 0 ) return( nfu_Okay );
 
    if( size < length ) size = length;
    if( size > ptwXY->allocatedSize ) {
        if( ptwXY_reallocatePoints( ptwXY, size, forceSmallerResize ) != nfu_Okay ) return( ptwXY->status );
    }
    pointsFrom = &(ptwXY->points[ptwXY_getNonOverflowLength( ptwXY ) - 1]);
    pointsTo = &(ptwXY->points[length - 1]);
    while( last != &(ptwXY->overflowHeader) ) {
        addNewPoint = 0;
        if( newPoint != NULL ) {
            if( ( pointsFrom >= ptwXY->points ) && ( pointsFrom->x > last->point.x ) ) {
                if( newPoint->x > pointsFrom->x ) addNewPoint = 1; }
            else {
                if( newPoint->x > last->point.x ) addNewPoint = 1;
            }
            if( addNewPoint == 1 ) {
                *pointsTo = *newPoint;
                newPoint = NULL;
            }
        }
        if( addNewPoint == 0 ) {
            if( ( pointsFrom >= ptwXY->points ) && ( pointsFrom->x > last->point.x ) ) {
                *pointsTo = *pointsFrom;
                pointsFrom--; }
            else {
                *pointsTo = last->point;
                last = last->prior;
            }
        }
        pointsTo--;
    }  // Loop checking, 11.06.2015, T. Koi
    while( ( newPoint != NULL ) && ( pointsFrom >= ptwXY->points ) ) { 
        if( newPoint->x > pointsFrom->x ) {
            *pointsTo = *newPoint;
             newPoint = NULL; }
         else {
            *pointsTo = *pointsFrom;
            pointsFrom--;
         }
         pointsTo--;
    } // Loop checking, 11.06.2015, T. Koi
    if( newPoint != NULL ) *pointsTo = *newPoint;
    ptwXY->overflowHeader.prior = &(ptwXY->overflowHeader);
    ptwXY->overflowHeader.next = &(ptwXY->overflowHeader);
    ptwXY->length = length;
    ptwXY->overflowLength = 0;
    return( nfu_Okay );
}

Referenced by ptwXY_appendXY(), ptwXY_createGaussianCenteredSigma1(), ptwXY_reallocateOverflowPoints(), ptwXY_setValueAtX_overrideIfClose(), ptwXY_simpleCoalescePoints(), and ptwXY_xSlice().

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_copy()

nfu_status ptwXY_copy	(	ptwXYPoints *	dest,
		ptwXYPoints *	src )

Definition at line 148 of file ptwXY_core.cc.

                                                             {
 
    int64_t i, nonOverflowLength = ptwXY_getNonOverflowLength( src );
    ptwXYPoint *pointFrom, *pointTo;
    ptwXYOverflowPoint *o, *overflowHeader = &(src->overflowHeader);
 
    if( dest->status != nfu_Okay ) return( dest->status );
    if( src->status != nfu_Okay ) return( src->status );
 
    ptwXY_clear( dest );
    if( dest->interpolation == ptwXY_interpolationOther ) {
        if( dest->interpolationOtherInfo.interpolationString != NULL ) {
            dest->interpolationOtherInfo.interpolationString = (char const *) nfu_free( (void *) dest->interpolationOtherInfo.interpolationString );
        }
    }
    dest->interpolation = ptwXY_interpolationLinLin; /* This and prior lines are in case interpolation is 'other' and ptwXY_reallocatePoints fails. */
    if( dest->allocatedSize < src->length ) ptwXY_reallocatePoints( dest, src->length, 0 );
    if( dest->status != nfu_Okay ) return( dest->status );
    dest->interpolation = src->interpolation;
    if( dest->interpolation == ptwXY_interpolationOther ) {
        if( src->interpolationOtherInfo.interpolationString != NULL ) {
            if( ( dest->interpolationOtherInfo.interpolationString = strdup( src->interpolationOtherInfo.interpolationString ) ) == NULL ) 
                return( dest->status = nfu_mallocError );
        } }
    else {
        dest->interpolationOtherInfo.interpolationString = src->interpolationOtherInfo.interpolationString;
    }
    dest->interpolationOtherInfo.getValueFunc = src->interpolationOtherInfo.getValueFunc;
    dest->interpolationOtherInfo.argList = src->interpolationOtherInfo.argList;
    dest->userFlag = src->userFlag;
    dest->biSectionMax = src->biSectionMax;
    dest->accuracy = src->accuracy;
    dest->minFractional_dx = src->minFractional_dx;
    pointFrom = src->points;
    o = src->overflowHeader.next;
    pointTo = dest->points;
    i = 0;
    while( o != overflowHeader ) {
        if( i < nonOverflowLength ) {
            if( pointFrom->x < o->point.x ) {
                *pointTo = *pointFrom;
                i++;
                pointFrom++; }
            else {
                *pointTo = o->point;
                o = o->next;
            } }
        else {
            *pointTo = o->point;
            o = o->next;
        }
        pointTo++;
    } // Loop checking, 11.06.2015, T. Koi
    for( ; i < nonOverflowLength; i++, pointFrom++, pointTo++ ) *pointTo = *pointFrom;
    dest->length = src->length;
    return( dest->status );
}

ptwXY_copyToC_XY()

nfu_status ptwXY_copyToC_XY	(	ptwXYPoints *	ptwXY,
		int64_t	index1,
		int64_t	index2,
		int64_t	allocatedSize,
		int64_t *	numberOfPoints,
		double *	xy )

Definition at line 424 of file ptwXY_convenient.cc.

                                                                                                                                               {
 
    int64_t i;
    double *d = xys;
    nfu_status status;
    ptwXYPoint *pointFrom;
 
    if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
    if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status );
 
    if( index1 < 0 ) index1 = 0;
    if( index2 > ptwXY->length ) index2 = ptwXY->length;
    if( index2 < index1 ) index2 = index1;
    *numberOfPoints = index2 - index1;
    if( allocatedSize < ( index2 - index1 ) ) return( nfu_insufficientMemory );
 
    for( i = index1, pointFrom = ptwXY->points; i < index2; i++, pointFrom++ ) {
        *(d++) = pointFrom->x;
        *(d++) = pointFrom->y;
    }
 
    return( nfu_Okay );
}

ptwXY_create()

ptwXYPoints * ptwXY_create	(	ptwXY_interpolation	interpolation,
		ptwXY_interpolationOtherInfo const *	interpolationOtherInfo,
		double	biSectionMax,
		double	accuracy,
		int64_t	primarySize,
		int64_t	secondarySize,
		int64_t	length,
		double const *	xy,
		nfu_status *	status,
		int	userFlag )

Definition at line 108 of file ptwXY_core.cc.

                                       {
 
    ptwXYPoints *ptwXY;
 
    if( primarySize < length ) primarySize = length;
    if( ( ptwXY = ptwXY_new( interpolation, interpolationOtherInfo, biSectionMax, accuracy, primarySize, 
            secondarySize, status, userFlag ) ) != NULL ) {
        if( ( *status = ptwXY_setXYData( ptwXY, length, xy ) ) != nfu_Okay ) {
            ptwXY = ptwXY_free( ptwXY );
        }
    }
    return( ptwXY );
}

Referenced by MCGIDI_target_heated_read(), and ptwXY_fromString().

ptwXY_createFrom_Xs_Ys()

ptwXYPoints * ptwXY_createFrom_Xs_Ys	(	ptwXY_interpolation	interpolation,
		ptwXY_interpolationOtherInfo const *	interpolationOtherInfo,
		double	biSectionMax,
		double	accuracy,
		int64_t	primarySize,
		int64_t	secondarySize,
		int64_t	length,
		double const *	Xs,
		double const *	Ys,
		nfu_status *	status,
		int	userFlag )

Definition at line 126 of file ptwXY_core.cc.

                                                         {
 
    int i;
    ptwXYPoints *ptwXY;
 
    if( primarySize < length ) primarySize = length;
    if( ( ptwXY = ptwXY_new( interpolation, interpolationOtherInfo, biSectionMax, accuracy, primarySize, 
            secondarySize, status, userFlag ) ) != NULL ) {
        for( i = 0; i < length; i++ ) {
            ptwXY->points[i].x = Xs[i];
            ptwXY->points[i].y = Ys[i];
        }
        ptwXY->length = length;
    }
 
    return( ptwXY );
}

ptwXY_createFromFunction()

ptwXYPoints * ptwXY_createFromFunction	(	int	n,
		double *	xs,
		ptwXY_createFromFunction_callback	func,
		void *	argList,
		double	accuracy,
		int	checkForRoots,
		int	biSectionMax,
		nfu_status *	status )

Definition at line 40 of file ptwXY_misc.cc.

                                           {
 
    int64_t i;
    double x1, y1, x2 = 0., y2, eps = ClosestAllowXFactor * DBL_EPSILON;
    ptwXYPoints *ptwXY;
    ptwXYPoint *p1, *p2;
 
    *status = nfu_Okay;
    if( n < 2 ) { *status = nfu_tooFewPoints; return( NULL ); }
    for( i = 1; i < n; i++ ) {
        if( xs[i-1] >= xs[i] ) *status = nfu_XNotAscending;
    }
    if( *status == nfu_XNotAscending ) return( NULL );
 
    x1 = xs[0];
    if( ( *status = func( x1, &y1, argList ) ) != nfu_Okay ) return( NULL );
    if( ( ptwXY = ptwXY_new( ptwXY_interpolationLinLin, NULL, biSectionMax, accuracy, 500, 50, status, 0 ) ) == NULL ) return( NULL );
    for( i = 1; i < n; i++ ) {
        if( ( *status = ptwXY_setValueAtX_overrideIfClose( ptwXY, x1, y1, eps, 0 ) ) != nfu_Okay ) goto err;
        x2 = xs[i];
        if( ( *status = func( x2, &y2, argList ) ) != nfu_Okay ) goto err;
        if( ( *status = ptwXY_createFromFunctionBisect( ptwXY, x1, y1, x2, y2, func, argList, 0, checkForRoots, eps ) ) != nfu_Okay ) goto err;
        x1 = x2;
        y1 = y2;
    }
    if( ( *status = ptwXY_setValueAtX_overrideIfClose( ptwXY, x2, y2, eps, 1 ) ) != nfu_Okay ) goto err;
 
    if( checkForRoots ) {
        if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) goto err;
        for( i = ptwXY->length - 1, p2 = NULL; i >= 0; i--, p2 = p1 ) { /* Work backward so lower points are still valid if a new point is added. */
            p1 = &(ptwXY->points[i]);
            if( p2 != NULL ) {
                if( ( p1->y * p2->y ) < 0. ) {
                    if( ( *status = ptwXY_createFromFunctionZeroCrossing( ptwXY, p1->x, p1->y, p2->x, p2->y, func, argList, eps ) ) != nfu_Okay ) goto err;
                }
            }
        }
    }
 
    return( ptwXY );
 
err:
    ptwXY_free( ptwXY );
    return( NULL );
}

Referenced by nf_Legendre_to_ptwXY(), and ptwXY_createFromFunction2().

ptwXY_createFromFunction2()

ptwXYPoints * ptwXY_createFromFunction2	(	ptwXPoints *	xs,
		ptwXY_createFromFunction_callback	func,
		void *	argList,
		double	accuracy,
		int	checkForRoots,
		int	biSectionMax,
		nfu_status *	status )

Definition at line 89 of file ptwXY_misc.cc.

                                           {
 
    return( ptwXY_createFromFunction( (int) xs->length, xs->points, func, argList, accuracy, checkForRoots, biSectionMax, status ) );
}

ptwXY_createGaussian()

ptwXYPoints * ptwXY_createGaussian	(	double	accuracy,
		double	xCenter,
		double	sigma,
		double	amplitude,
		double	xMin,
		double	xMax,
		double	dullEps,
		nfu_status *	status )

Definition at line 566 of file ptwXY_convenient.cc.

                                       {
 
    int64_t i;
    ptwXYPoints *gaussian, *sliced;
    ptwXYPoint *point;
 
    if( ( gaussian = ptwXY_createGaussianCenteredSigma1( accuracy, status ) ) == NULL ) return( NULL );
    for( i = 0, point = gaussian->points; i < gaussian->length; i++, point++ ) {
        point->x = point->x * sigma + xCenter;
        point->y *= amplitude;
    }
    if( ( gaussian->points[0].x < xMin ) || ( gaussian->points[gaussian->length - 1].x > xMax ) ) {
        if( ( sliced = ptwXY_xSlice( gaussian, xMin, xMax, 10, 1, status ) ) == NULL ) goto Err;
        ptwXY_free( gaussian );
        gaussian = sliced;
    }
 
    return( gaussian );
 
Err:
    ptwXY_free( gaussian );
    return( NULL );
}

ptwXY_createGaussianCenteredSigma1()

ptwXYPoints * ptwXY_createGaussianCenteredSigma1	(	double	accuracy,
		nfu_status *	status )

Definition at line 492 of file ptwXY_convenient.cc.

                                                                                       {
 
    int64_t i, n;
    ptwXYPoint *pm, *pp;
    double x1, y1, x2, y2, accuracy2, yMin = 1e-10;
    ptwXYPoints *gaussian;
 
    if( accuracy < 1e-5 ) accuracy = 1e-5;
    if( accuracy > 1e-1 ) accuracy = 1e-1;
    if( ( gaussian = ptwXY_new( ptwXY_interpolationLinLin, NULL, 1., accuracy, 200, 100, status, 0 ) ) == NULL ) return( NULL );
    accuracy2 = accuracy = gaussian->accuracy;
    if( accuracy2 > 5e-3 ) accuracy2 = 5e-3;
 
    x1 = -std::sqrt( -2. * G4Log( yMin ) );
    y1 = yMin;
    x2 = -5.2;
    y2 = G4Exp( -0.5 * x2 * x2 );
    if( ( *status = ptwXY_setValueAtX( gaussian, x1, y1 ) ) != nfu_Okay ) goto Err;
    gaussian->accuracy = 20 * accuracy2;
    if( ( *status = ptwXY_createGaussianCenteredSigma1_2( gaussian, x1, y1, x2, y2, 1 ) ) != nfu_Okay ) goto Err;
    x1 = x2;
    y1 = y2;
    x2 = -4.;
    y2 = G4Exp( -0.5 * x2 * x2 );
    gaussian->accuracy = 5 * accuracy2;
    if( ( *status = ptwXY_createGaussianCenteredSigma1_2( gaussian, x1, y1, x2, y2, 1 ) ) != nfu_Okay ) goto Err;
    x1 = x2;
    y1 = y2;
    x2 = -1;
    y2 = G4Exp( -0.5 * x2 * x2 );
    gaussian->accuracy = accuracy;
    if( ( *status = ptwXY_createGaussianCenteredSigma1_2( gaussian, x1, y1, x2, y2, 1 ) ) != nfu_Okay ) goto Err;
    x1 = x2;
    y1 = y2;
    x2 =  0;
    y2 = G4Exp( -0.5 * x2 * x2 );
    if( ( *status = ptwXY_createGaussianCenteredSigma1_2( gaussian, x1, y1, x2, y2, 1 ) ) != nfu_Okay ) goto Err;
 
    n = gaussian->length;
    if( ( *status = ptwXY_coalescePoints( gaussian, 2 * n + 1, NULL, 0 ) ) != nfu_Okay ) goto Err;
    if( ( *status = ptwXY_setValueAtX( gaussian, 0., 1. ) ) != nfu_Okay ) goto Err;
    pp = &(gaussian->points[gaussian->length]);
    for( i = 0, pm = pp - 2; i < n; i++, pp++, pm-- ) {
        *pp = *pm;
        pp->x *= -1;
    }
    gaussian->length = 2 * n + 1;
 
    return( gaussian );
 
Err:
    ptwXY_free( gaussian );
    return( NULL );
}

Referenced by ptwXY_createGaussian().

ptwXY_deletePoints()

nfu_status ptwXY_deletePoints	(	ptwXYPoints *	ptwXY,
		int64_t	i1,
		int64_t	i2 )

Definition at line 660 of file ptwXY_core.cc.

                                                                            {
 
    int64_t n = ptwXY->length - ( i2 - i1 );
 
    if( ( ptwXY->status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( ptwXY->status );
    if( ( i1 < 0 ) || ( i1 > i2 ) || ( i2 > ptwXY->length ) ) return( nfu_badIndex );
    if( i1 != i2 ) {
        for( ; i2 < ptwXY->length; i1++, i2++ ) ptwXY->points[i1] = ptwXY->points[i2];
        ptwXY->length = n;
    }
    return( ptwXY->status );
}

ptwXY_div_doubleFrom()

nfu_status ptwXY_div_doubleFrom	(	ptwXYPoints *	ptwXY,
		double	value )

Definition at line 44 of file ptwXY_binaryOperators.cc.

                                                                    { 
 
    if( value == 0. ) {
        ptwXY->status = nfu_divByZero; }
    else {
        ptwXY_slopeOffset( ptwXY, 1. / value, 0. ); 
    }
    return( ptwXY->status );
}

ptwXY_div_fromDouble()

nfu_status ptwXY_div_fromDouble	(	ptwXYPoints *	ptwXY,
		double	value )

Definition at line 53 of file ptwXY_binaryOperators.cc.

                                                                    {
/*
*   This does not do any infilling and it should?????????
*/
 
    int64_t i, nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
    ptwXYPoint *p;
    ptwXYOverflowPoint *o, *overflowHeader = &(ptwXY->overflowHeader);
 
    if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
    if( ptwXY->interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
 
    for( i = 0, p = ptwXY->points; i < nonOverflowLength; i++, p++ ) if( p->y == 0. ) ptwXY->status = nfu_divByZero;
    for( o = overflowHeader->next; o != overflowHeader; o = o->next ) if( o->point.y == 0. ) ptwXY->status = nfu_divByZero;
    if( ptwXY->status != nfu_divByZero ) {
        for( i = 0, p = ptwXY->points; i < nonOverflowLength; i++, p++ ) p->y = value / p->y;
        for( o = overflowHeader->next; o != overflowHeader; o = o->next ) o->point.y = value / o->point.y; 
    }
    return( ptwXY->status );
}

ptwXY_div_ptwXY()

ptwXYPoints * ptwXY_div_ptwXY	(	ptwXYPoints *	ptwXY1,
		ptwXYPoints *	ptwXY2,
		nfu_status *	status,
		int	safeDivide )

Definition at line 288 of file ptwXY_binaryOperators.cc.

                                                                                                             {
 
    int isNAN1, isNAN2;
    int64_t i, j, k, zeros = 0, length, iYs;
    double x1, x2, y1, y2, u1, u2, v1, v2, y, xz, nan = nfu_getNAN( ), s1, s2;
    ptwXYPoints *n = NULL;
    ptwXYPoint *p;
 
    if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
    if( ( *status = ptwXY_simpleCoalescePoints( ptwXY2 ) ) != nfu_Okay ) return( NULL );
    *status = nfu_otherInterpolation;
    if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
    if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( NULL );
    if( ( ptwXY1->interpolation == ptwXY_interpolationFlat ) || ( ptwXY1->interpolation == ptwXY_interpolationFlat ) )
        return( ptwXY_div_ptwXY_forFlats( ptwXY1, ptwXY2, status, safeDivide ) );
    
    if( ( *status = ptwXY_areDomainsMutual( ptwXY1, ptwXY2 ) ) != nfu_Okay ) return( NULL );
    if( ( n = ptwXY_union( ptwXY1, ptwXY2, status, ptwXY_union_fill | ptwXY_union_mergeClosePoints ) ) != NULL ) {
        for( i = 0, p = n->points; i < n->length; i++, p++ ) {
            if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, p->x, &y ) ) != nfu_Okay ) goto Err;
            if( y == 0. ) {
                if( p->y == 0. ) {
                    iYs = 0;
                    y1 = 0.;
                    y2 = 0.;
                    if( i > 0 ) {
                        if( ( *status = ptwXY_getSlopeAtX( ptwXY1, p->x, '-', &s1 ) ) != nfu_Okay ) {
                            if( *status != nfu_XOutsideDomain ) goto Err;
                            s1 = 0.;
                        }
                        if( ( *status = ptwXY_getSlopeAtX( ptwXY2, p->x, '-', &s2 ) ) != nfu_Okay ) goto Err;
                        if( s2 == 0. ) {
                            y1 = nan; }
                        else {
                            y1 = s1 / s2;
                        } 
                        iYs++;
                    }
                    if( i < ( n->length - 1 ) ) {
                        if( ( *status = ptwXY_getSlopeAtX( ptwXY1, p->x, '+', &s1 ) ) != nfu_Okay ) {
                            if( *status != nfu_XOutsideDomain ) goto Err;
                            s1 = 0.;
                        }
                        if( ( *status = ptwXY_getSlopeAtX( ptwXY2, p->x, '+', &s2 ) ) != nfu_Okay ) goto Err;
                        if( s2 == 0. ) {
                            y2 = nan; }
                        else {
                            y2 = s1 / s2;
                        }
                        iYs++;
                    }
                    p->y = ( y1 + y2 ) / iYs; 
                    if( nfu_isNAN( p->y ) ) zeros++; }
                else {
                    if( !safeDivide ) {
                        *status = nfu_divByZero;
                        goto Err;
                    }
                    zeros++;
                    p->y = nan;
                } }
            else {
                p->y /= y;
            }
        }
        length = n->length - 1;
        if( length > 0 ) {
            x2 = n->points[length].x;
            for( i = length - 1; i >= 0; i-- ) {             /* Find and add y zeros and NAN not currently in n's. */
                x1 = n->points[i].x;
                if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x1, &u1 ) ) != nfu_Okay ) goto Err;
                if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x2, &u2 ) ) != nfu_Okay ) goto Err;
                if( ( *status = ptwXY_getValueAtX( ptwXY2, x1, &v1 ) ) != nfu_Okay ) goto Err;
                if( ( *status = ptwXY_getValueAtX( ptwXY2, x2, &v2 ) ) != nfu_Okay ) goto Err;
                if( u1 * u2 < 0 ) {
                    xz = ( u1 * x2 - u2 * x1 ) / ( u1 - u2 );
                    if( ( *status = ptwXY_setValueAtX( n, xz, 0. ) ) != nfu_Okay ) goto Err;
                }
                if( v1 * v2 < 0 ) {
                    if( !safeDivide ) {
                        *status = nfu_divByZero;
                        goto Err;
                    }
                    zeros++;
                    xz = ( v1 * x2 - v2 * x1 ) / ( v1 - v2 );
                    if( ( *status = ptwXY_setValueAtX( n, xz, nan ) ) != nfu_Okay ) goto Err;
                }
                x2 = x1;
            }
            if( ( *status = ptwXY_simpleCoalescePoints( n ) ) != nfu_Okay ) goto Err;
            length = n->length;
            x2 = n->points[n->length-1].x;
            y2 = n->points[n->length-1].y;
            isNAN2 = nfu_isNAN( y2 );
            for( i = n->length - 2; i >= 0; i-- ) {             /* Make interpolation fit accuracy. Work backwards so new points will not mess up loop. */
                x1 = n->points[i].x;
                y1 = n->points[i].y;
                isNAN1 = nfu_isNAN( y1 );
                if( !isNAN1 || !isNAN2 ) {
                    if( ( *status = ptwXY_div_s_ptwXY( n, ptwXY1, ptwXY2, x1, y1, x2, y2, 0, isNAN1, isNAN2 ) ) != nfu_Okay ) goto Err;
                }
                x2 = x1;
                y2 = y1;
                isNAN2 = isNAN1;
            }
            ptwXY_update_biSectionMax( n, (double) length );
            if( zeros ) {
                if( ( *status = ptwXY_simpleCoalescePoints( n ) ) != nfu_Okay ) goto Err;
                for( i = 0; i < n->length; i++ ) if( !nfu_isNAN( n->points[i].y ) ) break;
                if( nfu_isNAN( n->points[0].y ) ) {                     /* Special case for first point. */
                    if( i == n->length ) {                              /* They are all nan's, what now? */
                        zeros = 0;
                        for( i = 0; i < n->length; i++ ) n->points[i].y = 0.; }
                    else {
                         n->points[0].y = 2. * n->points[i].y;
                        zeros--;
                    }
                }
                for( i = n->length - 1; i > 0; i-- ) if( !nfu_isNAN( n->points[i].y ) ) break;
                if( nfu_isNAN( n->points[n->length - 1].y ) ) {         /* Special case for last point. */
                    n->points[n->length - 1].y = 2. * n->points[i].y;
                    zeros--;
                }
                if( zeros ) {
                    for( i = 0; i < n->length; i++ ) if( nfu_isNAN( n->points[i].y ) ) break;
                    for( k = i + 1, j = i; k < n->length; k++ ) {
                    if( nfu_isNAN( n->points[k].y ) ) continue;
                        n->points[j] = n->points[k];
                        j++;
                    }
                    n->length = j;
                }
            }
        }
    }
    return( n );
 
Err:
    if( n ) ptwXY_free( n );
    return( NULL );
}

ptwXY_dullEdges()

nfu_status ptwXY_dullEdges	(	ptwXYPoints *	ptwXY,
		double	lowerEps,
		double	upperEps,
		int	positiveXOnly )

Definition at line 42 of file ptwXY_convenient.cc.

                                                                                                      {
 
#define minEps 5e-16
 
    nfu_status status;
    double xm, xp, dx, y, x1, y1, x2, y2, sign;
    ptwXYPoint *p;
 
/* This routine can only be used for linear interpolation for the y-axes since for log interpolation, y cannot be 0. 
This needs to be fixed and documented. */
    if( ( status = ptwXY->status ) != nfu_Okay ) return( status );
    if( ptwXY->interpolation == ptwXY_interpolationFlat ) return( nfu_invalidInterpolation );
    if( ptwXY->interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
 
    if( ptwXY->length < 2 ) return( nfu_Okay );
 
    if( lowerEps != 0. ) {
        if( std::fabs( lowerEps ) < minEps ) {
            sign = 1;
            if( lowerEps < 0. ) sign = -1;
            lowerEps = sign * minEps;
        }
 
        p = ptwXY_getPointAtIndex_Unsafely( ptwXY, 0 );
        x1 = p->x;
        y1 = p->y;
        p = ptwXY_getPointAtIndex_Unsafely( ptwXY, 1 );
        x2 = p->x;
        y2 = p->y;
 
        if( y1 != 0. ) {
            dx = std::fabs( x1 * lowerEps );
            if( x1 == 0 ) dx = std::fabs( lowerEps );
            xm = x1 - dx;
            xp = x1 + dx;
            if( ( xp + dx ) < x2 ) {
                if( ( status = ptwXY_getValueAtX( ptwXY, xp, &y ) ) != nfu_Okay ) return( status );
                if( ( status = ptwXY_setValueAtX( ptwXY, xp,  y ) ) != nfu_Okay ) return( status ); }
            else {
                xp = x2;
                y = y2;
            }
            if( lowerEps > 0 ) {
                if( ( status = ptwXY_setValueAtX( ptwXY, x1, 0. ) ) != nfu_Okay ) return( status ); }
            else {
                if( ( xm < 0. ) && ( x1 >= 0. ) && positiveXOnly ) {
                    if( ( status = ptwXY_setValueAtX( ptwXY, x1, 0. ) ) != nfu_Okay ) return( status ); }
                else {
                    if( ( status = ptwXY_setValueAtX( ptwXY, xm, 0. ) ) != nfu_Okay ) return( status );
                    if( ( status = ptwXY_interpolatePoint( ptwXY->interpolation, x1, &y, xm, 0., xp, y )  ) != nfu_Okay ) return( status );
                    if( ( status = ptwXY_setValueAtX( ptwXY, x1, y ) ) != nfu_Okay ) return( status );
                }
            }
        }
    }
 
    if( upperEps != 0. ) {
        if( std::fabs( upperEps ) < minEps ) {
            sign = 1;
            if( upperEps < 0. ) sign = -1;
            upperEps = sign * minEps;
        }
 
        p = ptwXY_getPointAtIndex_Unsafely( ptwXY, ptwXY->length - 2 );
        x1 = p->x;
        y1 = p->y;
        p = ptwXY_getPointAtIndex_Unsafely( ptwXY, ptwXY->length - 1 );
        x2 = p->x;
        y2 = p->y;
 
        if( y2 != 0. ) {
            dx = std::fabs( x2 * upperEps );
            if( x2 == 0 ) dx = std::fabs( upperEps );
            xm = x2 - dx;
            xp = x2 + dx;
            if( ( xm - dx ) > x1 ) {
                if( ( status = ptwXY_getValueAtX( ptwXY, xm, &y ) ) != nfu_Okay ) return( status );
                if( ( status = ptwXY_setValueAtX( ptwXY, xm,  y ) ) != nfu_Okay ) return( status ); }
            else {
                xm = x1;
                y = y1;
            }
            if( upperEps < 0 ) {
                if( ( status = ptwXY_setValueAtX( ptwXY, x2, 0. ) ) != nfu_Okay ) return( status ); }
            else {
                if( ( status = ptwXY_setValueAtX( ptwXY, xp, 0. ) ) != nfu_Okay ) return( status );
                if( ( status = ptwXY_interpolatePoint( ptwXY->interpolation, x2, &y, xm, y, xp, 0. )  ) != nfu_Okay ) return( status );
                if( ( status = ptwXY_setValueAtX( ptwXY, x2, y ) ) != nfu_Okay ) return( status );
            }
        }
    }
 
    return( ptwXY->status );
 
#undef minEps
}

Referenced by MCGIDI_reaction_fixDomains(), and ptwXY_mutualifyDomains().

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_f_integrate()

nfu_status ptwXY_f_integrate	(	ptwXY_interpolation	interpolation,
		double	x1,
		double	y1,
		double	x2,
		double	y2,
		double *	value )

Definition at line 34 of file ptwXY_integration.cc.

                                                                                                                             {
 
    nfu_status status = nfu_Okay;
    double r;
 
    *value = 0.;
    switch( interpolation ) {
    case ptwXY_interpolationLinLin :                            /* x linear, y linear */
        *value = 0.5 * ( y1 + y2 ) * ( x2 - x1 );
        break;
    case ptwXY_interpolationLinLog :                            /* x linear, y log */
        if( ( y1 <= 0. ) || ( y2 <= 0. ) ) {
            status = nfu_badIntegrationInput; }
        else {
            r = y2 / y1;
            if( std::fabs( r - 1. ) < 1e-4 ) {
                r = r - 1.;
                *value = y1 * ( x2 - x1 ) / ( 1. + r * ( -0.5 + r * ( 1. / 3. + r * ( -0.25 + .2 * r ) ) ) ); }
            else {
                *value = ( y2 - y1 ) * ( x2 - x1 ) / G4Log( r );
            }
        }
        break;
    case ptwXY_interpolationLogLin :                            /* x log, y linear */
        if( ( x1 <= 0. ) || ( x2 <= 0. ) ) {
            status = nfu_badIntegrationInput; }
        else {
            r = x2 / x1;
            if( std::fabs( r - 1. ) < 1e-4 ) {
                r = r - 1.;
                r = r * ( -0.5 + r * ( 1. / 3. + r * ( -0.25 + .2 * r ) ) );
                *value = x1 * ( y2 - y1 ) * r / ( 1. + r ) + y2 * ( x2 - x1 ); }
            else {
                *value = ( y1 - y2 ) * ( x2 - x1 ) / G4Log( r ) + x2 * y2 - x1 * y1;
            }
        }
        break;
    case ptwXY_interpolationLogLog :                            /* x log, y log */
        if( ( x1 <= 0. ) || ( x2 <= 0. ) || ( y1 <= 0. ) || ( y2 <= 0. ) ) {
            status = nfu_badIntegrationInput; }
        else {
            int i, n;
            double a, z, lx, ly, s, f;
 
            r = y2 / y1;
            if( std::fabs( r - 1. ) < 1e-4 ) {
                ly = ( y2 - y1 ) / y1;
                ly = ly * ( 1. + ly * ( -0.5 + ly * ( 1. / 3. - 0.25 * ly ) ) ); }
            else {
                ly = G4Log( r );
            }
            r = x2 / x1;
            if( std::fabs( r - 1. ) < 1e-4 ) {
                lx = ( x2 - x1 ) / x1;
                lx = lx * ( 1 + lx * ( -0.5 + lx * ( 1. / 3. - 0.25 * lx ) ) ); }
            else {
                lx = G4Log( r );
            }
            a = ly / lx;
            if( std::fabs( r - 1. ) < 1e-3 ) {
                z = ( x2 - x1 ) / x1;
                n = (int) a;
                if( n > 10 ) n = 12;
                if( n < 4 ) n = 6;
                a = a - n + 1;
                f = n + 1.;
                for( i = 0, s = 0.; i < n; i++, a++, f-- ) s = ( 1. + s ) * a * z / f;
                *value = y1 * ( x2 - x1 ) * ( 1. + s ); }
            else {
                *value = y1 * x1 * ( G4Pow::GetInstance()->powA( r, a + 1. ) - 1. ) / ( a + 1. );
            }
        }
        break;
    case ptwXY_interpolationFlat :                            /* x ?, y flat */
        *value = y1 * ( x2 - x1 );
        break;
    case ptwXY_interpolationOther :
        status = nfu_otherInterpolation;
    }
    return( status );
}

Referenced by ptwXY_integrate(), and ptwXY_runningIntegral().

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_free()

ptwXYPoints * ptwXY_free

(

ptwXYPoints *

ptwXY

)

Definition at line 574 of file ptwXY_core.cc.

                                              {
 
    if( ptwXY != NULL ) ptwXY_release( ptwXY );
    nfu_free( ptwXY );
    return( (ptwXYPoints *) NULL );
}

Referenced by GIDI_settings_processedFlux::GIDI_settings_processedFlux(), GIDI_settings_processedFlux::groupFunction(), MCGIDI_angular_parseFromTOM(), MCGIDI_energy_release(), MCGIDI_fromTOM_pdfsOfXGivenW(), MCGIDI_product_parseFromTOM(), MCGIDI_product_release(), MCGIDI_reaction_release(), MCGIDI_target_heated_read(), MCGIDI_target_heated_release(), GIDI_settings_processedFlux::operator=(), ptwXY_binary_ptwXY(), ptwXY_clip(), ptwXY_convolution(), ptwXY_create(), ptwXY_createFromFunction(), ptwXY_createGaussian(), ptwXY_createGaussianCenteredSigma1(), ptwXY_div_ptwXY(), ptwXY_flatInterpolationToLinear(), ptwXY_groupOneFunction(), ptwXY_groupThreeFunctions(), ptwXY_groupTwoFunctions(), ptwXY_intersectionWith_ptwX(), ptwXY_mul2_ptwXY(), ptwXY_sub_ptwXY(), ptwXY_thin(), ptwXY_toOtherInterpolation(), ptwXY_union(), ptwXY_unitbaseInterpolate(), ptwXY_xSlice(), and GIDI_settings_processedFlux::~GIDI_settings_processedFlux().

ptwXY_fromString()

ptwXYPoints * ptwXY_fromString	(	char const *	str,
		ptwXY_interpolation	interpolation,
		ptwXY_interpolationOtherInfo const *	interpolationOtherInfo,
		double	biSectionMax,
		double	accuracy,
		char **	endCharacter,
		nfu_status *	status )

Definition at line 236 of file ptwXY_misc.cc.

                                                                                    {
 
    int64_t numberConverted;
    double  *doublePtr;
    ptwXYPoints *ptwXY = NULL;
 
    if( ( *status = nfu_stringToListOfDoubles( str, &numberConverted, &doublePtr, endCharacter ) ) != nfu_Okay ) return( NULL );
    *status = nfu_oddNumberOfValues;
    if( ( numberConverted % 2 ) == 0 )
        ptwXY = ptwXY_create( interpolation, interpolationOtherInfo, biSectionMax, accuracy, numberConverted, 10, numberConverted / 2, doublePtr, status, 0 );
    nfu_free( doublePtr );
    return( ptwXY );
}

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_getAccuracy()

double ptwXY_getAccuracy

(

ptwXYPoints *

ptwXY

)

Definition at line 372 of file ptwXY_core.cc.

                                               {
 
    return( ptwXY->accuracy );
}

ptwXY_getBiSectionMax()

double ptwXY_getBiSectionMax

(

ptwXYPoints *

ptwXY

)

Definition at line 390 of file ptwXY_core.cc.

                                                   {
 
    return( ptwXY->biSectionMax );
}

ptwXY_getInterpolation()

ptwXY_interpolation ptwXY_getInterpolation

(

ptwXYPoints *

ptwXY

)

Definition at line 337 of file ptwXY_core.cc.

                                                                 {
 
    return( ptwXY->interpolation );
}

ptwXY_getInterpolationString()

char const * ptwXY_getInterpolationString

(

ptwXYPoints *

ptwXY

)

Definition at line 344 of file ptwXY_core.cc.

                                                               {
 
    return( ptwXY->interpolationOtherInfo.interpolationString );
}

ptwXY_getNonOverflowLength()

int64_t ptwXY_getNonOverflowLength

(

ptwXYPoints const *

ptwXY

)

Definition at line 590 of file ptwXY_core.cc.

                                                               {
 
    return( ptwXY->length - ptwXY->overflowLength );
}

Referenced by ptwXY_abs(), ptwXY_appendXY(), ptwXY_coalescePoints(), ptwXY_copy(), ptwXY_div_fromDouble(), ptwXY_getPointsAroundX_closeIsEqual(), ptwXY_getXMaxAndFrom(), ptwXY_getXMinAndFrom(), ptwXY_getYMax(), ptwXY_getYMin(), ptwXY_mod(), ptwXY_neg(), ptwXY_setValueAtX_overrideIfClose(), ptwXY_showInteralStructure(), and ptwXY_slopeOffset().

ptwXY_getPointAtIndex()

ptwXYPoint * ptwXY_getPointAtIndex	(	ptwXYPoints *	ptwXY,
		int64_t	index )

Definition at line 675 of file ptwXY_core.cc.

                                                                       {
 
    if( ptwXY->status != nfu_Okay ) return( NULL );
    if( ( index < 0 ) || ( index >= ptwXY->length ) ) return( NULL );
    return( ptwXY_getPointAtIndex_Unsafely( ptwXY, index ) );
}

Referenced by ptwXY_getXYPairAtIndex(), ptwXY_showInteralStructure(), and ptwXY_simpleWrite().

ptwXY_getPointAtIndex_Unsafely()

ptwXYPoint * ptwXY_getPointAtIndex_Unsafely	(	ptwXYPoints *	ptwXY,
		int64_t	index )

Definition at line 684 of file ptwXY_core.cc.

                                                                                {
 
    int64_t i;
    ptwXYOverflowPoint *overflowPoint;
 
    for( overflowPoint = ptwXY->overflowHeader.next, i = 0; overflowPoint != &(ptwXY->overflowHeader); overflowPoint = overflowPoint->next, i++ ) {
        if( overflowPoint->index == index ) return( &(overflowPoint->point) );
        if( overflowPoint->index > index ) break;
    }
    return( &(ptwXY->points[index - i]) );
}

Referenced by MCGIDI_angular_parseFromTOM(), MCGIDI_fromTOM_pdfOfX(), ptwXY_areDomainsMutual(), ptwXY_clip(), ptwXY_dullEdges(), ptwXY_getPointAtIndex(), ptwXY_getSlopeAtX(), ptwXY_mutualifyDomains(), and ptwXY_tweakDomainsToMutualify().

ptwXY_getPointsAroundX()

ptwXY_lessEqualGreaterX ptwXY_getPointsAroundX	(	ptwXYPoints *	ptwXY,
		double	x,
		ptwXYOverflowPoint *	lessThanEqualXPoint,
		ptwXYOverflowPoint *	greaterThanXPoint )

Definition at line 710 of file ptwXY_core.cc.

                                                                                                                                                               {
 
    int closeIsEqual;
    ptwXYPoint *closePoint;
 
    return( ptwXY_getPointsAroundX_closeIsEqual( ptwXY, x, lessThanEqualXPoint, greaterThanXPoint, 0, &closeIsEqual, &closePoint ) );
}

Referenced by ptwXY_getSlopeAtX(), and ptwXY_getValueAtX().

ptwXY_getPointsAroundX_closeIsEqual()

ptwXY_lessEqualGreaterX ptwXY_getPointsAroundX_closeIsEqual	(	ptwXYPoints *	ptwXY,
		double	x,
		ptwXYOverflowPoint *	lessThanEqualXPoint,
		ptwXYOverflowPoint *	greaterThanXPoint,
		double	eps,
		int *	closeIsEqual,
		ptwXYPoint **	closePoint )

Definition at line 720 of file ptwXY_core.cc.

                                                                                                        {
 
    int64_t overflowIndex, nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
    int64_t indexMin, indexMid, indexMax;
    ptwXY_dataFrom xMinFrom, xMaxFrom;
    double xMin = ptwXY_getXMinAndFrom( ptwXY, &xMinFrom ), xMax = ptwXY_getXMaxAndFrom( ptwXY, &xMaxFrom );
    ptwXYOverflowPoint *overflowPoint, *overflowHeader = &(ptwXY->overflowHeader);
    ptwXY_lessEqualGreaterX status = ptwXY_lessEqualGreaterX_empty;
    ptwXYPoint *lowerPoint = NULL, *upperPoint = NULL;
 
    ptwXY_initialOverflowPoint( lessThanEqualXPoint, overflowHeader, NULL );
    ptwXY_initialOverflowPoint( greaterThanXPoint, overflowHeader, NULL );
    if( ptwXY->length != 0 ) {
        if( x < xMin ) {
            status = ptwXY_lessEqualGreaterX_lessThan;
            if( xMinFrom == ptwXY_dataFrom_Points ) {
                greaterThanXPoint->prior = overflowHeader;
                greaterThanXPoint->index = 0;
                greaterThanXPoint->point = ptwXY->points[0];
                *closePoint = &(ptwXY->points[0]); }
            else {
                *greaterThanXPoint = *(overflowHeader->next);
                *closePoint = &(overflowHeader->next->point);
            } }
        else if( x > xMax ) {
            status = ptwXY_lessEqualGreaterX_greater;
            if( xMaxFrom == ptwXY_dataFrom_Points ) {
                lessThanEqualXPoint->prior = overflowHeader->prior;
                lessThanEqualXPoint->index = nonOverflowLength - 1;
                lessThanEqualXPoint->point = ptwXY->points[lessThanEqualXPoint->index];
                *closePoint = &(ptwXY->points[lessThanEqualXPoint->index]); }
            else {
                *lessThanEqualXPoint = *(overflowHeader->prior);
                *closePoint = &(overflowHeader->prior->point);
            } }
        else {                                                  /* xMin <= x <= xMax */
            status = ptwXY_lessEqualGreaterX_between;           /* Default for this condition, can only be between or equal. */
            for( overflowPoint = overflowHeader->next, overflowIndex = 0; overflowPoint != overflowHeader; 
                overflowPoint = overflowPoint->next, overflowIndex++ ) if( overflowPoint->point.x > x ) break;
            overflowPoint = overflowPoint->prior;
            if( ( overflowPoint != overflowHeader ) && ( overflowPoint->point.x == x ) ) {
                status = ptwXY_lessEqualGreaterX_equal;
                *lessThanEqualXPoint = *overflowPoint; }
            else if( ptwXY->length == 1 ) {                    /* If here and length = 1, then ptwXY->points[0].x == x. */
                status = ptwXY_lessEqualGreaterX_equal;
                lessThanEqualXPoint->index = 0;
                lessThanEqualXPoint->point = ptwXY->points[0]; }
            else {                                              /* ptwXY->length > 1 */
                indexMin = 0;
                indexMax = nonOverflowLength - 1;
                indexMid = ( indexMin + indexMax ) >> 1;
                while( ( indexMin != indexMid ) && ( indexMid != indexMax ) ) {
                    if( ptwXY->points[indexMid].x > x ) {
                        indexMax = indexMid; }
                    else {
                        indexMin = indexMid;
                    }
                    indexMid = ( indexMin + indexMax ) >> 1;
                } // Loop checking, 11.06.2015, T. Koi
                if( ptwXY->points[indexMin].x == x ) {
                    status = ptwXY_lessEqualGreaterX_equal;
                    lessThanEqualXPoint->index = indexMin;
                    lessThanEqualXPoint->point = ptwXY->points[indexMin]; }
                else if( ptwXY->points[indexMax].x == x ) {
                    status = ptwXY_lessEqualGreaterX_equal;
                    lessThanEqualXPoint->index = indexMax;
                    lessThanEqualXPoint->point = ptwXY->points[indexMax]; }
                else {
                    if( ptwXY->points[indexMin].x > x ) indexMax = 0;
                    if( ptwXY->points[indexMax].x < x ) indexMin = indexMax;
                    if( ( overflowPoint == overflowHeader ) ||     /* x < xMin of overflow points. */
                            ( ( ptwXY->points[indexMin].x > overflowPoint->point.x ) && ( ptwXY->points[indexMin].x < x ) ) ) {
                        if( overflowPoint != overflowHeader ) lessThanEqualXPoint->prior = overflowPoint;
                        lowerPoint = &(ptwXY->points[indexMin]);
                        lessThanEqualXPoint->index = indexMin;
                        lessThanEqualXPoint->point = ptwXY->points[indexMin]; }
                    else {
                        lowerPoint = &(overflowPoint->point);
                        *lessThanEqualXPoint = *overflowPoint;
                    }
                    if( ( overflowPoint->next == overflowHeader ) ||     /* x > xMax of overflow points. */
                            ( ( ptwXY->points[indexMax].x < overflowPoint->next->point.x ) && ( ptwXY->points[indexMax].x > x ) ) ) {
                        upperPoint = &(ptwXY->points[indexMax]);
                        greaterThanXPoint->index = indexMax;
                        greaterThanXPoint->point = ptwXY->points[indexMax]; }
                    else {
                        upperPoint = &(overflowPoint->next->point);
                        *greaterThanXPoint = *(overflowPoint->next);
                    }
                }
            }
        }
    }
 
    *closeIsEqual = 0;
    if( eps > 0 ) {
        double absX = std::fabs( x );
 
        if( status == ptwXY_lessEqualGreaterX_lessThan ) {
            if( absX < std::fabs( greaterThanXPoint->point.x ) ) absX = std::fabs( greaterThanXPoint->point.x );
            if( ( greaterThanXPoint->point.x - x ) < eps * absX ) *closeIsEqual = 1; }
        else if( status == ptwXY_lessEqualGreaterX_greater ) {
            if( absX < std::fabs( lessThanEqualXPoint->point.x ) ) absX = std::fabs( lessThanEqualXPoint->point.x );
            if( ( x - lessThanEqualXPoint->point.x ) < eps * absX ) *closeIsEqual = -1; }
        else if( status == ptwXY_lessEqualGreaterX_between ) {
            if( ( x - lessThanEqualXPoint->point.x ) < ( greaterThanXPoint->point.x - x ) ) {   /* x is closer to lower point. */
                *closePoint = lowerPoint;
                if( absX < std::fabs( lessThanEqualXPoint->point.x ) ) absX = std::fabs( lessThanEqualXPoint->point.x );
                if( ( x - lessThanEqualXPoint->point.x ) < eps * absX ) *closeIsEqual = -1; }
            else {                                                                              /* x is closer to upper point. */
                *closePoint = upperPoint;
                if( absX < std::fabs( greaterThanXPoint->point.x ) ) absX = std::fabs( greaterThanXPoint->point.x );
                if( ( greaterThanXPoint->point.x - x ) < eps * absX ) *closeIsEqual = 1;
            } }
        else if( status == ptwXY_lessEqualGreaterX_equal ) {
            *closeIsEqual = 1;
        }
    }
    return( status );
}

Referenced by ptwXY_getPointsAroundX(), and ptwXY_setValueAtX_overrideIfClose().

ptwXY_getSlopeAtX()

nfu_status ptwXY_getSlopeAtX	(	ptwXYPoints *	ptwXY,
		double	x,
		const char	side,
		double *	slope )

Definition at line 1139 of file ptwXY_core.cc.

                                                                                             {
 
    nfu_status status  = nfu_Okay;
    ptwXYOverflowPoint lessThanEqualXPoint, greaterThanXPoint;
    ptwXY_lessEqualGreaterX legx = ptwXY_getPointsAroundX( ptwXY, x, &lessThanEqualXPoint, &greaterThanXPoint );
    ptwXYPoint *point;
 
    *slope = 0.;
    if( ( side != '-' ) && ( side != '+' ) ) return( nfu_badInput );
 
    switch( legx ) {
    case ptwXY_lessEqualGreaterX_empty :
    case ptwXY_lessEqualGreaterX_lessThan :
    case ptwXY_lessEqualGreaterX_greater :
        status = nfu_XOutsideDomain;
        break;
    case ptwXY_lessEqualGreaterX_between :
        *slope = ( greaterThanXPoint.point.y - lessThanEqualXPoint.point.y ) / 
            ( greaterThanXPoint.point.x - lessThanEqualXPoint.point.x );
        break;
    case ptwXY_lessEqualGreaterX_equal :
        if( side == '-' ) {
            if( lessThanEqualXPoint.index == 0 ) {
                status = nfu_XOutsideDomain; }
            else {
                point = ptwXY_getPointAtIndex_Unsafely( ptwXY, lessThanEqualXPoint.index - 1 );
                *slope = ( lessThanEqualXPoint.point.y - point->y ) / ( lessThanEqualXPoint.point.x - point->x );
            } }
        else {
            if( lessThanEqualXPoint.index == ( ptwXY->length - 1 ) ) {
                status = nfu_XOutsideDomain; }
            else {
                point = ptwXY_getPointAtIndex_Unsafely( ptwXY, lessThanEqualXPoint.index + 1 );
                *slope = ( point->y - lessThanEqualXPoint.point.y ) / ( point->x - lessThanEqualXPoint.point.x );
            }
        }
    }
 
    return( status );
}

Referenced by ptwXY_div_ptwXY().

ptwXY_getStatus()

nfu_status ptwXY_getStatus

(

ptwXYPoints *

ptwXY

)

Definition at line 351 of file ptwXY_core.cc.

                                                 {
 
    return( ptwXY->status );
}

Referenced by nf_Legendre_from_ptwXY().

ptwXY_getUserFlag()

int ptwXY_getUserFlag

(

ptwXYPoints *

ptwXY

)

Definition at line 358 of file ptwXY_core.cc.

                                            {
 
    return( ptwXY->userFlag );
}

ptwXY_getValueAtX()

nfu_status ptwXY_getValueAtX	(	ptwXYPoints *	ptwXY,
		double	x,
		double *	y )

Definition at line 844 of file ptwXY_core.cc.

                                                                        {
 
    nfu_status status = nfu_XOutsideDomain;
    ptwXYOverflowPoint lessThanEqualXPoint, greaterThanXPoint;
    ptwXY_lessEqualGreaterX legx = ptwXY_getPointsAroundX( ptwXY, x, &lessThanEqualXPoint, &greaterThanXPoint );
 
    *y = 0.;
    if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
    switch( legx ) {
    case ptwXY_lessEqualGreaterX_empty :
    case ptwXY_lessEqualGreaterX_lessThan :
    case ptwXY_lessEqualGreaterX_greater :
        break;
    case ptwXY_lessEqualGreaterX_equal :
        status = nfu_Okay;
        *y = lessThanEqualXPoint.point.y;
        break;
    case ptwXY_lessEqualGreaterX_between :
        if( ptwXY->interpolationOtherInfo.getValueFunc != NULL ) {
            status = ptwXY->interpolationOtherInfo.getValueFunc( ptwXY->interpolationOtherInfo.argList, x, y, 
                lessThanEqualXPoint.point.x, lessThanEqualXPoint.point.y, greaterThanXPoint.point.x, greaterThanXPoint.point.y ); }
        else {
            status = ptwXY_interpolatePoint( ptwXY->interpolation, x, y, lessThanEqualXPoint.point.x, lessThanEqualXPoint.point.y, 
                greaterThanXPoint.point.x, greaterThanXPoint.point.y );
        }
        break;
    }
    return( status );
}

Referenced by MCGIDI_reaction_getCrossSectionAtE(), MCGIDI_sampling_ptwXY_getValueAtX(), MCGIDI_target_heated_getTotalCrossSectionAtE(), ptwXY_div_ptwXY(), ptwXY_dullEdges(), ptwXY_intersectionWith_ptwX(), and ptwXY_xSlice().

ptwXY_getXArray()

ptwXPoints * ptwXY_getXArray	(	ptwXYPoints *	ptwXY,
		nfu_status *	status )

Definition at line 24 of file ptwXY_convenient.cc.

                                                                      {
 
    int64_t i, n;
    ptwXPoints *xArray;
 
    if( ( *status = ptwXY->status ) != nfu_Okay ) return( NULL );
    n = ptwXY->length;
 
    if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( NULL );
    if( ( xArray = ptwX_new( n, status ) ) == NULL ) return( NULL );
    for( i = 0; i < n; i++ ) xArray->points[i] = ptwXY->points[i].x;
    xArray->length = n;
 
    return( xArray );
}

ptwXY_getXMax()

double ptwXY_getXMax

(

ptwXYPoints *

ptwXY

)

Definition at line 1239 of file ptwXY_core.cc.

                                           {
 
    ptwXY_dataFrom dataFrom;
 
    return( ptwXY_getXMaxAndFrom( ptwXY, &dataFrom ) );
}

Referenced by GIDI_settings_processedFlux::groupFunction(), MCGIDI_product_sampleMultiplicity(), MCGIDI_reaction_parseFromTOM(), MCGIDI_sampling_ptwXY_getValueAtX(), ptwXY_integrateDomain(), ptwXY_integrateDomainWithWeight_sqrt_x(), ptwXY_integrateDomainWithWeight_x(), ptwXY_xMinSlice(), and ptwXY_xSlice().

ptwXY_getXMaxAndFrom()

double ptwXY_getXMaxAndFrom	(	ptwXYPoints *	ptwXY,
		ptwXY_dataFrom *	dataFrom )

Definition at line 1215 of file ptwXY_core.cc.

                                                                            {
 
    int64_t nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
    double xMax = nfu_getNAN( );
 
    *dataFrom = ptwXY_dataFrom_Unknown;
    if( ptwXY->overflowLength > 0 ) {
        *dataFrom = ptwXY_dataFrom_Overflow;
        xMax = ptwXY->overflowHeader.prior->point.x;
        if( ( nonOverflowLength > 0 ) ) {
            if( xMax < ptwXY->points[nonOverflowLength-1].x ) {
                *dataFrom = ptwXY_dataFrom_Points;
                xMax = ptwXY->points[nonOverflowLength-1].x;
            }
        } }
    else if( ptwXY->length > 0 ) {
        *dataFrom = ptwXY_dataFrom_Points;
        xMax = ptwXY->points[nonOverflowLength-1].x;
    }
    return( xMax );
}

Referenced by ptwXY_appendXY(), ptwXY_getPointsAroundX_closeIsEqual(), and ptwXY_getXMax().

ptwXY_getXMin()

double ptwXY_getXMin

(

ptwXYPoints *

ptwXY

)

Definition at line 1206 of file ptwXY_core.cc.

                                           {
 
    ptwXY_dataFrom dataFrom;
 
    return( ptwXY_getXMinAndFrom( ptwXY, &dataFrom ) );
}

Referenced by GIDI_settings_processedFlux::groupFunction(), MCGIDI_reaction_parseFromTOM(), MCGIDI_sampling_ptwXY_getValueAtX(), ptwXY_integrateDomain(), ptwXY_integrateDomainWithWeight_sqrt_x(), ptwXY_integrateDomainWithWeight_x(), ptwXY_xMaxSlice(), and ptwXY_xSlice().

ptwXY_getXMinAndFrom()

double ptwXY_getXMinAndFrom	(	ptwXYPoints *	ptwXY,
		ptwXY_dataFrom *	dataFrom )

Definition at line 1182 of file ptwXY_core.cc.

                                                                            {
 
    int64_t nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
    double xMin = nfu_getNAN( );
 
    *dataFrom = ptwXY_dataFrom_Unknown;
    if( ptwXY->overflowLength > 0 ) {
        *dataFrom = ptwXY_dataFrom_Overflow;
        xMin = ptwXY->overflowHeader.next->point.x;
        if( nonOverflowLength >= 0 ) {
            if( xMin > ptwXY->points[0].x ) {
                *dataFrom = ptwXY_dataFrom_Points;
                xMin = ptwXY->points[0].x;
            }
        } }
    else if( nonOverflowLength > 0 ) {
        *dataFrom = ptwXY_dataFrom_Points;
        xMin = ptwXY->points[0].x;
    }
    return( xMin );
}

Referenced by ptwXY_getPointsAroundX_closeIsEqual(), and ptwXY_getXMin().

ptwXY_getXYPairAtIndex()

nfu_status ptwXY_getXYPairAtIndex	(	ptwXYPoints *	ptwXY,
		int64_t	index,
		double *	x,
		double *	y )

Definition at line 698 of file ptwXY_core.cc.

                                                                                             {
 
    ptwXYPoint *p = ptwXY_getPointAtIndex( ptwXY, index );
 
    if( p == NULL ) return( nfu_badIndex );
    *x = p->x;
    *y = p->y;
    return( nfu_Okay );
}

Referenced by nf_Legendre_from_ptwXY().

ptwXY_getYMax()

double ptwXY_getYMax

(

ptwXYPoints *

ptwXY

)

Definition at line 1269 of file ptwXY_core.cc.

                                           {
 
    int64_t i, n = ptwXY_getNonOverflowLength( ptwXY  );
    ptwXYPoint *p = ptwXY->points;
    ptwXYOverflowPoint *overflowPoint = ptwXY->overflowHeader.next;
    double yMax;
 
    if( ptwXY->length == 0 ) return( 0. );
    if( n > 0 ) {
        yMax = p->y;
        for( i = 1, p++; i < n; i++, p++ ) yMax = ( ( yMax > p->y ) ? yMax : p->y ); }
    else {
        yMax = overflowPoint->point.y;
    }
    for( ; overflowPoint != &(ptwXY->overflowHeader); overflowPoint = overflowPoint->next )
        yMax = ( ( yMax > overflowPoint->point.y ) ? yMax : overflowPoint->point.y );
    return( yMax );
}

Referenced by MCGIDI_product_parseFromTOM(), and ptwXY_clip().

ptwXY_getYMin()

double ptwXY_getYMin

(

ptwXYPoints *

ptwXY

)

Definition at line 1248 of file ptwXY_core.cc.

                                           {
 
    int64_t i, n = ptwXY_getNonOverflowLength( ptwXY  );
    ptwXYPoint *p = ptwXY->points;
    ptwXYOverflowPoint *overflowPoint = ptwXY->overflowHeader.next;
    double yMin;
 
    if( ptwXY->length == 0 ) return( 0. );
    if( n > 0 ) {
        yMin = p->y;
        for( i = 1, p++; i < n; i++, p++ ) yMin = ( ( yMin < p->y ) ? yMin : p->y ); }
    else {
        yMin = overflowPoint->point.y;
    }
    for( ; overflowPoint != &(ptwXY->overflowHeader); overflowPoint = overflowPoint->next ) 
        yMin = ( ( yMin < overflowPoint->point.y ) ? yMin : overflowPoint->point.y );
    return( yMin );
}

Referenced by MCGIDI_product_parseFromTOM(), and ptwXY_clip().

ptwXY_groupOneFunction()

ptwXPoints * ptwXY_groupOneFunction	(	ptwXYPoints *	ptwXY,
		ptwXPoints *	groupBoundaries,
		ptwXY_group_normType	normType,
		ptwXPoints *	ptwX_norm,
		nfu_status *	status )

Definition at line 363 of file ptwXY_integration.cc.

                                                                                                                                                                {
 
    int64_t i, igs, ngs;
    double x1, y1, x2, y2, y2p, xg1, xg2, sum;
    ptwXYPoints *f;
    ptwXPoints *groupedData = NULL;
 
    if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( NULL );
    if( ( *status = groupBoundaries->status ) != nfu_Okay ) return( NULL );
    *status = nfu_otherInterpolation;
    if( ptwXY->interpolation == ptwXY_interpolationOther ) return( NULL );
 
    ngs = ptwX_length( groupBoundaries ) - 1;
    if( normType == ptwXY_group_normType_norm ) {
        if( ptwX_norm == NULL ) {
            *status = nfu_badNorm;
            return( NULL );
        }
        *status = ptwX_norm->status;
        if( ptwX_norm->status != nfu_Okay ) return( NULL );
        if( ptwX_length( ptwX_norm ) != ngs ) {
            *status = nfu_badNorm;
            return( NULL );
        }
    }
 
    if( ( f = ptwXY_intersectionWith_ptwX( ptwXY, groupBoundaries, status ) ) == NULL ) return( NULL );
    if( f->length == 0 ) return( ptwX_createLine( ngs, ngs, 0, 0, status ) );
 
    if( ( groupedData = ptwX_new( ngs, status ) ) == NULL ) goto err;
    xg1 = groupBoundaries->points[0];
    x1 = f->points[0].x;
    y1 = f->points[0].y;
    for( igs = 0, i = 1; igs < ngs; igs++ ) {
        xg2 = groupBoundaries->points[igs+1];
        sum = 0;
        if( xg2 > x1 ) {
            for( ; i < f->length; i++, x1 = x2, y1 = y2 ) {
                x2 = f->points[i].x;
                if( x2 > xg2 ) break;
                y2p = y2 = f->points[i].y;
                if( f->interpolation == ptwXY_interpolationFlat ) y2p = y1;
                sum += ( y1 + y2p ) * ( x2 - x1 );
            }
        }
        if( sum != 0. ) {
            if( normType == ptwXY_group_normType_dx ) {
                sum /= ( xg2 - xg1 ); }
            else if( normType == ptwXY_group_normType_norm ) {
                if( ptwX_norm->points[igs] == 0. ) {
                    *status = nfu_divByZero;
                    goto err;
                }
                sum /= ptwX_norm->points[igs];
            }
        }
        groupedData->points[igs] = 0.5 * sum;
        groupedData->length++;
        xg1 = xg2;
    }
 
    ptwXY_free( f );
    return( groupedData );
 
err:
    ptwXY_free( f );
    if( groupedData != NULL ) ptwX_free( groupedData );
    return( NULL );
}

ptwXY_groupThreeFunctions()

ptwXPoints * ptwXY_groupThreeFunctions	(	ptwXYPoints *	ptwXY1,
		ptwXYPoints *	ptwXY2,
		ptwXYPoints *	ptwXY3,
		ptwXPoints *	groupBoundaries,
		ptwXY_group_normType	normType,
		ptwXPoints *	ptwX_norm,
		nfu_status *	status )

Definition at line 528 of file ptwXY_integration.cc.

                                                                                   {
 
    int64_t i, igs, ngs;
    double x1, fy1, gy1, hy1, x2, fy2, gy2, hy2, fy2p, gy2p, hy2p, xg1, xg2, sum;
    ptwXYPoints *f = NULL, *ff, *fff = NULL, *g = NULL, *gg = NULL, *h = NULL, *hh = NULL;
    ptwXPoints *groupedData = NULL;
 
    if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
    if( ( *status = ptwXY_simpleCoalescePoints( ptwXY2 ) ) != nfu_Okay ) return( NULL );
    if( ( *status = ptwXY_simpleCoalescePoints( ptwXY3 ) ) != nfu_Okay ) return( NULL );
    if( ( *status = groupBoundaries->status ) != nfu_Okay ) return( NULL );
    *status = nfu_otherInterpolation;
    if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
    if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( NULL );
    if( ptwXY3->interpolation == ptwXY_interpolationOther ) return( NULL );
 
    ngs = ptwX_length( groupBoundaries ) - 1;
    if( normType == ptwXY_group_normType_norm ) {
        if( ptwX_norm == NULL ) {
            *status = nfu_badNorm;
            return( NULL );
        }
        if( ( *status = ptwX_norm->status ) != nfu_Okay ) return( NULL );
        if( ptwX_length( ptwX_norm ) != ngs ) {
            *status = nfu_badNorm;
            return( NULL );
        }
    }
 
    if( ( ff = ptwXY_intersectionWith_ptwX( ptwXY1, groupBoundaries, status ) ) == NULL ) return( NULL );
    if( ( gg = ptwXY_intersectionWith_ptwX( ptwXY2, groupBoundaries, status ) ) == NULL ) goto err;
    if( ( hh = ptwXY_intersectionWith_ptwX( ptwXY3, groupBoundaries, status ) ) == NULL ) goto err;
    if( ( ff->length == 0 ) || ( gg->length == 0 ) || ( hh->length == 0 ) ) return( ptwX_createLine( ngs, ngs, 0, 0, status ) );
 
    if( ( *status = ptwXY_tweakDomainsToMutualify( ff, gg, 4, 0 ) ) != nfu_Okay ) goto err;
    if( ( *status = ptwXY_tweakDomainsToMutualify( ff, hh, 4, 0 ) ) != nfu_Okay ) goto err;
    if( ( *status = ptwXY_tweakDomainsToMutualify( gg, hh, 4, 0 ) ) != nfu_Okay ) goto err;
    if( ( fff = ptwXY_union(  ff,  gg, status, ptwXY_union_fill ) ) == NULL ) goto err;
    if( (   h = ptwXY_union(  hh, fff, status, ptwXY_union_fill ) ) == NULL ) goto err;
    if( (   f = ptwXY_union( fff,   h, status, ptwXY_union_fill ) ) == NULL ) goto err;
    if( (   g = ptwXY_union(  gg,   h, status, ptwXY_union_fill ) ) == NULL ) goto err;
 
    if( ( groupedData = ptwX_new( ngs, status ) ) == NULL ) goto err;
    xg1 = groupBoundaries->points[0];
    x1 = f->points[0].x;
    fy1 = f->points[0].y;
    gy1 = g->points[0].y;
    hy1 = h->points[0].y;
    for( igs = 0, i = 1; igs < ngs; igs++ ) {
        xg2 = groupBoundaries->points[igs+1];
        sum = 0;
        if( xg2 > x1 ) {
            for( ; i < f->length; i++, x1 = x2, fy1 = fy2, gy1 = gy2, hy1 = hy2 ) {
                x2 = f->points[i].x;
                if( x2 > xg2 ) break;
                fy2p = fy2 = f->points[i].y;
                if( f->interpolation == ptwXY_interpolationFlat ) fy2p = fy1;
                gy2p = gy2 = g->points[i].y;
                if( g->interpolation == ptwXY_interpolationFlat ) gy2p = gy1;
                hy2p = hy2 = h->points[i].y;
                if( h->interpolation == ptwXY_interpolationFlat ) hy2p = hy1;
                sum += ( ( fy1 + fy2p ) * ( gy1 + gy2p ) * ( hy1 + hy2p ) + 2 * fy1 * gy1 * hy1 + 2 * fy2p * gy2p * hy2p ) * ( x2 - x1 );
            }
        }
        if( sum != 0. ) {
            if( normType == ptwXY_group_normType_dx ) {
                sum /= ( xg2 - xg1 ); }
            else if( normType == ptwXY_group_normType_norm ) {
                if( ptwX_norm->points[igs] == 0. ) {
                    *status = nfu_divByZero;
                    goto err;
                }
                sum /= ptwX_norm->points[igs];
            }
        }
        groupedData->points[igs] = sum / 12.;
        groupedData->length++;
        xg1 = xg2;
    }
 
    ptwXY_free( f );
    ptwXY_free( g );
    ptwXY_free( h );
    ptwXY_free( ff );
    ptwXY_free( gg );
    ptwXY_free( hh );
    ptwXY_free( fff );
    return( groupedData );
 
err:
    ptwXY_free( ff );
    if( fff != NULL ) ptwXY_free( fff );
    if( gg != NULL ) ptwXY_free( gg );
    if( hh != NULL ) ptwXY_free( hh );
    if( f != NULL ) ptwXY_free( f );
    if( g != NULL ) ptwXY_free( g );
    if( h != NULL ) ptwXY_free( h );
    if( groupedData != NULL ) ptwX_free( groupedData );
    return( NULL );
}

ptwXY_groupTwoFunctions()

ptwXPoints * ptwXY_groupTwoFunctions	(	ptwXYPoints *	ptwXY1,
		ptwXYPoints *	ptwXY2,
		ptwXPoints *	groupBoundaries,
		ptwXY_group_normType	normType,
		ptwXPoints *	ptwX_norm,
		nfu_status *	status )

Definition at line 435 of file ptwXY_integration.cc.

                                                    {
 
    int64_t i, igs, ngs;
    double x1, fy1, gy1, x2, fy2, gy2, fy2p, gy2p, xg1, xg2, sum;
    ptwXYPoints *f = NULL, *ff, *g = NULL, *gg = NULL;
    ptwXPoints *groupedData = NULL;
 
    if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
    if( ( *status = ptwXY_simpleCoalescePoints( ptwXY2 ) ) != nfu_Okay ) return( NULL );
    if( ( *status = groupBoundaries->status ) != nfu_Okay ) return( NULL );
    *status = nfu_otherInterpolation;
    if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
    if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( NULL );
 
    ngs = ptwX_length( groupBoundaries ) - 1;
    if( normType == ptwXY_group_normType_norm ) {
        if( ptwX_norm == NULL ) {
            *status = nfu_badNorm;
            return( NULL );
        }
        if( ( *status = ptwX_norm->status ) != nfu_Okay ) return( NULL );
        if( ptwX_length( ptwX_norm ) != ngs ) {
            *status = nfu_badNorm;
            return( NULL );
        }
    }
 
    if( ( ff = ptwXY_intersectionWith_ptwX( ptwXY1, groupBoundaries, status ) ) == NULL ) return( NULL );
    if( ( gg = ptwXY_intersectionWith_ptwX( ptwXY2, groupBoundaries, status ) ) == NULL ) goto err;
    if( ( ff->length == 0 ) || ( gg->length == 0 ) ) {
        ptwXY_free( ff );
        ptwXY_free( gg );
        return( ptwX_createLine( ngs, ngs, 0, 0, status ) );
    }
 
    if( ( *status = ptwXY_tweakDomainsToMutualify( ff, gg, 4, 0 ) ) != nfu_Okay ) goto err;
    if( ( f = ptwXY_union( ff, gg, status, ptwXY_union_fill ) ) == NULL ) goto err;
    if( ( g = ptwXY_union( gg, f, status, ptwXY_union_fill ) ) == NULL ) goto err;
 
    if( ( groupedData = ptwX_new( ngs, status ) ) == NULL ) goto err;
    xg1 = groupBoundaries->points[0];
    x1 = f->points[0].x;
    fy1 = f->points[0].y;
    gy1 = g->points[0].y;
    for( igs = 0, i = 1; igs < ngs; igs++ ) {
        xg2 = groupBoundaries->points[igs+1];
        sum = 0;
        if( xg2 > x1 ) {
            for( ; i < f->length; i++, x1 = x2, fy1 = fy2, gy1 = gy2 ) {
                x2 = f->points[i].x;
                if( x2 > xg2 ) break;
                fy2p = fy2 = f->points[i].y;
                if( f->interpolation == ptwXY_interpolationFlat ) fy2p = fy1;
                gy2p = gy2 = g->points[i].y;
                if( g->interpolation == ptwXY_interpolationFlat ) gy2p = gy1;
                sum += ( ( fy1 + fy2p ) * ( gy1 + gy2p ) + fy1 * gy1 + fy2p * gy2p ) * ( x2 - x1 );
            }
        }
        if( sum != 0. ) {
            if( normType == ptwXY_group_normType_dx ) {
                sum /= ( xg2 - xg1 ); }
            else if( normType == ptwXY_group_normType_norm ) {
                if( ptwX_norm->points[igs] == 0. ) {
                    *status = nfu_divByZero;
                    goto err;
                }
                sum /= ptwX_norm->points[igs];
            }
        }
        groupedData->points[igs] = sum / 6.;
        groupedData->length++;
        xg1 = xg2;
    }
 
    ptwXY_free( f );
    ptwXY_free( g );
    ptwXY_free( ff );
    ptwXY_free( gg );
    return( groupedData );
 
err:
    ptwXY_free( ff );
    if( gg != NULL ) ptwXY_free( gg );
    // Coverity #63063
    if( f != NULL ) ptwXY_free( f );
    if( g != NULL ) ptwXY_free( g );
    if( groupedData != NULL ) ptwX_free( groupedData );
    return( NULL );
}

Referenced by GIDI_settings_processedFlux::groupFunction().

ptwXY_integrate()

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

Definition at line 118 of file ptwXY_integration.cc.

                                                                                           {
 
    int64_t i, n = ptwXY->length;
    double sum = 0., dSum, x, y, x1, x2, y1, y2, _sign = 1.;
    ptwXYPoint *point;
 
    if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
    *status = nfu_otherInterpolation;
    if( ptwXY->interpolation == ptwXY_interpolationOther ) return( 0. );
 
    if( xMax < xMin ) {
        x = xMin;
        xMin = xMax;
        xMax = x;
        _sign = -1.;
    }
    if( n < 2 ) return( 0. );
 
    if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( 0. );
    for( i = 0, point = ptwXY->points; i < n; i++, point++ ) {
        if( point->x >= xMin ) break;
    }
    if( i == n ) return( 0. );
    x2 = point->x;
    y2 = point->y;
    if( i > 0 ) {
        if( x2 > xMin ) {
            x1 = point[-1].x;
            y1 = point[-1].y;
            if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMin, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( 0. );
            if( x2 > xMax ) {
                double yMax;
 
                if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMax, &yMax, x1, y1, x2, y2 ) ) != nfu_Okay ) return( 0. );
                if( ( *status = ptwXY_f_integrate( ptwXY->interpolation, xMin, y, xMax, yMax, &sum ) ) != nfu_Okay ) return( 0. );
                return( sum ); }
            else {
                if( ( *status = ptwXY_f_integrate( ptwXY->interpolation, xMin, y, x2, y2, &sum ) ) != nfu_Okay ) return( 0. );
            }
        }
    }
    i++;
    point++;
    for( ; i < n; i++, point++ ) {
        x1 = x2;
        y1 = y2;
        x2 = point->x;
        y2 = point->y;
        if( x2 > xMax ) {
            if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMax, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( 0. );
            if( ( *status = ptwXY_f_integrate( ptwXY->interpolation, x1, y1, xMax, y, &dSum ) ) != nfu_Okay ) return( 0. );
            sum += dSum;
            break;
        }
        if( ( *status = ptwXY_f_integrate( ptwXY->interpolation, x1, y1, x2, y2, &dSum ) ) != nfu_Okay ) return( 0. );
        sum += dSum;
    }
 
    return( _sign * sum );
}

Referenced by ptwXY_integrateDomain().

ptwXY_integrateDomain()

double ptwXY_integrateDomain	(	ptwXYPoints *	ptwXY,
		nfu_status *	status )

Definition at line 181 of file ptwXY_integration.cc.

                                                                       {
 
    if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
    if( ptwXY->length > 0 ) return( ptwXY_integrate( ptwXY, ptwXY_getXMin( ptwXY ), ptwXY_getXMax( ptwXY ), status ) );
    return( 0. );
}

Referenced by ptwXY_normalize().

ptwXY_integrateDomainWithWeight_sqrt_x()

double ptwXY_integrateDomainWithWeight_sqrt_x	(	ptwXYPoints *	ptwXY,
		nfu_status *	status )

Definition at line 284 of file ptwXY_integration.cc.

                                                                                        {
 
    if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
    if( ptwXY->length < 2 ) return( 0. );
    return( ptwXY_integrateWithWeight_sqrt_x( ptwXY, ptwXY_getXMin( ptwXY ), ptwXY_getXMax( ptwXY ), status ) );
}

ptwXY_integrateDomainWithWeight_x()

double ptwXY_integrateDomainWithWeight_x	(	ptwXYPoints *	ptwXY,
		nfu_status *	status )

Definition at line 210 of file ptwXY_integration.cc.

                                                                                   {
 
    if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
    if( ptwXY->length < 2 ) return( 0. );
    return( ptwXY_integrateWithWeight_x( ptwXY, ptwXY_getXMin( ptwXY ), ptwXY_getXMax( ptwXY ), status ) );
}

ptwXY_integrateWithFunction()

double ptwXY_integrateWithFunction	(	ptwXYPoints *	ptwXY,
		ptwXY_createFromFunction_callback	func,
		void *	argList,
		double	xMin,
		double	xMax,
		int	degree,
		int	recursionLimit,
		double	tolerance,
		nfu_status *	status )

Definition at line 657 of file ptwXY_integration.cc.

                                                                                                         {
 
    int64_t i1, i2, n1 = ptwXY->length;
    long evaluations;
    double integral = 0., integral_, sign = -1., xa, xb;
    ptwXY_integrateWithFunctionInfo integrateWithFunctionInfo;
    ptwXYPoint *point;
 
    if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
 
    if( xMin == xMax ) return( 0. );
    if( n1 < 2 ) return( 0. );
 
    ptwXY_simpleCoalescePoints( ptwXY );
 
    if( xMin > xMax ) {
        sign = xMin;
        xMin = xMax;
        xMax = sign;
        sign = -1.;
    }
    if( xMin >= ptwXY->points[n1-1].x ) return( 0. );
    if( xMax <= ptwXY->points[0].x ) return( 0. );
 
    for( i1 = 0; i1 < ( n1 - 1 ); i1++ ) {
        if( ptwXY->points[i1+1].x > xMin ) break;
    }
    for( i2 = n1 - 1; i2 > i1; i2-- ) {
        if( ptwXY->points[i2-1].x < xMax ) break;
    }
    point = &(ptwXY->points[i1]);
 
    integrateWithFunctionInfo.degree = degree;
    integrateWithFunctionInfo.func = func;
    integrateWithFunctionInfo.argList = argList;
    integrateWithFunctionInfo.interpolation = ptwXY->interpolation;
    integrateWithFunctionInfo.x2 = point->x;
    integrateWithFunctionInfo.y2 = point->y;
 
    xa = xMin;
    for( ; i1 < i2; i1++ ) {
        integrateWithFunctionInfo.x1 = integrateWithFunctionInfo.x2;
        integrateWithFunctionInfo.y1 = integrateWithFunctionInfo.y2;
        ++point;
        integrateWithFunctionInfo.x2 = point->x;
        integrateWithFunctionInfo.y2 = point->y;
        xb = point->x;
        if( xb > xMax ) xb = xMax;
        *status = nf_GnG_adaptiveQuadrature( ptwXY_integrateWithFunction2, ptwXY_integrateWithFunction3, &integrateWithFunctionInfo,
            xa, xb, recursionLimit, tolerance, &integral_, &evaluations );
        if( *status != nfu_Okay ) return( 0. );
        integral += integral_;
        xa = xb;
    }
 
    return( integral );
}

ptwXY_integrateWithWeight_sqrt_x()

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

Definition at line 293 of file ptwXY_integration.cc.

                                                                                                            {
 
    int64_t i, n = ptwXY->length;
    double sum = 0., x, y, x1, x2, y1, y2, _sign = 1., sqrt_x1, sqrt_x2, inv_apb, c;
    ptwXYPoint *point;
 
    if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
    *status = nfu_unsupportedInterpolation;
    if( ( ptwXY->interpolation != ptwXY_interpolationLinLin ) &&
        ( ptwXY->interpolation != ptwXY_interpolationFlat ) ) return( 0. );
 
    if( n < 2 ) return( 0. );
    if( xMax < xMin ) {
        x = xMin;
        xMin = xMax;
        xMax = x;
        _sign = -1.;
    }
 
    if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( 0. );
    for( i = 0, point = ptwXY->points; i < n; ++i, ++point ) {
        if( point->x >= xMin ) break;
    }
    if( i == n ) return( 0. );
    x2 = point->x;
    y2 = point->y;
    if( i > 0 ) {
        if( x2 > xMin ) {
            if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMin, &y, point[-1].x, point[-1].y, x2, y2 ) ) != nfu_Okay ) return( 0. );
            x2 = xMin;
            y2 = y;
            --i;
            --point;
        }
    }
    ++i;
    ++point;
    sqrt_x2 = std::sqrt( x2 );
    for( ; i < n; ++i, ++point ) {
        x1 = x2;
        y1 = y2;
        sqrt_x1 = sqrt_x2;
        x2 = point->x;
        y2 = point->y;
        if( x2 > xMax ) {
            if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMax, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( 0. );
            x2 = xMax;
            y2 = y;
        }
        sqrt_x2 = std::sqrt( x2 );
        inv_apb = sqrt_x1 + sqrt_x2;
        c = 2. * ( sqrt_x1 * sqrt_x2 + x1 + x2 );
        switch( ptwXY->interpolation ) {
        case ptwXY_interpolationFlat :
            sum += ( sqrt_x2 - sqrt_x1 ) * y1 * 2.5 * c;
            break;
        case ptwXY_interpolationLinLin :
            sum += ( sqrt_x2 - sqrt_x1 ) * ( y1 * ( c + x1 * ( 1. + sqrt_x2 / inv_apb ) ) + y2 * ( c + x2 * ( 1. + sqrt_x1 / inv_apb ) ) );
            break;
        default :       /* Only to stop compilers from complaining. */
            break;
        }
        if( x2 == xMax ) break;
    }
 
    return( 2. / 15. * _sign * sum );
}

Referenced by ptwXY_integrateDomainWithWeight_sqrt_x().

ptwXY_integrateWithWeight_x()

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

Definition at line 219 of file ptwXY_integration.cc.

                                                                                                       {
 
    int64_t i, n = ptwXY->length;
    double sum = 0., x, y, x1, x2, y1, y2, _sign = 1.;
    ptwXYPoint *point;
 
    if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
    *status = nfu_unsupportedInterpolation;
    if( ( ptwXY->interpolation != ptwXY_interpolationLinLin ) && 
        ( ptwXY->interpolation != ptwXY_interpolationFlat ) ) return( 0. );
 
    if( n < 2 ) return( 0. );
    if( xMax < xMin ) {
        x = xMin;
        xMin = xMax;
        xMax = x;
        _sign = -1.;
    }
 
    if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( 0. );
    for( i = 0, point = ptwXY->points; i < n; ++i, ++point ) {
        if( point->x >= xMin ) break;
    }
    if( i == n ) return( 0. );
    x2 = point->x;
    y2 = point->y;
    if( i > 0 ) {
        if( x2 > xMin ) {
            if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMin, &y, point[-1].x, point[-1].y, x2, y2 ) ) != nfu_Okay ) return( 0. );
            x2 = xMin;
            y2 = y;
            --i;
            --point;
        }
    }
    ++i;
    ++point;
    for( ; i < n; ++i, ++point ) {
        x1 = x2;
        y1 = y2;
        x2 = point->x;
        y2 = point->y;
        if( x2 > xMax ) {
            if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMax, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( 0. );
            x2 = xMax;
            y2 = y;
        }
        switch( ptwXY->interpolation ) {
        case ptwXY_interpolationFlat :
            sum += ( x2 - x1 ) * y1 * 3 * ( x1 + x2 );
            break;
        case ptwXY_interpolationLinLin :
            sum += ( x2 - x1 ) * ( y1 * ( 2 * x1 + x2 ) + y2 * ( x1 + 2 * x2 ) );
            break;
        default :       /* Only to stop compilers from complaining. */
            break;
        }
        if( x2 == xMax ) break;
    }
 
    return( _sign * sum / 6 );
}

Referenced by ptwXY_integrateDomainWithWeight_x().

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_intersectionWith_ptwX()

ptwXYPoints * ptwXY_intersectionWith_ptwX	(	ptwXYPoints *	ptwXY,
		ptwXPoints *	ptwX,
		nfu_status *	status )

Definition at line 194 of file ptwXY_convenient.cc.

                                                                                                     {
 
    int64_t i, i1, i2, lengthX = ptwX_length( ptwX );
    double x, y, xMin, xMax;
    ptwXYPoints *n = NULL;
 
    if( ( *status = ptwXY->status ) != nfu_Okay ) return( NULL );
    if( ( *status = ptwX->status ) != nfu_Okay ) return( NULL );
    if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) goto Err;
    *status = nfu_otherInterpolation;
    if( ptwXY->interpolation == ptwXY_interpolationOther ) return( NULL );
    if( ( n = ptwXY_clone( ptwXY, status ) ) == NULL ) return( NULL );
    if( ptwXY->length == 0 ) return( n );
    xMin = ptwXY->points[0].x;
    xMax = ptwXY->points[ptwXY->length - 1].x;
 
    if( ( xMin >= ptwX->points[lengthX-1] ) || ( xMax <= ptwX->points[0] ) ) {  /* No overlap. */
        n->length = 0;
        return( n );
    }
 
    for( i = 0; i < lengthX; i++ ) {        /* Fill in ptwXY at x-points in ptwX. */
        x = ptwX->points[i];
        if( x <= xMin ) continue;
        if( x >= xMax ) break;
        if( ( *status = ptwXY_getValueAtX( ptwXY, x, &y ) ) != nfu_Okay ) goto Err;
        if( ( *status = ptwXY_setValueAtX( n, x, y ) ) != nfu_Okay ) goto Err;
    }
    if( ( *status = ptwXY_simpleCoalescePoints( n ) ) != nfu_Okay ) goto Err;
 
    i1 = 0;
    i2 = n->length - 1;
    if( lengthX > 0 ) {
        x = ptwX->points[0];
        if( x > n->points[i1].x ) {
            for( ; i1 < n->length; i1++ ) {
                if( n->points[i1].x == x ) break;
            }
        }
 
        x = ptwX->points[lengthX - 1];
        if( x < n->points[i2].x ) {
            for( ; i2 > i1; i2-- ) {
                if( n->points[i2].x == x ) break;
            }
        }
    }
    i2++;
 
    if( i1 != 0 ) {
        for( i = i1; i < i2; i++ ) n->points[i - i1] = n->points[i];
    }
    n->length = i2 - i1;
 
    return( n );
 
Err:
     ptwXY_free( n );
     return( NULL );
}

Referenced by ptwXY_groupOneFunction(), ptwXY_groupThreeFunctions(), and ptwXY_groupTwoFunctions().

ptwXY_length()

int64_t ptwXY_length

(

ptwXYPoints *

ptwXY

)

Definition at line 583 of file ptwXY_core.cc.

                                           {
 
    return( ptwXY->length );
}

Referenced by MCGIDI_angular_parseFromTOM(), MCGIDI_fromTOM_pdfOfX(), MCGIDI_product_parseFromTOM(), nf_Legendre_from_ptwXY(), and ptwXY_valueTo_ptwXAndY().

ptwXY_mergeClosePoints()

nfu_status ptwXY_mergeClosePoints	(	ptwXYPoints *	ptwXY,
		double	epsilon )

Definition at line 141 of file ptwXY_convenient.cc.

                                                                        {
 
    int64_t i, i1, j, k, n = ptwXY->length;
    double x, y;
    ptwXYPoint *p1, *p2;
 
    if( n < 2 ) return( ptwXY->status );
    if( epsilon < 4 * DBL_EPSILON ) epsilon = 4 * DBL_EPSILON;
    if( ptwXY_simpleCoalescePoints( ptwXY ) != nfu_Okay ) return( ptwXY->status );
 
    p2 = ptwXY->points;
    x = p2->x;
    for( i1 = 1, p2++; i1 < ( n - 1 ); i1++, p2++ ) {                 /* The first point shall remain the first point and all points close to it are deleted. */
        if( ( p2->x - x ) > 0.5 * epsilon * ( std::fabs( x ) + std::fabs( p2->x ) ) ) break;
    }
    if( i1 != 1 ) {
        for( i = i1, p1 = &(ptwXY->points[1]); i < n; i++, p1++, p2++ ) *p1 = *p2;
        n = ptwXY->length = ptwXY->length - i1 + 1;
    }
 
    p1 = &(ptwXY->points[n-1]);
    x = p1->x;
    for( i1 = n - 2, p1--; i1 > 0; i1--, p1-- ) {            /* The last point shall remain the last point and all points close to it are deleted. */
        if( x - p1->x > 0.5 * epsilon * ( std::fabs( x ) + std::fabs( p1->x ) ) ) break;
    }
    if( i1 != ( n - 2 ) ) {
        ptwXY->points[i1 + 1] = ptwXY->points[n - 1];
        n = ptwXY->length = i1 + 2;
    }
 
    for( i = 1; i < n - 1; i++ ) {
        p1 = &(ptwXY->points[i]);
        x = p1->x;
        y = p1->y;
        for( j = i + 1, p2 = &(ptwXY->points[i+1]); j < n - 1; j++, p2++ ) {
            if( ( p2->x - p1->x ) > 0.5 * epsilon * ( std::fabs( p2->x ) + std::fabs( p1->x ) ) ) break;
            x += p2->x;
            y += p2->y;
        }
        if( ( k = ( j - i ) ) > 1 ) {
            p1->x = x / k;
            p1->y = y / k;
            for( p1 = &(ptwXY->points[i+1]); j < n; j++, p1++, p2++ ) *p1 = *p2;
            n -= ( k - 1 );
        }
    }
    ptwXY->length = n;
 
    return( ptwXY->status );
}

Referenced by ptwXY_union().

ptwXY_mergeFromXsAndYs()

nfu_status ptwXY_mergeFromXsAndYs	(	ptwXYPoints *	ptwXY,
		int	length,
		double *	xs,
		double *	ys )

Definition at line 966 of file ptwXY_core.cc.

                                                                                            {
 
    return( ptwXY_mergeFrom( ptwXY, 1, length, xs, ys ) );
}

ptwXY_mergeFromXYs()

nfu_status ptwXY_mergeFromXYs	(	ptwXYPoints *	ptwXY,
		int	length,
		double *	xys )

Definition at line 973 of file ptwXY_core.cc.

                                                                             {
 
    int i;
    double *xs, *p1, *p2;
    nfu_status status;
 
    if( length < 0 ) return( nfu_badInput );
    if( length == 0 ) return( nfu_Okay );
    if( ( xs = (double *) nfu_malloc( length * sizeof( double ) ) ) == NULL ) return( nfu_mallocError );
    for( i = 0, p1 = xs, p2 = xys; i < length; i++, p1++, p2 += 2 ) *p1 = *p2;
    status = ptwXY_mergeFrom( ptwXY, 2, length, xs, xys );
    nfu_free( xs );
 
    return( status );
}

ptwXY_mod()

nfu_status ptwXY_mod	(	ptwXYPoints *	ptwXY,
		double	m,
		int	pythonMod )

Definition at line 76 of file ptwXY_binaryOperators.cc.

                                                                    { 
 
    int64_t i, nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
    ptwXYPoint *p;
    ptwXYOverflowPoint *o, *overflowHeader = &(ptwXY->overflowHeader);
 
    if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
    if( m == 0 ) return( ptwXY->status = nfu_divByZero );
 
    for( i = 0, p = ptwXY->points; i < nonOverflowLength; i++, p++ ) p->y = ptwXY_mod2( p->y, m, pythonMod );
    for( o = overflowHeader->next; o != overflowHeader; o = o->next ) o->point.y = ptwXY_mod2( o->point.y, m, pythonMod );
    return( ptwXY->status );
}

ptwXY_mul2_ptwXY()

ptwXYPoints * ptwXY_mul2_ptwXY	(	ptwXYPoints *	ptwXY1,
		ptwXYPoints *	ptwXY2,
		nfu_status *	status )

Definition at line 187 of file ptwXY_binaryOperators.cc.

                                                                                              {
 
    int64_t i, length;
    ptwXYPoints *n = NULL;
    int found;
    double x1, y1, x2, y2, u1, u2, v1, v2, xz1 = 0, xz2 = 0, x;
 
    *status = nfu_otherInterpolation;
    if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
    if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( NULL );
    if( ( n = ptwXY_mul_ptwXY( ptwXY1, ptwXY2, status ) ) == NULL ) return( n );
    if( ptwXY1->interpolation == ptwXY_interpolationFlat ) return( n );
    if( ptwXY2->interpolation == ptwXY_interpolationFlat ) return( n );
    length = n->length - 1;
    if( length > 0 ) {
        x2 = n->points[length].x;
        for( i = length - 1; i >= 0; i-- ) {             /* Find and add y zeros not currently in n's. */
            x1 = n->points[i].x;
            if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x1, &u1 ) ) != nfu_Okay ) goto Err;
            if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x2, &u2 ) ) != nfu_Okay ) goto Err;
            if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, x1, &v1 ) ) != nfu_Okay ) goto Err;
            if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, x2, &v2 ) ) != nfu_Okay ) goto Err;
            found = 0;
            if( u1 * u2 < 0 ) {
                xz1 = ( u1 * x2 - u2 * x1 ) / ( u1 - u2 );
                if( ( *status = ptwXY_setValueAtX( n, xz1, 0. ) ) != nfu_Okay ) goto Err;
                found = 1;
            }
            if( v1 * v2 < 0 ) {
                xz2 = ( v1 * x2 - v2 * x1 ) / ( v1 - v2 );
                if( ( *status = ptwXY_setValueAtX( n, xz2, 0. ) ) != nfu_Okay ) goto Err;
                found += 1;
            }
            if( found > 1 ) {
                x = 0.5 * ( xz1 + xz2 );
                if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x, &u1 ) ) != nfu_Okay ) goto Err;
                if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, x, &v1 ) ) != nfu_Okay ) goto Err;
                if( ( *status = ptwXY_setValueAtX( n, x, u1 * v1 ) ) != nfu_Okay ) goto Err;
            }
            x2 = x1;
        }
 
        if( ( *status = ptwXY_simpleCoalescePoints( n ) ) != nfu_Okay ) goto Err;
        length = n->length;
        x2 = n->points[n->length-1].x;
        y2 = n->points[n->length-1].y;
        for( i = n->length - 2; i >= 0; i-- ) {             /* Make interpolation fit accuracy. Work backwards so new points will not mess up loop. */
            x1 = n->points[i].x;
            y1 = n->points[i].y;
            if( ( *status = ptwXY_mul2_s_ptwXY( n, ptwXY1, ptwXY2, x1, y1, x2, y2, 0 ) ) != nfu_Okay ) goto Err;
            x2 = x1;
            y2 = y1;
        }
        ptwXY_update_biSectionMax( n, (double) length );
    }
    return( n );
 
Err:
    if( n ) ptwXY_free( n );
    return( NULL );
}

ptwXY_mul_double()

nfu_status ptwXY_mul_double	(	ptwXYPoints *	ptwXY,
		double	value )

Definition at line 43 of file ptwXY_binaryOperators.cc.

{ return( ptwXY_slopeOffset( ptwXY, value, 0. ) ); }

ptwXY_mul_ptwXY()

ptwXYPoints * ptwXY_mul_ptwXY	(	ptwXYPoints *	ptwXY1,
		ptwXYPoints *	ptwXY2,
		nfu_status *	status )

Definition at line 171 of file ptwXY_binaryOperators.cc.

                                                                                             {
 
    ptwXYPoints *mul;
 
    if( ptwXY1->length == 0 ) {
        mul = ptwXY_clone( ptwXY1, status ); }
    else if( ptwXY2->length == 0 ) {
        mul = ptwXY_clone( ptwXY2, status ); }
    else {
        mul = ptwXY_binary_ptwXY( ptwXY1, ptwXY2, 0., 0., 1., status );
    }
    return( mul );
}

Referenced by ptwXY_mul2_ptwXY().

ptwXY_mutualifyDomains()

nfu_status ptwXY_mutualifyDomains	(	ptwXYPoints *	ptwXY1,
		double	lowerEps1,
		double	upperEps1,
		int	positiveXOnly1,
		ptwXYPoints *	ptwXY2,
		double	lowerEps2,
		double	upperEps2,
		int	positiveXOnly2 )

Definition at line 368 of file ptwXY_convenient.cc.

                                                                                                                        {
 
    nfu_status status;
    int64_t n1 = ptwXY1->length, n2 = ptwXY2->length;
    ptwXYPoint *xy1, *xy2;
 
    switch( status = ptwXY_areDomainsMutual( ptwXY1, ptwXY2 ) ) {
    case nfu_Okay :
    case nfu_empty :
        return( nfu_Okay );
    case nfu_domainsNotMutual :
        break;
    default :
        return( status );
    }
 
    if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
    if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
    if( ptwXY1->interpolation == ptwXY_interpolationFlat ) return( nfu_invalidInterpolation );
    if( ptwXY2->interpolation == ptwXY_interpolationFlat ) return( nfu_invalidInterpolation );
 
    xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, 0 );
    xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, 0 );
    if( xy1->x < xy2->x ) {
        lowerEps1 = 0.;
        if( xy2->y == 0. ) lowerEps2 = 0.; }
    else if( xy1->x > xy2->x ) {
        lowerEps2 = 0.;
        if( xy1->y == 0. ) lowerEps1 = 0.; }
    else {
        lowerEps1 = lowerEps2 = 0.;
    }
 
    xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, n1 - 1 );
    xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, n2 - 1 );
    if( xy1->x < xy2->x ) {
        upperEps2 = 0.;
        if( xy1->y == 0. ) upperEps1 = 0.; }
    else if( xy1->x > xy2->x ) {
        upperEps1 = 0.;
        if( xy2->y == 0. ) upperEps2 = 0.; }
    else {
        upperEps1 = upperEps2 = 0.;
    }
 
    if( ( lowerEps1 != 0. ) || ( upperEps1 != 0. ) ) 
        if( ( status = ptwXY_dullEdges( ptwXY1, lowerEps1, upperEps1, positiveXOnly1 ) ) != nfu_Okay ) return( status );
    if( ( lowerEps2 != 0. ) || ( upperEps2 != 0. ) ) 
        if( ( status = ptwXY_dullEdges( ptwXY2, lowerEps2, upperEps2, positiveXOnly2 ) ) != nfu_Okay ) return( status );
 
    return( status );
}

ptwXY_neg()

nfu_status ptwXY_neg

(

ptwXYPoints *

ptwXY

)

Definition at line 34 of file ptwXY_unitaryOperators.cc.

                                           {
 
    int64_t i, nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
    ptwXYPoint *p;
    ptwXYOverflowPoint *o, *overflowHeader = &(ptwXY->overflowHeader);
 
    if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
 
    for( i = 0, p = ptwXY->points; i < nonOverflowLength; i++, p++ ) p->y = -p->y;
    for( o = overflowHeader->next; o != overflowHeader; o = o->next ) o->point.y = -o->point.y;
    return( ptwXY->status );
}

Referenced by ptwXY_sub_ptwXY().

ptwXY_new()

ptwXYPoints * ptwXY_new	(	ptwXY_interpolation	interpolation,
		ptwXY_interpolationOtherInfo const *	interpolationOtherInfo,
		double	biSectionMax,
		double	accuracy,
		int64_t	primarySize,
		int64_t	secondarySize,
		nfu_status *	status,
		int	userFlag )

Definition at line 29 of file ptwXY_core.cc.

                                                                                                    {
 
    ptwXYPoints *ptwXY = (ptwXYPoints *) nfu_calloc( sizeof( ptwXYPoints ), 1 );
 
    *status = nfu_mallocError;
    if( ptwXY == NULL ) return( NULL );
    ptwXY_setup( ptwXY, interpolation, interpolationOtherInfo, biSectionMax, accuracy, primarySize, 
        secondarySize, userFlag );
    if( ( *status = ptwXY->status ) != nfu_Okay ) {
        ptwXY = (ptwXYPoints *) nfu_free( ptwXY );
    }
    return( ptwXY );
}

Referenced by MCGIDI_product_parseFromTOM(), ptwXY_clip(), ptwXY_convolution(), ptwXY_create(), ptwXY_createFrom_Xs_Ys(), ptwXY_createFromFunction(), ptwXY_createGaussianCenteredSigma1(), ptwXY_flatInterpolationToLinear(), ptwXY_slice(), ptwXY_thin(), ptwXY_union(), ptwXY_valueTo_ptwXY(), and ptwXY_xSlice().

ptwXY_normalize()

nfu_status ptwXY_normalize

(

ptwXYPoints *

ptwXY1

)

Definition at line 190 of file ptwXY_integration.cc.

                                                 {
/*
*   This function assumes ptwXY_integrateDomain checks status and coalesces the points.
*/
 
    int64_t i;
    nfu_status status; 
    double sum = ptwXY_integrateDomain( ptwXY, &status );
 
    if( status != nfu_Okay ) return( status );
    if( sum == 0. ) {
        status = nfu_badNorm; }
    else {
        for( i = 0; i < ptwXY->length; i++ ) ptwXY->points[i].y /= sum;
    }
    return( status );
}

ptwXY_pow()

nfu_status ptwXY_pow	(	ptwXYPoints *	ptwXY,
		double	p )

Definition at line 24 of file ptwXY_functions.cc.

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

ptwXY_reallocateOverflowPoints()

nfu_status ptwXY_reallocateOverflowPoints	(	ptwXYPoints *	ptwXY,
		int64_t	size )

Definition at line 439 of file ptwXY_core.cc.

                                                                              {
/*
*   This is for allocating/reallocating the secondary data memory.
*/
    nfu_status status = nfu_Okay;
 
    if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
 
    if( size < ptwXY_minimumOverflowSize ) size = ptwXY_minimumOverflowSize;      /* ptwXY_minimumOverflowSize must be > 0. */
    if( size < ptwXY->overflowLength ) status = ptwXY_coalescePoints( ptwXY, ptwXY->length + ptwXY->overflowAllocatedSize, NULL, 0 );
    if( status == nfu_Okay ) {
        if( size != ptwXY->overflowAllocatedSize ) {
            ptwXY->overflowPoints = (ptwXYOverflowPoint *) nfu_realloc( (size_t) size * sizeof( ptwXYOverflowPoint ), ptwXY->overflowPoints );
            if( ptwXY->overflowPoints == NULL ) {
                ptwXY->length = 0;
                ptwXY->overflowLength = 0;
                ptwXY->mallocFailedSize = size;
                size = 0;
                ptwXY->status = nfu_mallocError;
            }
        }
        ptwXY->overflowAllocatedSize = size; }
    else {
        ptwXY->status = status;
    }
    return( ptwXY->status );
}

Referenced by ptwXY_setup().

ptwXY_reallocatePoints()

nfu_status ptwXY_reallocatePoints	(	ptwXYPoints *	ptwXY,
		int64_t	size,
		int	forceSmallerResize )

Definition at line 410 of file ptwXY_core.cc.

                                                                                              {
/*
*   This is for allocating/reallocating the primary data memory.
*/
    if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
 
    if( size < ptwXY_minimumSize ) size = ptwXY_minimumSize;                      /* ptwXY_minimumSize must be > 0. */
    if( size < ptwXY->length ) size = ptwXY->length;
    if( size != ptwXY->allocatedSize ) {
        if( size > ptwXY->allocatedSize ) {                                        /* Increase size of allocated points. */
            ptwXY->points = (ptwXYPoint *) nfu_realloc( (size_t) size * sizeof( ptwXYPoint ), ptwXY->points ); }
        else if( ( ptwXY->allocatedSize > 2 * size ) || forceSmallerResize ) {     /* Decrease size, if at least 1/2 size reduction or if forced to. */
            ptwXY->points = (ptwXYPoint *) nfu_realloc( (size_t) size * sizeof( ptwXYPoint ), ptwXY->points ); }
        else {
            size = ptwXY->allocatedSize;                                           /* Size is < ptwXY->allocatedSize, but realloc not called. */
        }
        if( ptwXY->points == NULL ) {
            ptwXY->length = 0;
            ptwXY->mallocFailedSize = size;
            size = 0;
            ptwXY->status = nfu_mallocError;
        }
        ptwXY->allocatedSize = size;
    }
    return( ptwXY->status );
}

Referenced by ptwXY_coalescePoints(), ptwXY_copy(), ptwXY_setup(), ptwXY_setXYData(), and ptwXY_setXYDataFromXsAndYs().

ptwXY_release()

nfu_status ptwXY_release

(

ptwXYPoints *

ptwXY

)

Definition at line 549 of file ptwXY_core.cc.

                                               {
/*
*   Note, this routine does not free ptwXY (i.e., it does not undo all of ptwXY_new).
*/
 
    if( ptwXY->interpolation == ptwXY_interpolationOther ) {
        if( ptwXY->interpolationOtherInfo.interpolationString != NULL ) 
            ptwXY->interpolationOtherInfo.interpolationString = (char const *) nfu_free( (void *) ptwXY->interpolationOtherInfo.interpolationString );
    }
    ptwXY->interpolation = ptwXY_interpolationLinLin;
    ptwXY->interpolationOtherInfo.getValueFunc = NULL;
    ptwXY->interpolationOtherInfo.argList = NULL;
    ptwXY->length = 0;
    ptwXY->allocatedSize = 0;
    ptwXY->points = (ptwXYPoint *) nfu_free( ptwXY->points );
 
    ptwXY->overflowLength = 0;
    ptwXY->overflowAllocatedSize = 0;
    ptwXY->overflowPoints = (ptwXYOverflowPoint *) nfu_free( ptwXY->overflowPoints );
 
    return( nfu_Okay );
}

Referenced by ptwXY_free(), and ptwXY_setup().

ptwXY_runningIntegral()

ptwXPoints * ptwXY_runningIntegral	(	ptwXYPoints *	ptwXY,
		nfu_status *	status )

Definition at line 632 of file ptwXY_integration.cc.

                                                                            {
 
    int i;
    ptwXPoints *runningIntegral = NULL;
    double integral = 0., sum;
 
    if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( NULL );
    if( ( runningIntegral = ptwX_new( ptwXY->length, status ) ) == NULL ) goto err;
 
    if( ( *status = ptwX_setPointAtIndex( runningIntegral, 0, 0. ) ) != nfu_Okay ) goto err;
    for( i = 1; i < ptwXY->length; i++ ) {
        if( ( *status = ptwXY_f_integrate( ptwXY->interpolation, ptwXY->points[i-1].x, ptwXY->points[i-1].y, 
            ptwXY->points[i].x, ptwXY->points[i].y, &sum ) ) != nfu_Okay ) goto err;
        integral += sum;
        if( ( *status = ptwX_setPointAtIndex( runningIntegral, i, integral ) ) != nfu_Okay ) goto err;
    }
    return( runningIntegral );
 
err:
    if( runningIntegral != NULL ) ptwX_free( runningIntegral );
    return( NULL );
}

Referenced by MCGIDI_angular_parseFromTOM(), and MCGIDI_fromTOM_pdfOfX().

ptwXY_scaleOffsetXAndY()

nfu_status ptwXY_scaleOffsetXAndY	(	ptwXYPoints *	ptwXY,
		double	xScale,
		double	xOffset,
		double	yScale,
		double	yOffset )

Definition at line 481 of file ptwXY_methods.cc.

                                                                                                                      {
 
    int64_t i1, length = ptwXY->length;
    ptwXYPoint *p1;
    nfu_status status;
 
    if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
    if( xScale == 0 ) return( nfu_XNotAscending );
 
    if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status );
 
    for( i1 = 0, p1 = ptwXY->points; i1 < length; i1++, p1++ ) {
        p1->x = xScale * p1->x + xOffset;
        p1->y = yScale * p1->y + yOffset;
    }
 
    if( xScale < 0 ) {
        int64_t length_2 = length / 2;
        ptwXYPoint tmp, *p2;
 
        for( i1 = 0, p1 = ptwXY->points, p2 = &(ptwXY->points[length-1]); i1 < length_2; i1++ ) {
            tmp = *p1;
            *p1 = *p2;
            *p2 = tmp;
        }
    }
 
    return( ptwXY->status );
}

ptwXY_setAccuracy()

double ptwXY_setAccuracy	(	ptwXYPoints *	ptwXY,
		double	accuracy )

Definition at line 379 of file ptwXY_core.cc.

                                                                {
 
    if( accuracy < ptwXY_minAccuracy ) accuracy = ptwXY_minAccuracy;
    if( accuracy < ptwXY->accuracy ) accuracy = ptwXY->accuracy;
    if( accuracy > 1 ) accuracy = 1.;
    ptwXY->accuracy = accuracy;
    return( ptwXY->accuracy );
}

Referenced by ptwXY_setup().

ptwXY_setBiSectionMax()

double ptwXY_setBiSectionMax	(	ptwXYPoints *	ptwXY,
		double	biSectionMax )

Definition at line 397 of file ptwXY_core.cc.

                                                                        {
 
    if( biSectionMax < 0 ) {
        biSectionMax = 0; }
    else if( biSectionMax > ptwXY_maxBiSectionMax ) {
        biSectionMax = ptwXY_maxBiSectionMax;
    }
    ptwXY->biSectionMax = biSectionMax;
    return( ptwXY->biSectionMax );
}

Referenced by ptwXY_setup().

ptwXY_setup()

nfu_status ptwXY_setup	(	ptwXYPoints *	ptwXY,
		ptwXY_interpolation	interpolation,
		ptwXY_interpolationOtherInfo const *	interpolationOtherInfo,
		double	biSectionMax,
		double	accuracy,
		int64_t	primarySize,
		int64_t	secondarySize,
		int	userFlag )

Definition at line 46 of file ptwXY_core.cc.

                                                                                                     {
 
    ptwXY->status = nfu_Okay;
    ptwXY->typeX =  ptwXY_sigma_none;
    ptwXY->typeY =  ptwXY_sigma_none;
    ptwXY->interpolation = interpolation;
    ptwXY->interpolationOtherInfo.interpolationString = NULL;
    ptwXY->interpolationOtherInfo.getValueFunc = NULL;
    ptwXY->interpolationOtherInfo.argList = NULL;
    switch( interpolation ) {
    case ptwXY_interpolationLinLin :
        ptwXY->interpolationOtherInfo.interpolationString = linLinInterpolationString; break;
    case ptwXY_interpolationLinLog :
        ptwXY->interpolationOtherInfo.interpolationString = linLogInterpolationString; break;
    case ptwXY_interpolationLogLin :
        ptwXY->interpolationOtherInfo.interpolationString = logLinInterpolationString; break;
    case ptwXY_interpolationLogLog :
        ptwXY->interpolationOtherInfo.interpolationString = logLogInterpolationString; break;
    case ptwXY_interpolationFlat :
        ptwXY->interpolationOtherInfo.interpolationString = flatInterpolationString; break;
    case ptwXY_interpolationOther :         /* For ptwXY_interpolationOther, interpolationOtherInfo and interpolationString must be defined. */
        if( interpolationOtherInfo == NULL ) {              
                ptwXY->status = nfu_otherInterpolation; }
        else {
            if( interpolationOtherInfo->interpolationString == NULL ) {
                ptwXY->status = nfu_otherInterpolation; }
            else {
                if( ( ptwXY->interpolationOtherInfo.interpolationString = strdup( interpolationOtherInfo->interpolationString ) ) == NULL ) {
                    ptwXY->status = nfu_mallocError;
                }
            }
            ptwXY->interpolationOtherInfo.getValueFunc = interpolationOtherInfo->getValueFunc;
            ptwXY->interpolationOtherInfo.argList = interpolationOtherInfo->argList;
        }
    }
    ptwXY->userFlag = 0;
    ptwXY_setUserFlag( ptwXY, userFlag );
    ptwXY->biSectionMax = ptwXY_maxBiSectionMax;
    ptwXY_setBiSectionMax( ptwXY, biSectionMax );
    ptwXY->accuracy = ptwXY_minAccuracy;
    ptwXY_setAccuracy( ptwXY, accuracy );
 
    ptwXY->length = 0;
    ptwXY->allocatedSize = 0;
    ptwXY->overflowLength = 0;
    ptwXY->overflowAllocatedSize = 0;
    ptwXY->mallocFailedSize = 0;
 
    ptwXY_initialOverflowPoint( &(ptwXY->overflowHeader), &(ptwXY->overflowHeader), &(ptwXY->overflowHeader) );
 
    ptwXY->points = NULL;
    ptwXY->overflowPoints = NULL;
 
    ptwXY_reallocatePoints( ptwXY, primarySize, 0 );
    ptwXY_reallocateOverflowPoints( ptwXY, secondarySize );
    if( ptwXY->status != nfu_Okay ) ptwXY_release( ptwXY );
    return( ptwXY->status );
}

Referenced by ptwXY_new().

ptwXY_setUserFlag()

void ptwXY_setUserFlag	(	ptwXYPoints *	ptwXY,
		int	userFlag )

Definition at line 365 of file ptwXY_core.cc.

                                                           {
 
    ptwXY->userFlag = userFlag;
}

Referenced by ptwXY_setup().

ptwXY_setValueAtX()

nfu_status ptwXY_setValueAtX	(	ptwXYPoints *	ptwXY,
		double	x,
		double	y )

Definition at line 876 of file ptwXY_core.cc.

                                                                       {
 
    return( ptwXY_setValueAtX_overrideIfClose( ptwXY, x, y, 0., 0 ) );
}

Referenced by MCGIDI_angular_parseFromTOM(), MCGIDI_fromTOM_pdfsOfXGivenW(), ptwXY_clip(), ptwXY_convolution(), ptwXY_createGaussianCenteredSigma1(), ptwXY_div_ptwXY(), ptwXY_dullEdges(), ptwXY_flatInterpolationToLinear(), ptwXY_intersectionWith_ptwX(), ptwXY_mul2_ptwXY(), ptwXY_thicken(), ptwXY_valueTo_ptwXY(), and ptwXY_xSlice().

ptwXY_setValueAtX_overrideIfClose()

nfu_status ptwXY_setValueAtX_overrideIfClose	(	ptwXYPoints *	ptwXY,
		double	x,
		double	y,
		double	eps,
		int	override )

Definition at line 883 of file ptwXY_core.cc.

                                                                                                                 {
 
    int closeIsEqual;
    int64_t nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY ), i;
    nfu_status status = nfu_Okay;
    ptwXY_lessEqualGreaterX legx;
    ptwXYPoint *point = NULL, newPoint = { x, y };
    ptwXYOverflowPoint *overflowPoint, *p, *overflowHeader = &(ptwXY->overflowHeader);
    ptwXYOverflowPoint lessThanEqualXPoint, greaterThanXPoint;
    ptwXYPoint *closePoint = NULL;
 
    if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
 
    legx = ptwXY_getPointsAroundX_closeIsEqual( ptwXY, x, &lessThanEqualXPoint, &greaterThanXPoint, eps, &closeIsEqual, &closePoint );
    switch( legx ) {
    case ptwXY_lessEqualGreaterX_lessThan :
    case ptwXY_lessEqualGreaterX_greater :
    case ptwXY_lessEqualGreaterX_between :
        if( closeIsEqual ) {
            if( !override ) return( status );
            point = closePoint;
            legx = ptwXY_lessEqualGreaterX_equal;
            x = point->x; }
        else {
            if( ( legx == ptwXY_lessEqualGreaterX_greater ) && ( nonOverflowLength < ptwXY->allocatedSize ) ) {
                point = &(ptwXY->points[nonOverflowLength]); }
            else {
                if( ptwXY->overflowLength == ptwXY->overflowAllocatedSize ) 
                    return( ptwXY_coalescePoints( ptwXY, ptwXY->length + ptwXY->overflowAllocatedSize, &newPoint, 0 ) );
                overflowPoint = &(ptwXY->overflowPoints[ptwXY->overflowLength]);
                if( legx == ptwXY_lessEqualGreaterX_lessThan ) {
                    overflowPoint->prior = greaterThanXPoint.prior;
                    overflowPoint->index = 0; }
                else {                                      /* Between or greater and must go in overflow area. */
                    if( legx == ptwXY_lessEqualGreaterX_greater ) {
                        overflowPoint->prior = overflowHeader->prior;
                        overflowPoint->index = ptwXY->length; }
                    else {
                        overflowPoint->prior = lessThanEqualXPoint.prior;
                        if( lessThanEqualXPoint.next != NULL ) {
                            if( lessThanEqualXPoint.point.x < x ) 
                                overflowPoint->prior = lessThanEqualXPoint.prior->next;
                            i = 1; }
                        else {
                            for( p = overflowHeader->next, i = 1; p != overflowHeader; p = p->next, i++ ) 
                                if( p->point.x > x ) break;
                        }
                        overflowPoint->index = lessThanEqualXPoint.index + i;
                    }
                }
                overflowPoint->next = overflowPoint->prior->next;
                overflowPoint->prior->next = overflowPoint;
                overflowPoint->next->prior = overflowPoint;
                point = &(overflowPoint->point);
                for( overflowPoint = overflowPoint->next; overflowPoint != overflowHeader; overflowPoint = overflowPoint->next ) {
                    overflowPoint->index++;
                }
                ptwXY->overflowLength++;
            }
        }
        break;
    case ptwXY_lessEqualGreaterX_empty :
        point = ptwXY->points;                 /* ptwXY_minimumSize must be > 0 so there is always space here. */
        break;
    case ptwXY_lessEqualGreaterX_equal :
        if( closeIsEqual && !override ) return( status );
        if( lessThanEqualXPoint.next == NULL ) {
            point = &(ptwXY->points[lessThanEqualXPoint.index]); }
        else {
            point = &(lessThanEqualXPoint.prior->next->point);
        }
        break;
    }
    if( status == nfu_Okay ) {
        point->x = x;
        point->y = y;
        if( legx != ptwXY_lessEqualGreaterX_equal ) ptwXY->length++;
    }
    return( status );
}

Referenced by ptwXY_createFromFunction(), and ptwXY_setValueAtX().

ptwXY_setXYData()

nfu_status ptwXY_setXYData	(	ptwXYPoints *	ptwXY,
		int64_t	length,
		double const *	xy )

Definition at line 597 of file ptwXY_core.cc.

                                                                                   {
 
    nfu_status status = nfu_Okay;
    int64_t i;
    ptwXYPoint *p;
    double const *d = xy;
    double xOld = 0.;
 
    if( length > ptwXY->allocatedSize ) {
        status = ptwXY_reallocatePoints( ptwXY, length, 0 );
        if( status != nfu_Okay ) return( status );
    }
    for( i = 0, p = ptwXY->points; i < length; i++, p++ ) {
        if( i != 0 ) {
            if( *d <= xOld ) {
                status = nfu_XNotAscending;
                ptwXY->status = nfu_XNotAscending;
                length = 0;
                break;
            }
        }
        xOld = *d;
        p->x = *(d++);
        p->y = *(d++);
    }
    ptwXY->overflowHeader.next = &(ptwXY->overflowHeader);
    ptwXY->overflowHeader.prior = &(ptwXY->overflowHeader);
    ptwXY->overflowLength = 0;
    ptwXY->length = length;
    return( ptwXY->status = status );
}

Referenced by ptwXY_create().

ptwXY_setXYDataFromXsAndYs()

nfu_status ptwXY_setXYDataFromXsAndYs	(	ptwXYPoints *	ptwXY,
		int64_t	length,
		double const *	x,
		double const *	y )

Definition at line 631 of file ptwXY_core.cc.

                                                                                                              {
 
    nfu_status status;
    int64_t i;
    ptwXYPoint *p;
    double xOld = 0.;
 
    if( ( status = ptwXY_clear( ptwXY ) ) != nfu_Okay ) return( status );
    if( length > ptwXY->allocatedSize ) {
        if( ( status = ptwXY_reallocatePoints( ptwXY, length, 0 ) ) != nfu_Okay ) return( status );
    }
    for( i = 0, p = ptwXY->points; i < length; i++, p++, x++, y++ ) {
        if( i != 0 ) {
            if( *x <= xOld ) {
                status = ptwXY->status = nfu_XNotAscending;
                length = 0;
                break;
            }
        }
        xOld = *x;
        p->x = *x;
        p->y = *y;
    }
    ptwXY->length = length;
    return( status );
}

ptwXY_setXYPairAtIndex()

nfu_status ptwXY_setXYPairAtIndex	(	ptwXYPoints *	ptwXY,
		int64_t	index,
		double	x,
		double	y )

Definition at line 1098 of file ptwXY_core.cc.

                                                                                           {
 
    int64_t i, ip1;
    ptwXYOverflowPoint *overflowPoint, *pm1, *pp1;
 
    if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
 
    if( ( index < 0 ) || ( index >= ptwXY->length ) ) return( nfu_badIndex );
    for( overflowPoint = ptwXY->overflowHeader.next, i = 0; overflowPoint != &(ptwXY->overflowHeader); overflowPoint = overflowPoint->next, i++ ) {
        if( overflowPoint->index >= index ) break;
    }
    ip1 = i;
    pm1 = pp1 = overflowPoint;
    if( overflowPoint->index == index ) {                                           /* Note, if overflowPoint is header, then its index = -1. */
        pp1 = overflowPoint->next;
        ip1++;
    }
    if( ( pp1 != &(ptwXY->overflowHeader) ) && ( pp1->index == ( index + 1 ) ) ) {     /* This if and else check that x < element[index+1]'s x values. */
        if( pp1->point.x <= x ) return( nfu_badIndexForX ); }
    else {
        if( ( ( index + 1 ) < ptwXY->length ) && ( ptwXY->points[index + 1 - ip1].x <= x ) ) return( nfu_badIndexForX );
    }
    if( overflowPoint != &(ptwXY->overflowHeader) ) pm1 = overflowPoint->prior;
    if( ( pm1 != &(ptwXY->overflowHeader) ) && ( pm1->index == ( index - 1 ) ) ) {     /* This if and else check that x > element[index-1]'s x values. */
        if( pm1->point.x >= x ) return( nfu_badIndexForX ); }
    else {
        if( ( ( index - 1 ) >= 0 ) && ( ptwXY->points[index - 1 - i].x >= x ) ) return( nfu_badIndexForX );
    }
    if( ( overflowPoint != &(ptwXY->overflowHeader) ) && ( overflowPoint->index == index ) ) {
        overflowPoint->point.x = x;
        overflowPoint->point.y = y; }
    else {
        index -= i;
        ptwXY->points[index].x = x;
        ptwXY->points[index].y = y;
    }
    return( nfu_Okay );
}

ptwXY_showInteralStructure()

void ptwXY_showInteralStructure	(	ptwXYPoints *	ptwXY,
		FILE *	f,
		int	printPointersAsNull )

Definition at line 253 of file ptwXY_misc.cc.

                                                                                        {
 
    int64_t i, n = ptwXY_getNonOverflowLength( ptwXY );
    ptwXYPoint *point = ptwXY->points;
    ptwXYOverflowPoint *overflowPoint;
 
    fprintf( f, "status = %d  interpolation = %d  length = %d  allocatedSize = %d\n", 
        (int) ptwXY->status, (int) ptwXY->interpolation, (int) ptwXY->length, (int) ptwXY->allocatedSize );
    fprintf( f, "userFlag = %d  biSectionMax = %.8e  accuracy = %.2e  minFractional_dx = %.6e\n", 
        ptwXY->userFlag, ptwXY->biSectionMax, ptwXY->accuracy, ptwXY->minFractional_dx );
    fprintf( f, "interpolationString = %s\n", ptwXY->interpolationOtherInfo.interpolationString );
    fprintf( f, "getValueFunc is NULL = %d. argList is NULL = %d.\n", 
        ptwXY->interpolationOtherInfo.getValueFunc == NULL, ptwXY->interpolationOtherInfo.argList == NULL );
    fprintf( f, "  overflowLength = %d  overflowAllocatedSize = %d  mallocFailedSize = %d\n", 
        (int) ptwXY->overflowLength, (int) ptwXY->overflowAllocatedSize, (int) ptwXY->mallocFailedSize );
    fprintf( f, "  Points data, points = %20p\n", ( printPointersAsNull ? NULL : (void*)ptwXY->points ) );
    for( i = 0; i < n; i++,  point++ ) fprintf( f, "    %14.7e %14.7e\n", point->x, point->y );
    fprintf( f, "  Overflow points data; %20p\n", ( printPointersAsNull ? NULL : (void*)&(ptwXY->overflowHeader) ) );
    for( overflowPoint = ptwXY->overflowHeader.next; overflowPoint != &(ptwXY->overflowHeader); overflowPoint = overflowPoint->next ) {
        fprintf( f, "    %14.7e %14.7e %8d %20p %20p %20p\n", overflowPoint->point.x, overflowPoint->point.y, (int) overflowPoint->index, 
    (void*) ( printPointersAsNull ? NULL : overflowPoint ), (void*) ( printPointersAsNull ? NULL : overflowPoint->prior ), 
    (void*) ( printPointersAsNull ? NULL : overflowPoint->next ) );
    }
    fprintf( f, "  Points in order\n" );
    for( i = 0; i < ptwXY->length; i++ ) {
        point = ptwXY_getPointAtIndex( ptwXY, i );
        fprintf( f, "    %14.7e %14.7e\n", point->x, point->y );
    }
}

ptwXY_simpleCoalescePoints()

nfu_status ptwXY_simpleCoalescePoints

(

ptwXYPoints *

ptwXY

)

Definition at line 529 of file ptwXY_core.cc.

                                                            {
 
    return( ptwXY_coalescePoints( ptwXY, ptwXY->length, NULL, 0 ) );
}

Referenced by MCGIDI_angular_parseFromTOM(), MCGIDI_fromTOM_pdfOfX(), MCGIDI_product_parseFromTOM(), ptwXY_applyFunction(), ptwXY_clip(), ptwXY_convolution(), ptwXY_copyToC_XY(), ptwXY_createFromFunction(), ptwXY_deletePoints(), ptwXY_div_ptwXY(), ptwXY_exp(), ptwXY_flatInterpolationToLinear(), ptwXY_getXArray(), ptwXY_groupOneFunction(), ptwXY_groupThreeFunctions(), ptwXY_groupTwoFunctions(), ptwXY_integrate(), ptwXY_integrateWithFunction(), ptwXY_integrateWithWeight_sqrt_x(), ptwXY_integrateWithWeight_x(), ptwXY_intersectionWith_ptwX(), ptwXY_mergeClosePoints(), ptwXY_mul2_ptwXY(), ptwXY_runningIntegral(), ptwXY_scaleOffsetXAndY(), ptwXY_slice(), ptwXY_thicken(), ptwXY_thin(), ptwXY_trim(), ptwXY_union(), and ptwXY_valueTo_ptwXAndY().

ptwXY_simplePrint()

void ptwXY_simplePrint	(	ptwXYPoints *	ptwXY,
		char *	format )

Definition at line 298 of file ptwXY_misc.cc.

                                                           {
 
    ptwXY_simpleWrite( ptwXY, stdout, format );
}

ptwXY_simpleWrite()

void ptwXY_simpleWrite	(	ptwXYPoints *	ptwXY,
		FILE *	f,
		char *	format )

Definition at line 285 of file ptwXY_misc.cc.

                                                                    {
 
    int64_t i;
    ptwXYPoint *point;
 
    for( i = 0; i < ptwXY->length; i++ ) {
        point = ptwXY_getPointAtIndex( ptwXY, i );
        fprintf( f, format, point->x, point->y );
    }
}

Referenced by ptwXY_simplePrint().

ptwXY_slice()

ptwXYPoints * ptwXY_slice	(	ptwXYPoints *	ptwXY,
		int64_t	index1,
		int64_t	index2,
		int64_t	secondarySize,
		nfu_status *	status )

Definition at line 248 of file ptwXY_core.cc.

                                                                                                                          {
 
    int64_t i, length;
    ptwXYPoints *n;
 
    *status = nfu_badSelf;
    if( ptwXY->status != nfu_Okay ) return( NULL );
 
    *status = nfu_badIndex;
    if( index2 < index1 ) return( NULL );
    if( index1 < 0 ) index1 = 0;
    if( index2 > ptwXY->length ) index2 = ptwXY->length;
 
    length = index2 - index1;
    if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( NULL );
    if( ( n = ptwXY_new( ptwXY->interpolation, &(ptwXY->interpolationOtherInfo), ptwXY->biSectionMax, 
        ptwXY->accuracy, length, secondarySize, status, ptwXY->userFlag ) ) == NULL ) return( NULL );
 
    *status = n->status = ptwXY->status;
    for( i = index1; i < index2; i++ ) n->points[i - index1] = ptwXY->points[i];
    n->length = length;
    return( n );
}

Referenced by ptwXY_clone().

ptwXY_slopeOffset()

nfu_status ptwXY_slopeOffset	(	ptwXYPoints *	ptwXY,
		double	slope,
		double	offset )

Definition at line 25 of file ptwXY_binaryOperators.cc.

                                                                                { 
 
    int64_t i, nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
    ptwXYPoint *p;
    ptwXYOverflowPoint *o, *overflowHeader = &(ptwXY->overflowHeader);
 
    if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
 
    for( i = 0, p = ptwXY->points; i < nonOverflowLength; i++, p++ ) p->y = slope * p->y + offset;
    for( o = overflowHeader->next; o != overflowHeader; o = o->next ) o->point.y = slope * o->point.y + offset; 
    return( ptwXY->status );
}   

Referenced by ptwXY_add_double(), ptwXY_div_doubleFrom(), ptwXY_mul_double(), ptwXY_sub_doubleFrom(), and ptwXY_sub_fromDouble().

ptwXY_sub_doubleFrom()

nfu_status ptwXY_sub_doubleFrom	(	ptwXYPoints *	ptwXY,
		double	value )

Definition at line 41 of file ptwXY_binaryOperators.cc.

{ return( ptwXY_slopeOffset( ptwXY,  1., -value ) ); }

ptwXY_sub_fromDouble()

nfu_status ptwXY_sub_fromDouble	(	ptwXYPoints *	ptwXY,
		double	value )

Definition at line 42 of file ptwXY_binaryOperators.cc.

{ return( ptwXY_slopeOffset( ptwXY, -1.,  value ) ); }

ptwXY_sub_ptwXY()

ptwXYPoints * ptwXY_sub_ptwXY	(	ptwXYPoints *	ptwXY1,
		ptwXYPoints *	ptwXY2,
		nfu_status *	status )

Definition at line 154 of file ptwXY_binaryOperators.cc.

                                                                                             {
 
    ptwXYPoints *diff;
 
    if( ptwXY1->length == 0 ) {
        diff = ptwXY_clone( ptwXY2, status );
        if( ( *status = ptwXY_neg( diff ) ) != nfu_Okay ) diff = ptwXY_free( diff ); }
    else if( ptwXY2->length == 0 ) {
        diff = ptwXY_clone( ptwXY1, status ); }
    else {
        diff = ptwXY_binary_ptwXY( ptwXY1, ptwXY2, 1., -1., 0., status );
    }
    return( diff );
}

ptwXY_thicken()

nfu_status ptwXY_thicken	(	ptwXYPoints *	ptwXY1,
		int	sectionSubdivideMax,
		double	dxMax,
		double	fxMax )

Definition at line 144 of file ptwXY_methods.cc.

                                                                                                     {
 
    double x1, x2 = 0., y1, y2 = 0., fx = 1.1, x, dx, dxp, lfx, y;    /* fx initialized so compilers want complain. */
    int64_t i, notFirstPass = 0;
    int nfx, nDone, doLinear;
    nfu_status status;
 
    if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
    if( ( sectionSubdivideMax < 1 ) || ( dxMax < 0. ) || ( fxMax < 1. ) ) return( nfu_badInput );
    if( sectionSubdivideMax > ptwXY_sectionSubdivideMax ) sectionSubdivideMax = ptwXY_sectionSubdivideMax;
    if( ( status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( status );
    for( i = ptwXY1->length - 1; i >= 0; i-- ) {
        x1 = ptwXY1->points[i].x;
        y1 = ptwXY1->points[i].y;
        if( notFirstPass ) {
            dx = ptwXY_thicken_linear_dx( sectionSubdivideMax, dxMax, x1, x2 );
 
            if( x1 == 0. ) {
                doLinear = 1; }
            else {
                fx = x2 / x1;
                if( fx > 0. ) {
                    lfx = G4Log( fx );
                    if( fxMax == 1. ) {
                        nfx = sectionSubdivideMax; }
                    else {
                        nfx = ( (int) ( lfx / G4Log( fxMax ) ) ) + 1;
                        if( nfx > sectionSubdivideMax ) nfx = sectionSubdivideMax;
                    }
                    if( nfx > 0 ) fx = G4Exp( lfx / nfx );
                    doLinear = 0;
                    if( dx < ( fx - 1 ) * x1 ) doLinear = 1; }
                else {
                    doLinear = 1;
                }
            }
            x = x1;
            dxp = dx;
            nDone = 0;
            while( 1 ) {
                if( doLinear ) {
                    x += dx; }
                else {
                    dx = ptwXY_thicken_linear_dx( sectionSubdivideMax - nDone, dxMax, x, x2 );
                    if( dx <= ( fx - 1 ) * x ) {
                        dxp = dx;
                        doLinear = 1;
                        continue;
                    }
                    dxp = ( fx - 1. ) * x;
                    x *= fx;
                }
                if( ( x2 - x ) < 0.05 * std::fabs( dxp ) ) break;
                if( ( status = ptwXY_interpolatePoint( ptwXY1->interpolation, x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( status );
                if( ( status = ptwXY_setValueAtX( ptwXY1, x, y ) ) != nfu_Okay ) return( status );
                nDone++;
            } // Loop checking, 11.06.2015, T. Koi
        } 
        notFirstPass = 1;
        x2 = x1;
        y2 = y1;
    }
    return( status );
}

ptwXY_thin()

ptwXYPoints * ptwXY_thin	(	ptwXYPoints *	ptwXY1,
		double	accuracy,
		nfu_status *	status )

Definition at line 211 of file ptwXY_methods.cc.

                                                                                    {
 
    int64_t i, j, length = ptwXY1->length;
    ptwXYPoints *thinned = NULL;
    double y1, y2, y3;
    char *thin = NULL;
 
    if( length < 3 ) return( ptwXY_clone( ptwXY1, status ) );   /* Logic below requires at least 2 points. */
    if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
    *status = nfu_otherInterpolation;
    if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
    if( accuracy < ptwXY1->accuracy ) accuracy = ptwXY1->accuracy;
    if( ( thinned = ptwXY_new( ptwXY1->interpolation, &(ptwXY1->interpolationOtherInfo), 
        ptwXY1->biSectionMax, accuracy, length, ptwXY1->overflowLength, status, ptwXY1->userFlag ) ) == NULL ) return( NULL );
 
    thinned->points[0] = ptwXY1->points[0];                     /* This sections removes middle point if surrounding points have the same y-value. */
    y1 = ptwXY1->points[0].y;
    y2 = ptwXY1->points[1].y;
    for( i = 2, j = 1; i < length; i++ ) {
        y3 = ptwXY1->points[i].y;
        if( ( y1 != y2 ) || ( y2 != y3 ) ) {
            thinned->points[j++] = ptwXY1->points[i - 1];
            y1 = y2;
            y2 = y3;
        }
    }
    thinned->points[j++] = ptwXY1->points[length - 1];
 
    if( ptwXY1->interpolation != ptwXY_interpolationFlat ) {    /* Now call ptwXY_thin2 for more thinning. */
        length = thinned->length = j;
        if( ( thin = (char *) nfu_calloc( 1, (size_t) length ) ) == NULL ) goto Err;
        if( ( *status = ptwXY_thin2( thinned, thin, accuracy, 0, length - 1 ) ) != nfu_Okay ) goto Err;
        for( j = 1; j < length; j++ ) if( thin[j] != 0 ) break;
        for( i = j + 1; i < length; i++ ) {
            if( thin[i] == 0 ) {
                thinned->points[j] = thinned->points[i];
                j++;
            }
        }
        nfu_free( thin );
    }
    thinned->length = j;
 
    return( thinned );
 
Err:
    ptwXY_free( thinned );
    if( thin != NULL ) nfu_free( thin );
    return( NULL );
}

ptwXY_toOtherInterpolation()

ptwXYPoints * ptwXY_toOtherInterpolation	(	ptwXYPoints *	ptwXY,
		ptwXY_interpolation	interpolation,
		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_trim()

nfu_status ptwXY_trim

(

ptwXYPoints *

ptwXY

)

Definition at line 315 of file ptwXY_methods.cc.

                                            {
/*
c   Remove extra zeros at beginning and end.
*/
 
    int64_t i, i1, i2;
    nfu_status status;
 
    if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
    if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status );
    for( i1 = 0; i1 < ptwXY->length; i1++ ) {
        if( ptwXY->points[i1].y != 0 ) break;
    }
    if( i1 > 0 ) i1--;
    for( i2 = ptwXY->length - 1; i2 >= 0; i2-- ) {
        if( ptwXY->points[i2].y != 0 ) break;
    }
    i2++;
    if( i2 < ptwXY->length ) i2++;
    if( i2 > i1 ) {
        if( i1 > 0 ) {
            for( i = i1; i < i2; i++ ) ptwXY->points[i - i1] = ptwXY->points[i];
        }
        ptwXY->length = i2 - i1; }
    else if( i2 < i1 ) {                /* Remove all zeros between endpoints. */
        ptwXY->points[1] = ptwXY->points[ptwXY->length - 1];
        ptwXY->length = 2;
    }
 
    return( nfu_Okay );
}

ptwXY_tweakDomainsToMutualify()

nfu_status ptwXY_tweakDomainsToMutualify	(	ptwXYPoints *	ptwXY1,
		ptwXYPoints *	ptwXY2,
		int	epsilonFactor,
		double	epsilon )

Definition at line 295 of file ptwXY_convenient.cc.

                                                                                                                        {
 
    nfu_status status;
    int64_t n1 = ptwXY1->length, n2 = ptwXY2->length;
    double sum, diff;
    ptwXYPoint *xy1, *xy2;
 
    epsilon = std::fabs( epsilon ) + std::fabs( epsilonFactor * DBL_EPSILON );
 
    if( ( status = ptwXY1->status ) != nfu_Okay ) return( status );
    if( ( status = ptwXY2->status ) != nfu_Okay ) return( status );
    if( n1 == 0 ) return( nfu_empty );
    if( n2 == 0 ) return( nfu_empty );
    if( n1 < 2 ) { 
        status = nfu_tooFewPoints; }
    else if( n2 < 2 ) { 
        status = nfu_tooFewPoints; }
    else {
        xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, 0 );
        xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, 0 );
        if( xy1->x < xy2->x ) {
            if( xy2->y != 0. ) {
                sum = std::fabs( xy1->x ) + std::fabs( xy2->x );
                diff = std::fabs( xy2->x - xy1->x );
                if( diff > epsilon * sum ) {
                    status = nfu_domainsNotMutual; }
                else {
                    xy1->x = xy2->x;
                }
            } }
        else if( xy1->x > xy2->x ) {
            if( xy1->y != 0. ) {
                sum = std::fabs( xy1->x ) + std::fabs( xy2->x );
                diff = std::fabs( xy2->x - xy1->x );
                if( diff > epsilon * sum ) {
                    status = nfu_domainsNotMutual; }
                else {
                    xy2->x = xy1->x;
                }
            }
        }
 
        if( status == nfu_Okay ) {
            xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, n1 - 1 );
            xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, n2 - 1 );
            if( xy1->x < xy2->x ) {
                if( xy1->y != 0. ) {
                    sum = std::fabs( xy1->x ) + std::fabs( xy2->x );
                    diff = std::fabs( xy2->x - xy1->x );
                    if( diff > epsilon * sum ) {
                        status = nfu_domainsNotMutual; }
                    else {
                        xy2->x = xy1->x;
                    }
                } }
            else if( xy1->x > xy2->x ) {
                if( xy2->y != 0. ) {
                    sum = std::fabs( xy1->x ) + std::fabs( xy2->x );
                    diff = std::fabs( xy2->x - xy1->x );
                    if( diff > epsilon * sum ) {
                        status = nfu_domainsNotMutual; }
                    else {
                        xy1->x = xy2->x;
                    }
                }
            }
        }
    }
    return( status );
}

Referenced by ptwXY_groupThreeFunctions(), and ptwXY_groupTwoFunctions().

ptwXY_union()

ptwXYPoints * ptwXY_union	(	ptwXYPoints *	ptwXY1,
		ptwXYPoints *	ptwXY2,
		nfu_status *	status,
		int	unionOptions )

Definition at line 349 of file ptwXY_methods.cc.

                                                                                                           {
 
    int64_t overflowSize, i, i1 = 0, i2 = 0, n1 = ptwXY1->length, n2 = ptwXY2->length, length;
    int fillWithFirst = unionOptions & ptwXY_union_fill, trim = unionOptions & ptwXY_union_trim;
    ptwXYPoints *n;
    double x1 = 0., x2 = 0., y1 = 0., y2 = 0., y, biSectionMax, accuracy;
 
    if( ( *status = ptwXY1->status ) != nfu_Okay ) return( NULL );
    if( ( *status = ptwXY2->status ) != nfu_Okay ) return( NULL );
    *status = nfu_otherInterpolation;
    if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
/*
*   Many other routines use the fact that ptwXY_union calls ptwXY_coalescePoints for ptwXY1 and ptwXY2 so do not change it.
*/
    if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
    if( ( *status = ptwXY_simpleCoalescePoints( ptwXY2 ) ) != nfu_Okay ) return( NULL );
 
    if( ( n1 == 1 ) || ( n2 == 1 ) ) {
        *status = nfu_tooFewPoints;
        return( NULL );
    }
    if( trim ) {
        if( n1 > 0 ) {
            if( n2 > 0 ) {
                if( ptwXY1->points[0].x < ptwXY2->points[0].x ) {
                    while( i1 < n1 ) { // Loop checking, 11.05.2015, T. Koi
                        if( ptwXY1->points[i1].x >= ptwXY2->points[0].x ) break;
                        if( fillWithFirst ) {
                            if( i1 < ( ptwXY1->length - 1 ) ) {
                                x1 = ptwXY1->points[i1].x;
                                y1 = ptwXY1->points[i1].y;
                                x2 = ptwXY1->points[i1+1].x;
                                y2 = ptwXY1->points[i1+1].y; 
                            }
                        }
                        i1++;
                    } }
                else {
                    while( i2 < n2 ) { // Loop checking, 11.06.2015, T. Koi
                        if( ptwXY2->points[i2].x >= ptwXY1->points[0].x ) break;
                        i2++;
                    }
                }
                if( ptwXY1->points[n1-1].x > ptwXY2->points[n2-1].x ) {
                    while( i1 < n1 ) { // Loop checking, 11.06.2015, T. Koi
                        if( ptwXY1->points[n1-1].x <= ptwXY2->points[n2-1].x ) break;
                        n1--;
                    } }
                else {
                    while( i2 < n2 ) { // Loop checking, 11.06.2015, T. Koi
                        if( ptwXY2->points[n2-1].x <= ptwXY1->points[n1-1].x ) break;
                        n2--;
                    }
                } }
            else {
                n1 = 0;
            } }
        else {
            n2 = 0;
        }
    }
    overflowSize = ptwXY1->overflowAllocatedSize;
    if( overflowSize < ptwXY2->overflowAllocatedSize ) overflowSize = ptwXY2->overflowAllocatedSize;
    length = ( n1 - i1 ) + ( n2 - i2 );
    if( length == 0 ) length = ptwXY_minimumSize;
    biSectionMax = ptwXY1->biSectionMax;
    if( biSectionMax < ptwXY2->biSectionMax ) biSectionMax = ptwXY2->biSectionMax;
    accuracy = ptwXY1->accuracy;
    if( accuracy < ptwXY2->accuracy ) accuracy = ptwXY2->accuracy;
    n = ptwXY_new( ptwXY1->interpolation, NULL, biSectionMax, accuracy, length, overflowSize, status, ptwXY1->userFlag );
    if( n == NULL ) return( NULL );
 
    for( i = 0; ( i1 < n1 ) && ( i2 < n2 ); i++ ) {
        y = 0.;
        if( ptwXY1->points[i1].x <= ptwXY2->points[i2].x ) {
            n->points[i].x = ptwXY1->points[i1].x;
            if( fillWithFirst ) {
                y = ptwXY1->points[i1].y;
                if( i1 < ( ptwXY1->length - 1 ) ) {
                    x1 = ptwXY1->points[i1].x;
                    y1 = ptwXY1->points[i1].y;
                    x2 = ptwXY1->points[i1+1].x;
                    y2 = ptwXY1->points[i1+1].y; }
                else {
                    y1 = 0.;
                    y2 = 0.;
                }
            }
            if( ptwXY1->points[i1].x == ptwXY2->points[i2].x ) i2++;
            i1++; }
        else {
            n->points[i].x = ptwXY2->points[i2].x;
            if( fillWithFirst && ( ( y1 != 0. ) || ( y2 != 0. ) ) ) {
                if( ( *status = ptwXY_interpolatePoint( ptwXY1->interpolation, ptwXY2->points[i2].x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) {
                    ptwXY_free( n );
                    return( NULL );
                }
            }
            i2++;
        }
        n->points[i].y = y;
    }
 
    y = 0.;
    for( ; i1 < n1; i1++, i++ ) {
        n->points[i].x = ptwXY1->points[i1].x;
        if( fillWithFirst ) y = ptwXY1->points[i1].y;
        n->points[i].y = y;
    }
    for( ; i2 < n2; i2++, i++ ) {
        n->points[i].x = ptwXY2->points[i2].x;
        if( fillWithFirst && trim && ( n->points[i].x <= x2 ) ) {
            if( ( *status = ptwXY_interpolatePoint( ptwXY1->interpolation, n->points[i].x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) {
                ptwXY_free( n );
                return( NULL );
            }
        }
        n->points[i].y = y;
    }
    n->length = i;
 
    if( unionOptions & ptwXY_union_mergeClosePoints ) {
        if( ( *status = ptwXY_mergeClosePoints( n, 4 * DBL_EPSILON ) ) != nfu_Okay ) {
            ptwXY_free( n );
            return( NULL );
        }
    }
    return( n );
}

Referenced by ptwXY_binary_ptwXY(), ptwXY_div_ptwXY(), ptwXY_groupThreeFunctions(), and ptwXY_groupTwoFunctions().

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 );
}

ptwXY_update_biSectionMax()

void ptwXY_update_biSectionMax	(	ptwXYPoints *	ptwXY1,
		double	oldLength )

Definition at line 31 of file ptwXY_misc.cc.

                                                                        {
 
    ptwXY1->biSectionMax = ptwXY1->biSectionMax - 1.442695 * G4Log( ptwXY1->length / oldLength ); /* 1.442695 = 1 / std::log( 2. ) */
    if( ptwXY1->biSectionMax < 0 ) ptwXY1->biSectionMax = 0;
    if( ptwXY1->biSectionMax > ptwXY_maxBiSectionMax ) ptwXY1->biSectionMax = ptwXY_maxBiSectionMax;
}

Referenced by ptwXY_applyFunction(), ptwXY_div_ptwXY(), and ptwXY_mul2_ptwXY().

ptwXY_valueTo_ptwXAndY()

nfu_status ptwXY_valueTo_ptwXAndY	(	ptwXYPoints *	ptwXY,
		double **	xs,
		double **	ys )

Definition at line 450 of file ptwXY_convenient.cc.

                                                                                  {
 
    int64_t i1, length = ptwXY_length( ptwXY );
    double *xps, *yps;
    ptwXYPoint *pointFrom;
    nfu_status status;
 
    if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
    if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status );
 
    if( ( *xs = (double *) malloc( length * sizeof( double ) ) ) == NULL ) return( nfu_mallocError );
    if( ( *ys = (double *) malloc( length * sizeof( double ) ) ) == NULL ) {
        free( *xs );
        *xs = NULL;
        return( nfu_mallocError );
    }
 
    for( i1 = 0, pointFrom = ptwXY->points, xps = *xs, yps = *ys; i1 < length; ++i1, ++pointFrom, ++xps, ++yps ) {
        *xps = pointFrom->x;
        *yps = pointFrom->y;
    }
 
    return( nfu_Okay );
}

ptwXY_valueTo_ptwXY()

ptwXYPoints * ptwXY_valueTo_ptwXY	(	double	x1,
		double	x2,
		double	y,
		nfu_status *	status )

Definition at line 477 of file ptwXY_convenient.cc.

                                                                                       {
 
    ptwXYPoints *n;
 
    *status = nfu_XNotAscending;
    if( x1 >= x2 ) return( NULL );
    *status = nfu_Okay;
    if( ( n = ptwXY_new( ptwXY_interpolationLinLin, NULL, ptwXY_maxBiSectionMax, ptwXY_minAccuracy, 2, 0, status, 0 ) ) == NULL ) return( NULL );
    ptwXY_setValueAtX( n, x1, y );
    ptwXY_setValueAtX( n, x2, y );
    return( n );
}

ptwXY_xMaxSlice()

ptwXYPoints * ptwXY_xMaxSlice	(	ptwXYPoints *	ptwXY,
		double	xMax,
		int64_t	secondarySize,
		int	fill,
		nfu_status *	status )

Definition at line 326 of file ptwXY_core.cc.

                                                                                                                     {
 
    double xMin = 0.9 * xMax - 1;
 
    if( xMax < 0 ) xMin = 1.1 * xMax - 1;
    if( ptwXY->length > 0 ) xMin = ptwXY_getXMin( ptwXY );
    return( ptwXY_xSlice( ptwXY, xMin, xMax, secondarySize, fill, status ) );
}

ptwXY_xMinSlice()

ptwXYPoints * ptwXY_xMinSlice	(	ptwXYPoints *	ptwXY,
		double	xMin,
		int64_t	secondarySize,
		int	fill,
		nfu_status *	status )

Definition at line 315 of file ptwXY_core.cc.

                                                                                                                     {
 
    double xMax = 1.1 * xMin + 1;
 
    if( xMin < 0 ) xMax = 0.9 * xMin + 1;
    if( ptwXY->length > 0 ) xMax = ptwXY_getXMax( ptwXY );
    return( ptwXY_xSlice( ptwXY, xMin, xMax, secondarySize, fill, status ) );
}

ptwXY_xSlice()

ptwXYPoints * ptwXY_xSlice	(	ptwXYPoints *	ptwXY,
		double	xMin,
		double	xMax,
		int64_t	secondarySize,
		int	fill,
		nfu_status *	status )

Definition at line 274 of file ptwXY_core.cc.

                                                                                                                               {
 
    int64_t i, i1, i2;
    double y;
    ptwXYPoints *n = NULL;
 
    if( ( *status = ptwXY->status ) != nfu_Okay ) return( NULL );
 
    if( ( ptwXY->length == 0 ) || ( ptwXY_getXMin( ptwXY ) >= xMax ) || ( ptwXY_getXMax( ptwXY ) <= xMin ) ) {
        n = ptwXY_new( ptwXY->interpolation, &(ptwXY->interpolationOtherInfo), ptwXY->biSectionMax, 
            ptwXY->accuracy, 0, secondarySize, status, ptwXY->userFlag ); }
    else {
        if( ( n = ptwXY_clone( ptwXY, status ) ) == NULL ) return( NULL );
        if( ( n->points[0].x < xMin ) || ( n->points[n->length - 1].x > xMax ) ) {
            if( fill && ( n->points[n->length - 1].x > xMax ) ) {
                if( ( *status = ptwXY_getValueAtX( n, xMax, &y ) ) != nfu_Okay ) goto Err;
                if( ( *status = ptwXY_setValueAtX( n, xMax,  y ) ) != nfu_Okay ) goto Err;
            }
            if( fill && ( n->points[0].x < xMin ) ) {
                if( ( *status = ptwXY_getValueAtX( n, xMin, &y ) ) != nfu_Okay ) goto Err;
                if( ( *status = ptwXY_setValueAtX( n, xMin,  y ) ) != nfu_Okay ) goto Err;
            }
            ptwXY_coalescePoints( n, n->length + n->overflowAllocatedSize, NULL, 0 );
            for( i1 = 0; i1 < n->length; i1++ ) if( n->points[i1].x >= xMin ) break;
            for( i2 = n->length - 1; i2 > 0; i2-- ) if( n->points[i2].x <= xMax ) break;
            i2++;
            if( i1 > 0 ) {
                for( i = i1; i < i2; i++ ) n->points[i- i1] = n->points[i];
            }
            n->length = i2 - i1;
        }
    }
    return( n );
 
Err:
    if( n != NULL ) ptwXY_free( n );
    return( NULL );
}

Referenced by GIDI_settings_processedFlux::groupFunction(), ptwXY_createGaussian(), ptwXY_xMaxSlice(), and ptwXY_xMinSlice().