Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ptwXY.h File Reference
#include <stdio.h>
#include <stdint.h>
#include <nf_utilities.h>
#include <ptwX.h>
Include dependency graph for ptwXY.h:
This graph shows which files directly or indirectly include this file:

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

ptwXYPointsptwXY_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)
 
ptwXYPointsptwXY_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)
 
ptwXYPointsptwXY_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)
 
ptwXYPointsptwXY_clone (ptwXYPoints *ptwXY, nfu_status *status)
 
ptwXYPointsptwXY_cloneToInterpolation (ptwXYPoints *ptwXY, ptwXY_interpolation interpolationTo, nfu_status *status)
 
ptwXYPointsptwXY_slice (ptwXYPoints *ptwXY, int64_t index1, int64_t index2, int64_t secondarySize, nfu_status *status)
 
ptwXYPointsptwXY_xSlice (ptwXYPoints *ptwXY, double xMin, double xMax, int64_t secondarySize, int fill, nfu_status *status)
 
ptwXYPointsptwXY_xMinSlice (ptwXYPoints *ptwXY, double xMin, int64_t secondarySize, int fill, nfu_status *status)
 
ptwXYPointsptwXY_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)
 
ptwXYPointsptwXY_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)
 
ptwXYPointptwXY_getPointAtIndex (ptwXYPoints *ptwXY, int64_t index)
 
ptwXYPointptwXY_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)
 
ptwXYPointsptwXY_thin (ptwXYPoints *ptwXY1, double accuracy, nfu_status *status)
 
nfu_status ptwXY_trim (ptwXYPoints *ptwXY)
 
ptwXYPointsptwXY_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)
 
ptwXYPointsptwXY_binary_ptwXY (ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, double v1, double v2, double v1v2, nfu_status *status)
 
ptwXYPointsptwXY_add_ptwXY (ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status)
 
ptwXYPointsptwXY_sub_ptwXY (ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status)
 
ptwXYPointsptwXY_mul_ptwXY (ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status)
 
ptwXYPointsptwXY_mul2_ptwXY (ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status)
 
ptwXYPointsptwXY_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)
 
ptwXYPointsptwXY_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)
 
ptwXYPointsptwXY_flatInterpolationToLinear (ptwXYPoints *ptwXY, double lowerEps, double upperEps, nfu_status *status)
 
ptwXYPointsptwXY_toOtherInterpolation (ptwXYPoints *ptwXY, ptwXY_interpolation interpolation, double accuracy, nfu_status *status)
 
ptwXYPointsptwXY_unitbaseInterpolate (double w, double w1, ptwXYPoints *ptwXY1, double w2, ptwXYPoints *ptwXY2, nfu_status *status)
 
ptwXYPointsptwXY_toUnitbase (ptwXYPoints *ptwXY, nfu_status *status)
 
ptwXYPointsptwXY_fromUnitbase (ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status)
 
ptwXPointsptwXY_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)
 
ptwXYPointsptwXY_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)
 
ptwXYPointsptwXY_valueTo_ptwXY (double x1, double x2, double y, nfu_status *status)
 
ptwXYPointsptwXY_createGaussianCenteredSigma1 (double accuracy, nfu_status *status)
 
ptwXYPointsptwXY_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)
 
ptwXYPointsptwXY_createFromFunction (int n, double *xs, ptwXY_createFromFunction_callback func, void *argList, double accuracy, int checkForRoots, int biSectionMax, nfu_status *status)
 
ptwXYPointsptwXY_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)
 
ptwXYPointsptwXY_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)
 
ptwXPointsptwXY_groupOneFunction (ptwXYPoints *ptwXY, ptwXPoints *groupBoundaries, ptwXY_group_normType normType, ptwXPoints *ptwX_norm, nfu_status *status)
 
ptwXPointsptwXY_groupTwoFunctions (ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, ptwXPoints *groupBoundaries, ptwXY_group_normType normType, ptwXPoints *ptwX_norm, nfu_status *status)
 
ptwXPointsptwXY_groupThreeFunctions (ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, ptwXYPoints *ptwXY3, ptwXPoints *groupBoundaries, ptwXY_group_normType normType, ptwXPoints *ptwX_norm, nfu_status *status)
 
ptwXPointsptwXY_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

#define ClosestAllowXFactor   10

Definition at line 25 of file ptwXY.h.

#define ptwXY_maxBiSectionMax   20

Definition at line 22 of file ptwXY.h.

#define ptwXY_minAccuracy   1e-14

Definition at line 23 of file ptwXY.h.

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

Definition at line 21 of file ptwXY.h.

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

Definition at line 20 of file ptwXY.h.

#define ptwXY_sectionSubdivideMax   1 << 16

Definition at line 24 of file ptwXY.h.

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

Definition at line 31 of file ptwXY.h.

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

Definition at line 33 of file ptwXY.h.

#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.

Typedef Documentation

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

Definition at line 66 of file ptwXY.h.

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

Definition at line 65 of file ptwXY.h.

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.

typedef enum ptwXY_sigma_e ptwXY_sigma
typedef struct ptwXYPoint_s ptwXYPoint
typedef struct ptwXYPoints_s ptwXYPoints

Enumeration Type Documentation

Enumerator
ptwXY_dataFrom_Unknown 
ptwXY_dataFrom_Points 
ptwXY_dataFrom_Overflow 

Definition at line 27 of file ptwXY.h.

Enumerator
ptwXY_group_normType_none 
ptwXY_group_normType_dx 
ptwXY_group_normType_norm 

Definition at line 28 of file ptwXY.h.

Enumerator
ptwXY_interpolationLinLin 
ptwXY_interpolationLinLog 
ptwXY_interpolationLogLin 
ptwXY_interpolationLogLog 
ptwXY_interpolationFlat 
ptwXY_interpolationOther 

Definition at line 35 of file ptwXY.h.

Enumerator
ptwXY_lessEqualGreaterX_empty 
ptwXY_lessEqualGreaterX_lessThan 
ptwXY_lessEqualGreaterX_equal 
ptwXY_lessEqualGreaterX_between 
ptwXY_lessEqualGreaterX_greater 

Definition at line 57 of file ptwXY.h.

Enumerator
ptwXY_sigma_none 
ptwXY_sigma_plusMinus 
ptwXY_sigma_Minus 
ptwXY_sigma_plus 

Definition at line 34 of file ptwXY.h.

Function Documentation

nfu_status ptwXY_abs ( ptwXYPoints ptwXY)

Definition at line 19 of file ptwXY_unitaryOperators.cc.

19  {
20 
21  int64_t i, nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
22  ptwXYPoint *p;
23  ptwXYOverflowPoint *o, *overflowHeader = &(ptwXY->overflowHeader);
24 
25  if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
26 
27  for( i = 0, p = ptwXY->points; i < nonOverflowLength; i++, p++ ) p->y = std::fabs( p->y );
28  for( o = overflowHeader->next; o != overflowHeader; o = o->next ) o->point.y = std::fabs( o->point.y );
29  return( ptwXY->status );
30 }
ptwXYPoint point
Definition: ptwXY.h:80
ptwXYPoint * points
Definition: ptwXY.h:99
const char * p
Definition: xmltok.h:285
int64_t ptwXY_getNonOverflowLength(ptwXYPoints const *ptwXY)
Definition: ptwXY_core.cc:590
double y
Definition: ptwXY.h:62
struct ptwXYOverflowPoint_s * next
Definition: ptwXY.h:78
nfu_status status
Definition: ptwXY.h:85
ptwXYOverflowPoint overflowHeader
Definition: ptwXY.h:98

Here is the call graph for this function:

nfu_status ptwXY_add_double ( ptwXYPoints ptwXY,
double  value 
)

Definition at line 40 of file ptwXY_binaryOperators.cc.

40 { return( ptwXY_slopeOffset( ptwXY, 1., value ) ); }
const XML_Char int const XML_Char * value
Definition: expat.h:331
nfu_status ptwXY_slopeOffset(ptwXYPoints *ptwXY, double slope, double offset)

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 138 of file ptwXY_binaryOperators.cc.

138  {
139 
140  ptwXYPoints *sum;
141 
142  if( ptwXY1->length == 0 ) {
143  sum = ptwXY_clone( ptwXY2, status ); }
144  else if( ptwXY2->length == 0 ) {
145  sum = ptwXY_clone( ptwXY1, status ); }
146  else {
147  sum = ptwXY_binary_ptwXY( ptwXY1, ptwXY2, 1., 1., 0., status );
148  }
149  return( sum );
150 }
ptwXYPoints * ptwXY_clone(ptwXYPoints *ptwXY, nfu_status *status)
Definition: ptwXY_core.cc:208
int64_t length
Definition: ptwXY.h:93
ptwXYPoints * ptwXY_binary_ptwXY(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, double v1, double v2, double v1v2, nfu_status *status)

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 1062 of file ptwXY_core.cc.

1062  {
1063 
1064  int64_t nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
1065  ptwXY_dataFrom dataFrom;
1066 
1067  if( ptwXY->length != 0 ) {
1068  double xMax = ptwXY_getXMaxAndFrom( ptwXY, &dataFrom );
1069  if( xMax >= x ) return( nfu_XNotAscending );
1070  }
1071 
1072  if( nonOverflowLength < ptwXY->allocatedSize ) { /* Room at end of points. Also handles the case when length = 0. */
1073  ptwXY->points[nonOverflowLength].x = x;
1074  ptwXY->points[nonOverflowLength].y = y; }
1075  else {
1076  if( ptwXY->overflowLength == ptwXY->overflowAllocatedSize ) {
1077  ptwXYPoint newPoint = { x, y };
1078  return( ptwXY_coalescePoints( ptwXY, ptwXY->length + ptwXY->overflowAllocatedSize, &newPoint, 0 ) ); }
1079  else { /* Add to end of overflow. */
1080  ptwXYOverflowPoint *overflowPoint = &(ptwXY->overflowPoints[ptwXY->overflowLength]);
1081 
1082  overflowPoint->prior = ptwXY->overflowHeader.prior;
1083  overflowPoint->next = overflowPoint->prior->next;
1084  overflowPoint->index = ptwXY->length;
1085  overflowPoint->prior->next = overflowPoint;
1086  overflowPoint->next->prior = overflowPoint;
1087  overflowPoint->point.x = x;
1088  overflowPoint->point.y = y;
1089  ptwXY->overflowLength++;
1090  }
1091  }
1092  ptwXY->length++;
1093  return( nfu_Okay );
1094 }
nfu_status ptwXY_coalescePoints(ptwXYPoints *ptwXY, int64_t size, ptwXYPoint *newPoint, int forceSmallerResize)
Definition: ptwXY_core.cc:469
ptwXYPoint point
Definition: ptwXY.h:80
ptwXYPoint * points
Definition: ptwXY.h:99
double ptwXY_getXMaxAndFrom(ptwXYPoints *ptwXY, ptwXY_dataFrom *dataFrom)
Definition: ptwXY_core.cc:1215
tuple x
Definition: test.py:50
struct ptwXYOverflowPoint_s * prior
Definition: ptwXY.h:77
int64_t length
Definition: ptwXY.h:93
int64_t ptwXY_getNonOverflowLength(ptwXYPoints const *ptwXY)
Definition: ptwXY_core.cc:590
enum ptwXY_dataFrom_e ptwXY_dataFrom
double y
Definition: ptwXY.h:62
struct ptwXYOverflowPoint_s * next
Definition: ptwXY.h:78
double x
Definition: ptwXY.h:62
int64_t overflowAllocatedSize
Definition: ptwXY.h:96
int64_t overflowLength
Definition: ptwXY.h:95
ptwXYOverflowPoint overflowHeader
Definition: ptwXY.h:98
int64_t index
Definition: ptwXY.h:79
ptwXYOverflowPoint * overflowPoints
Definition: ptwXY.h:100

Here is the call graph for this function:

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

Definition at line 146 of file ptwXY_misc.cc.

146  {
147 
148  int64_t i, originalLength = ptwXY1->length, notFirstPass = 0;
149  double y1, y2 = 0;
150  nfu_status status;
151  ptwXYPoint p1, p2;
152 
153  checkForRoots = checkForRoots && ptwXY1->biSectionMax;
154  if( ptwXY1->status != nfu_Okay ) return( ptwXY1->status );
157  if( ( status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( status );
158  for( i = originalLength - 1; i >= 0; i-- ) {
159  y1 = ptwXY1->points[i].y;
160  if( ( status = func( &(ptwXY1->points[i]), argList ) ) != nfu_Okay ) return( status );
161  p1 = ptwXY1->points[i];
162  if( notFirstPass ) {
163  if( ( status = ptwXY_applyFunction2( ptwXY1, y1, y2, &p1, &p2, func, argList, 0, checkForRoots ) ) != nfu_Okay ) return( status );
164  }
165  notFirstPass = 1;
166  p2 = p1;
167  y2 = y1;
168  }
169  ptwXY_update_biSectionMax( ptwXY1, (double) originalLength );
170  return( status );
171 }
static nfu_status ptwXY_applyFunction2(ptwXYPoints *ptwXY1, double y1, double y2, ptwXYPoint *p1, ptwXYPoint *p2, ptwXY_applyFunction_callback func, void *argList, int level, int checkForRoots)
Definition: ptwXY_misc.cc:175
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
ptwXYPoint * points
Definition: ptwXY.h:99
double biSectionMax
Definition: ptwXY.h:90
int64_t length
Definition: ptwXY.h:93
enum nfu_status_e nfu_status
double y
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
nfu_status status
Definition: ptwXY.h:85
void ptwXY_update_biSectionMax(ptwXYPoints *ptwXY1, double oldLength)
Definition: ptwXY_misc.cc:31

Here is the call graph for this function:

Here is the caller graph for this function:

nfu_status ptwXY_areDomainsMutual ( ptwXYPoints ptwXY1,
ptwXYPoints ptwXY2 
)

Definition at line 257 of file ptwXY_convenient.cc.

257  {
258 
259  nfu_status status;
260  int64_t n1 = ptwXY1->length, n2 = ptwXY2->length;
261  ptwXYPoint *xy1, *xy2;
262 
263  if( ( status = ptwXY1->status ) != nfu_Okay ) return( status );
264  if( ( status = ptwXY2->status ) != nfu_Okay ) return( status );
265  if( n1 == 0 ) return( nfu_empty );
266  if( n2 == 0 ) return( nfu_empty );
267  if( n1 < 2 ) {
268  status = nfu_tooFewPoints; }
269  else if( n2 < 2 ) {
270  status = nfu_tooFewPoints; }
271  else {
272  xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, 0 );
273  xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, 0 );
274  if( xy1->x < xy2->x ) {
275  if( xy2->y != 0. ) status = nfu_domainsNotMutual; }
276  else if( xy1->x > xy2->x ) {
277  if( xy1->y != 0. ) status = nfu_domainsNotMutual;
278  }
279 
280  if( status == nfu_Okay ) {
281  xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, n1 - 1 );
282  xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, n2 - 1 );
283  if( xy1->x < xy2->x ) {
284  if( xy1->y != 0. ) status = nfu_domainsNotMutual; }
285  else if( xy1->x > xy2->x ) {
286  if( xy2->y != 0. ) status = nfu_domainsNotMutual;
287  }
288  }
289  }
290  return( status );
291 }
ptwXYPoint * ptwXY_getPointAtIndex_Unsafely(ptwXYPoints *ptwXY, int64_t index)
Definition: ptwXY_core.cc:684
int64_t length
Definition: ptwXY.h:93
enum nfu_status_e nfu_status
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
nfu_status status
Definition: ptwXY.h:85

Here is the call graph for this function:

Here is the caller graph for this function:

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.

108  {
109 
110  int64_t i;
111  int unionOptions = ptwXY_union_fill | ptwXY_union_mergeClosePoints;
112  double y;
113  ptwXYPoints *n;
114  ptwXYPoint *p;
115 
116  *status = nfu_otherInterpolation;
117  if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
118  if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( NULL );
119  if( ( *status = ptwXY_areDomainsMutual( ptwXY1, ptwXY2 ) ) != nfu_Okay ) return( NULL );
120  if( ( ptwXY1->interpolation == ptwXY_interpolationFlat ) || ( ptwXY2->interpolation == ptwXY_interpolationFlat ) ) {
121  *status = nfu_invalidInterpolation;
122  if( ( ptwXY1->interpolation != ptwXY2->interpolation ) ) return( NULL );
123  }
124  if( ( n = ptwXY_union( ptwXY1, ptwXY2, status, unionOptions ) ) != NULL ) {
125  for( i = 0, p = n->points; i < n->length; i++, p++ ) {
126  if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, p->x, &y ) ) != nfu_Okay ) goto Err;
127  p->y = v1 * p->y + v2 * y + v1v2 * y * p->y;
128  }
129  }
130  return( n );
131 Err:
132  if( n ) ptwXY_free( n );
133  return( NULL );
134 }
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
ptwXYPoint * points
Definition: ptwXY.h:99
const char * p
Definition: xmltok.h:285
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
#define ptwXY_union_fill
Definition: ptwXY.h:31
ptwXYPoints * ptwXY_union(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int unionOptions)
static nfu_status ptwXY_getValueAtX_ignore_XOutsideDomainError(ptwXYPoints *ptwXY1, double x, double *y)
const G4int n
nfu_status ptwXY_areDomainsMutual(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2)
#define ptwXY_union_mergeClosePoints
Definition: ptwXY.h:33
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62

Here is the call graph for this function:

Here is the caller graph for this function:

nfu_status ptwXY_clear ( ptwXYPoints ptwXY)

Definition at line 536 of file ptwXY_core.cc.

536  {
537 
538  if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
539 
540  ptwXY->length = 0;
541  ptwXY->overflowLength = 0;
542  ptwXY->overflowHeader.prior = &(ptwXY->overflowHeader);
543  ptwXY->overflowHeader.next = &(ptwXY->overflowHeader);
544  return( nfu_Okay );
545 }
struct ptwXYOverflowPoint_s * prior
Definition: ptwXY.h:77
int64_t length
Definition: ptwXY.h:93
struct ptwXYOverflowPoint_s * next
Definition: ptwXY.h:78
int64_t overflowLength
Definition: ptwXY.h:95
nfu_status status
Definition: ptwXY.h:85
ptwXYOverflowPoint overflowHeader
Definition: ptwXY.h:98

Here is the caller graph for this function:

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

Definition at line 25 of file ptwXY_methods.cc.

25  {
26 /*
27  This function acts oddly for xy = [ [ 1, 0 ], [ 3, -2 ], [ 4, 1 ] ] and yMin = 0.2, why???????
28  This function probably only works for linear, linear interpolation (mainly because of ptwXY_clip2).
29 */
30  int64_t i, j, n;
31  double x2, y2;
32  nfu_status status;
33  ptwXYPoints *clipped;
34  ptwXYPoint *points;
35 
36  if( ( status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( status );
38  n = ptwXY1->length;
39  if( n > 0 ) {
40  i = 0;
41  if( ptwXY_getYMax( ptwXY1 ) < yMin ) i = 1;
42  if( ptwXY_getYMin( ptwXY1 ) > yMax ) i = 1;
43  if( i == 1 ) return( ptwXY_clear( ptwXY1 ) );
44  }
45  if( n == 1 ) {
46  y2 = ptwXY1->points[0].y;
47  if( y2 < yMin ) {
48  ptwXY1->points[0].y = yMin; }
49  else if( y2 > yMax ) {
50  ptwXY1->points[0].y = yMax;
51  } }
52  else if( n > 1 ) {
53  if( ( clipped = ptwXY_new( ptwXY1->interpolation, &(ptwXY1->interpolationOtherInfo),
54  ptwXY1->biSectionMax, ptwXY1->accuracy, n, 10, &status, ptwXY1->userFlag ) ) == NULL )
55  return( ptwXY1->status = status );
56  for( i = 0; i < n; i++ ) {
57  x2 = ptwXY1->points[i].x;
58  y2 = ptwXY1->points[i].y;
59  if( y2 < yMin ) {
60  if( i > 0 ) {
61  points = ptwXY_getPointAtIndex_Unsafely( clipped, clipped->length - 1 );
62  if( points->y > yMin ) {
63  if( ( status = ptwXY_clip2( clipped, yMin, points->x, points->y, x2, y2 ) ) != nfu_Okay ) goto Err;
64  }
65  }
66  if( ( status = ptwXY_setValueAtX( clipped, x2, yMin ) ) != nfu_Okay ) goto Err;
67  j = i;
68  for( i++; i < n; i++ ) if( !( ptwXY1->points[i].y < yMin ) ) break;
69  if( i < n ) {
70  x2 = ptwXY1->points[i].x;
71  y2 = ptwXY1->points[i].y;
72  if( ( status = ptwXY_clip2( clipped, yMin, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) goto Err;
73  if( y2 > yMax ) {
74  if( ( status = ptwXY_clip2( clipped, yMax, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) goto Err;
75  } }
76  else if( j != n - 1 ) {
77  if( ( status = ptwXY_setValueAtX( clipped, ptwXY1->points[n - 1].x, yMin ) ) != nfu_Okay ) goto Err;
78  }
79  i--; }
80  else if( y2 > yMax ) {
81  if( i > 0 ) {
82  points = ptwXY_getPointAtIndex_Unsafely( clipped, clipped->length - 1 );
83  if( points->y < yMax ) {
84  if( ( status = ptwXY_clip2( clipped, yMax, points->x, points->y, x2, y2 ) ) != nfu_Okay ) goto Err;
85  }
86  }
87  if( ( status = ptwXY_setValueAtX( clipped, x2, yMax ) ) != nfu_Okay ) goto Err;
88  j = i;
89  for( i++; i < n; i++ ) if( !( ptwXY1->points[i].y > yMax ) ) break;
90  if( i < n ) {
91  x2 = ptwXY1->points[i].x;
92  y2 = ptwXY1->points[i].y;
93  if( ( status = ptwXY_clip2( clipped, yMax, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) goto Err;
94  if( y2 < yMin ) {
95  if( ( status = ptwXY_clip2( clipped, yMin, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) goto Err;
96  } }
97  else if( j != n - 1 ) {
98  if( ( status = ptwXY_setValueAtX( clipped, ptwXY1->points[n - 1].x, yMax ) ) != nfu_Okay ) goto Err;
99  }
100  i--; }
101  else {
102  if( ( status = ptwXY_setValueAtX( clipped, x2, y2 ) ) != nfu_Okay ) goto Err;
103  }
104  }
105  if( ( status = ptwXY_simpleCoalescePoints( clipped ) ) != nfu_Okay ) goto Err;
106  ptwXY1->length = clipped->length; /* The squeamish may want to skip the next few lines. */
107  clipped->length = n;
108  n = ptwXY1->allocatedSize;
109  ptwXY1->allocatedSize = clipped->allocatedSize;
110  clipped->allocatedSize = n;
111  points = clipped->points;
112  clipped->points = ptwXY1->points;
113  ptwXY1->points = points;
114  ptwXY_free( clipped );
115  }
116 
117  return( ptwXY1->status );
118 
119 Err:
120  ptwXY_free( clipped );
121  return( ptwXY1->status = status );
122 }
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
int userFlag
Definition: ptwXY.h:89
ptwXYPoint * points
Definition: ptwXY.h:99
ptwXYPoint * ptwXY_getPointAtIndex_Unsafely(ptwXYPoints *ptwXY, int64_t index)
Definition: ptwXY_core.cc:684
double ptwXY_getYMin(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:1248
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
Definition: ptwXY_core.cc:876
double biSectionMax
Definition: ptwXY.h:90
int64_t length
Definition: ptwXY.h:93
double ptwXY_getYMax(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:1269
static nfu_status ptwXY_clip2(ptwXYPoints *ptwXY1, double y, double x1, double y1, double x2, double y2)
enum nfu_status_e nfu_status
const G4int n
ptwXY_interpolationOtherInfo interpolationOtherInfo
Definition: ptwXY.h:88
int64_t allocatedSize
Definition: ptwXY.h:94
nfu_status ptwXY_clear(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:536
double y
Definition: ptwXY.h:62
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: ptwXY_core.cc:29
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
nfu_status status
Definition: ptwXY.h:85
double accuracy
Definition: ptwXY.h:91

Here is the call graph for this function:

ptwXYPoints* ptwXY_clone ( ptwXYPoints ptwXY,
nfu_status status 
)

Definition at line 208 of file ptwXY_core.cc.

208  {
209 
210  return( ptwXY_slice( ptwXY, 0, ptwXY->length, ptwXY->overflowAllocatedSize, status ) );
211 }
int64_t length
Definition: ptwXY.h:93
int64_t overflowAllocatedSize
Definition: ptwXY.h:96
ptwXYPoints * ptwXY_slice(ptwXYPoints *ptwXY, int64_t index1, int64_t index2, int64_t secondarySize, nfu_status *status)
Definition: ptwXY_core.cc:248

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 215 of file ptwXY_core.cc.

215  {
216 
217  ptwXYPoints *n1;
218 
219  if( interpolationTo == ptwXY_interpolationOther ) {
220  *status = nfu_otherInterpolation;
221  return( NULL );
222  }
223  if( ( n1 = ptwXY_clone( ptwXY, status ) ) != NULL ) {
225  n1->interpolation = interpolationTo;
226  switch( interpolationTo ) {
237  case ptwXY_interpolationOther : /* Does not happen, but needed to stop compilers from complaining. */
238  break;
239  }
241  n1->interpolationOtherInfo.argList = NULL;
242  }
243  return( n1 );
244 }
static char const logLinInterpolationString[]
Definition: ptwXY_core.cc:19
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
static char const logLogInterpolationString[]
Definition: ptwXY_core.cc:20
ptwXYPoints * ptwXY_clone(ptwXYPoints *ptwXY, nfu_status *status)
Definition: ptwXY_core.cc:208
ptwXY_getValue_callback getValueFunc
Definition: ptwXY.h:71
static char const linLinInterpolationString[]
Definition: ptwXY_core.cc:17
void * nfu_free(void *p)
ptwXY_interpolationOtherInfo interpolationOtherInfo
Definition: ptwXY.h:88
char const * interpolationString
Definition: ptwXY.h:70
static char const flatInterpolationString[]
Definition: ptwXY_core.cc:21
static char const linLogInterpolationString[]
Definition: ptwXY_core.cc:18

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 469 of file ptwXY_core.cc.

469  {
470 
471  int addNewPoint;
472  int64_t length = ptwXY->length + ( ( newPoint != NULL ) ? 1 : 0 );
473  ptwXYOverflowPoint *last = ptwXY->overflowHeader.prior;
474  ptwXYPoint *pointsFrom, *pointsTo;
475 
476  if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
477  if( ptwXY->overflowLength == 0 ) return( nfu_Okay );
478 
479  if( size < length ) size = length;
480  if( size > ptwXY->allocatedSize ) {
481  if( ptwXY_reallocatePoints( ptwXY, size, forceSmallerResize ) != nfu_Okay ) return( ptwXY->status );
482  }
483  pointsFrom = &(ptwXY->points[ptwXY_getNonOverflowLength( ptwXY ) - 1]);
484  pointsTo = &(ptwXY->points[length - 1]);
485  while( last != &(ptwXY->overflowHeader) ) {
486  addNewPoint = 0;
487  if( newPoint != NULL ) {
488  if( ( pointsFrom >= ptwXY->points ) && ( pointsFrom->x > last->point.x ) ) {
489  if( newPoint->x > pointsFrom->x ) addNewPoint = 1; }
490  else {
491  if( newPoint->x > last->point.x ) addNewPoint = 1;
492  }
493  if( addNewPoint == 1 ) {
494  *pointsTo = *newPoint;
495  newPoint = NULL;
496  }
497  }
498  if( addNewPoint == 0 ) {
499  if( ( pointsFrom >= ptwXY->points ) && ( pointsFrom->x > last->point.x ) ) {
500  *pointsTo = *pointsFrom;
501  pointsFrom--; }
502  else {
503  *pointsTo = last->point;
504  last = last->prior;
505  }
506  }
507  pointsTo--;
508  } // Loop checking, 11.06.2015, T. Koi
509  while( ( newPoint != NULL ) && ( pointsFrom >= ptwXY->points ) ) {
510  if( newPoint->x > pointsFrom->x ) {
511  *pointsTo = *newPoint;
512  newPoint = NULL; }
513  else {
514  *pointsTo = *pointsFrom;
515  pointsFrom--;
516  }
517  pointsTo--;
518  } // Loop checking, 11.06.2015, T. Koi
519  if( newPoint != NULL ) *pointsTo = *newPoint;
520  ptwXY->overflowHeader.prior = &(ptwXY->overflowHeader);
521  ptwXY->overflowHeader.next = &(ptwXY->overflowHeader);
522  ptwXY->length = length;
523  ptwXY->overflowLength = 0;
524  return( nfu_Okay );
525 }
ptwXYPoint point
Definition: ptwXY.h:80
ptwXYPoint * points
Definition: ptwXY.h:99
struct ptwXYOverflowPoint_s * prior
Definition: ptwXY.h:77
int64_t length
Definition: ptwXY.h:93
int64_t ptwXY_getNonOverflowLength(ptwXYPoints const *ptwXY)
Definition: ptwXY_core.cc:590
int64_t allocatedSize
Definition: ptwXY.h:94
struct ptwXYOverflowPoint_s * next
Definition: ptwXY.h:78
double x
Definition: ptwXY.h:62
int64_t overflowLength
Definition: ptwXY.h:95
nfu_status status
Definition: ptwXY.h:85
ptwXYOverflowPoint overflowHeader
Definition: ptwXY.h:98
nfu_status ptwXY_reallocatePoints(ptwXYPoints *ptwXY, int64_t size, int forceSmallerResize)
Definition: ptwXY_core.cc:410

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 96 of file ptwXY_functions.cc.

96  {
97 /*
98 * Currently, only supports linear-linear interpolation.
99 *
100 * This function calculates c(y) = integral dx f1(x) * f2(y-x)
101 *
102 */
103  int64_t i1, i2, n1, n2, n;
104  ptwXYPoints *f1 = ptwXY1, *f2 = ptwXY2, *convolute;
105  double accuracy = ptwXY1->accuracy, yMin, yMax, c, y, dy;
106 
107  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
108  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY2 ) ) != nfu_Okay ) return( NULL );
109 
111  if( ( ptwXY1->interpolation != ptwXY_interpolationLinLin ) || ( ptwXY2->interpolation != ptwXY_interpolationLinLin ) ) return( NULL );
112  *status = nfu_Okay;
113 
114  n1 = f1->length;
115  n2 = f2->length;
116 
117  if( ( n1 == 0 ) || ( n2 == 0 ) ) {
118  convolute = ptwXY_new( ptwXY_interpolationLinLin, NULL, 1., accuracy, 0, 0, status, 0 );
119  return( convolute );
120  }
121 
122  if( ( n1 == 1 ) || ( n2 == 1 ) ) {
123  *status = nfu_tooFewPoints;
124  return( NULL );
125  }
126 
127  if( accuracy < ptwXY2->accuracy ) accuracy = ptwXY2->accuracy;
128  n = n1 * n2;
129  if( mode == 0 ) {
130  mode = 1;
131  if( n > 1000 ) mode = -1;
132  }
133  if( n > 100000 ) mode = -1;
134  if( ( convolute = ptwXY_new( ptwXY_interpolationLinLin, NULL, 1., accuracy, 400, 40, status, 0 ) ) == NULL ) return( NULL );
135 
136  yMin = f1->points[0].x + f2->points[0].x;
137  yMax = f1->points[n1 - 1].x + f2->points[n2 - 1].x;
138 
139  if( ( *status = ptwXY_setValueAtX( convolute, yMin, 0. ) ) != nfu_Okay ) goto Err;
140 
141  if( mode < 0 ) {
142  dy = ( yMax - yMin ) / 400;
143  for( y = yMin + dy; y < yMax; y += dy ) {
144  if( ( *status = ptwXY_convolution2( f1, f2, y, yMin, &c ) ) != nfu_Okay ) goto Err;
145  if( ( *status = ptwXY_setValueAtX( convolute, y, c ) ) != nfu_Okay ) goto Err;
146  } }
147  else {
148  for( i1 = 0; i1 < n1; i1++ ) {
149  for( i2 = 0; i2 < n2; i2++ ) {
150  y = yMin + ( f1->points[i1].x - f1->points[0].x ) + ( f2->points[i2].x - f2->points[0].x );
151  if( y <= yMin ) continue;
152  if( y >= yMax ) continue;
153  if( ( *status = ptwXY_convolution2( f1, f2, y, yMin, &c ) ) != nfu_Okay ) goto Err;
154  if( ( *status = ptwXY_setValueAtX( convolute, y, c ) ) != nfu_Okay ) goto Err;
155  }
156  }
157  }
158  if( ( *status = ptwXY_setValueAtX( convolute, yMax, 0. ) ) != nfu_Okay ) goto Err;
159  if( ( *status = ptwXY_simpleCoalescePoints( convolute ) ) != nfu_Okay ) goto Err;
160  for( i1 = convolute->length - 1; i1 > 0; i1-- ) {
161  if( ( *status = ptwXY_convolution3( convolute, f1, f2, convolute->points[i1 - 1].x, convolute->points[i1 - 1].y,
162  convolute->points[i1].x, convolute->points[i1].y, yMin ) ) != nfu_Okay ) goto Err;
163  }
164 
165  return( convolute );
166 
167 Err:
168  ptwXY_free( convolute );
169  return( NULL );
170 }
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
static nfu_status ptwXY_convolution3(ptwXYPoints *convolute, ptwXYPoints *f1, ptwXYPoints *f2, double y1, double c1, double y2, double c2, double yMin)
ptwXYPoint * points
Definition: ptwXY.h:99
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
Definition: ptwXY_core.cc:876
int64_t length
Definition: ptwXY.h:93
const G4int n
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: ptwXY_core.cc:29
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
tuple c
Definition: test.py:13
double accuracy
Definition: ptwXY.h:91
static nfu_status ptwXY_convolution2(ptwXYPoints *f1, ptwXYPoints *f2, double y, double yMin, double *c)

Here is the call graph for this function:

nfu_status ptwXY_copy ( ptwXYPoints dest,
ptwXYPoints src 
)

Definition at line 148 of file ptwXY_core.cc.

148  {
149 
150  int64_t i, nonOverflowLength = ptwXY_getNonOverflowLength( src );
151  ptwXYPoint *pointFrom, *pointTo;
152  ptwXYOverflowPoint *o, *overflowHeader = &(src->overflowHeader);
153 
154  if( dest->status != nfu_Okay ) return( dest->status );
155  if( src->status != nfu_Okay ) return( src->status );
156 
157  ptwXY_clear( dest );
158  if( dest->interpolation == ptwXY_interpolationOther ) {
159  if( dest->interpolationOtherInfo.interpolationString != NULL ) {
161  }
162  }
163  dest->interpolation = ptwXY_interpolationLinLin; /* This and prior lines are in case interpolation is 'other' and ptwXY_reallocatePoints fails. */
164  if( dest->allocatedSize < src->length ) ptwXY_reallocatePoints( dest, src->length, 0 );
165  if( dest->status != nfu_Okay ) return( dest->status );
166  dest->interpolation = src->interpolation;
167  if( dest->interpolation == ptwXY_interpolationOther ) {
168  if( src->interpolationOtherInfo.interpolationString != NULL ) {
170  return( dest->status = nfu_mallocError );
171  } }
172  else {
174  }
177  dest->userFlag = src->userFlag;
178  dest->biSectionMax = src->biSectionMax;
179  dest->accuracy = src->accuracy;
180  dest->minFractional_dx = src->minFractional_dx;
181  pointFrom = src->points;
182  o = src->overflowHeader.next;
183  pointTo = dest->points;
184  i = 0;
185  while( o != overflowHeader ) {
186  if( i < nonOverflowLength ) {
187  if( pointFrom->x < o->point.x ) {
188  *pointTo = *pointFrom;
189  i++;
190  pointFrom++; }
191  else {
192  *pointTo = o->point;
193  o = o->next;
194  } }
195  else {
196  *pointTo = o->point;
197  o = o->next;
198  }
199  pointTo++;
200  } // Loop checking, 11.06.2015, T. Koi
201  for( ; i < nonOverflowLength; i++, pointFrom++, pointTo++ ) *pointTo = *pointFrom;
202  dest->length = src->length;
203  return( dest->status );
204 }
double minFractional_dx
Definition: ptwXY.h:92
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
ptwXYPoint point
Definition: ptwXY.h:80
int userFlag
Definition: ptwXY.h:89
ptwXYPoint * points
Definition: ptwXY.h:99
ptwXY_getValue_callback getValueFunc
Definition: ptwXY.h:71
double biSectionMax
Definition: ptwXY.h:90
int64_t length
Definition: ptwXY.h:93
int64_t ptwXY_getNonOverflowLength(ptwXYPoints const *ptwXY)
Definition: ptwXY_core.cc:590
void * nfu_free(void *p)
ptwXY_interpolationOtherInfo interpolationOtherInfo
Definition: ptwXY.h:88
int64_t allocatedSize
Definition: ptwXY.h:94
nfu_status ptwXY_clear(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:536
struct ptwXYOverflowPoint_s * next
Definition: ptwXY.h:78
double x
Definition: ptwXY.h:62
nfu_status status
Definition: ptwXY.h:85
char const * interpolationString
Definition: ptwXY.h:70
ptwXYOverflowPoint overflowHeader
Definition: ptwXY.h:98
double accuracy
Definition: ptwXY.h:91
nfu_status ptwXY_reallocatePoints(ptwXYPoints *ptwXY, int64_t size, int forceSmallerResize)
Definition: ptwXY_core.cc:410

Here is the call graph for this function:

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.

424  {
425 
426  int64_t i;
427  double *d = xys;
428  nfu_status status;
429  ptwXYPoint *pointFrom;
430 
431  if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
432  if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status );
433 
434  if( index1 < 0 ) index1 = 0;
435  if( index2 > ptwXY->length ) index2 = ptwXY->length;
436  if( index2 < index1 ) index2 = index1;
437  *numberOfPoints = index2 - index1;
438  if( allocatedSize < ( index2 - index1 ) ) return( nfu_insufficientMemory );
439 
440  for( i = index1, pointFrom = ptwXY->points; i < index2; i++, pointFrom++ ) {
441  *(d++) = pointFrom->x;
442  *(d++) = pointFrom->y;
443  }
444 
445  return( nfu_Okay );
446 }
ptwXYPoint * points
Definition: ptwXY.h:99
int64_t length
Definition: ptwXY.h:93
enum nfu_status_e nfu_status
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
nfu_status status
Definition: ptwXY.h:85

Here is the call graph for this function:

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.

110  {
111 
112  ptwXYPoints *ptwXY;
113 
114  if( primarySize < length ) primarySize = length;
115  if( ( ptwXY = ptwXY_new( interpolation, interpolationOtherInfo, biSectionMax, accuracy, primarySize,
116  secondarySize, status, userFlag ) ) != NULL ) {
117  if( ( *status = ptwXY_setXYData( ptwXY, length, xy ) ) != nfu_Okay ) {
118  ptwXY = ptwXY_free( ptwXY );
119  }
120  }
121  return( ptwXY );
122 }
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
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: ptwXY_core.cc:29
nfu_status ptwXY_setXYData(ptwXYPoints *ptwXY, int64_t length, double const *xy)
Definition: ptwXY_core.cc:597

Here is the call graph for this function:

Here is the caller graph for this function:

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.

128  {
129 
130  int i;
131  ptwXYPoints *ptwXY;
132 
133  if( primarySize < length ) primarySize = length;
134  if( ( ptwXY = ptwXY_new( interpolation, interpolationOtherInfo, biSectionMax, accuracy, primarySize,
135  secondarySize, status, userFlag ) ) != NULL ) {
136  for( i = 0; i < length; i++ ) {
137  ptwXY->points[i].x = Xs[i];
138  ptwXY->points[i].y = Ys[i];
139  }
140  ptwXY->length = length;
141  }
142 
143  return( ptwXY );
144 }
ptwXYPoint * points
Definition: ptwXY.h:99
int64_t length
Definition: ptwXY.h:93
double y
Definition: ptwXY.h:62
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: ptwXY_core.cc:29
double x
Definition: ptwXY.h:62

Here is the call graph for this function:

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.

41  {
42 
43  int64_t i;
44  double x1, y1, x2, y2, eps = ClosestAllowXFactor * DBL_EPSILON;
45  ptwXYPoints *ptwXY;
46  ptwXYPoint *p1, *p2;
47 
48  *status = nfu_Okay;
49  if( n < 2 ) { *status = nfu_tooFewPoints; return( NULL ); }
50  for( i = 1; i < n; i++ ) {
51  if( xs[i-1] >= xs[i] ) *status = nfu_XNotAscending;
52  }
53  if( *status == nfu_XNotAscending ) return( NULL );
54 
55  x1 = xs[0];
56  if( ( *status = func( x1, &y1, argList ) ) != nfu_Okay ) return( NULL );
57  if( ( ptwXY = ptwXY_new( ptwXY_interpolationLinLin, NULL, biSectionMax, accuracy, 500, 50, status, 0 ) ) == NULL ) return( NULL );
58  for( i = 1; i < n; i++ ) {
59  if( ( *status = ptwXY_setValueAtX_overrideIfClose( ptwXY, x1, y1, eps, 0 ) ) != nfu_Okay ) goto err;
60  x2 = xs[i];
61  if( ( *status = func( x2, &y2, argList ) ) != nfu_Okay ) goto err;
62  if( ( *status = ptwXY_createFromFunctionBisect( ptwXY, x1, y1, x2, y2, func, argList, 0, checkForRoots, eps ) ) != nfu_Okay ) goto err;
63  x1 = x2;
64  y1 = y2;
65  }
66  if( ( *status = ptwXY_setValueAtX_overrideIfClose( ptwXY, x2, y2, eps, 1 ) ) != nfu_Okay ) goto err;
67 
68  if( checkForRoots ) {
69  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) goto err;
70  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. */
71  p1 = &(ptwXY->points[i]);
72  if( p2 != NULL ) {
73  if( ( p1->y * p2->y ) < 0. ) {
74  if( ( *status = ptwXY_createFromFunctionZeroCrossing( ptwXY, p1->x, p1->y, p2->x, p2->y, func, argList, eps ) ) != nfu_Okay ) goto err;
75  }
76  }
77  }
78  }
79 
80  return( ptwXY );
81 
82 err:
83  ptwXY_free( ptwXY );
84  return( NULL );
85 }
#define ClosestAllowXFactor
Definition: ptwXY.h:25
static nfu_status ptwXY_createFromFunctionZeroCrossing(ptwXYPoints *ptwXY, double x1, double y1, double x2, double y2, ptwXY_createFromFunction_callback func, void *argList, double eps)
Definition: ptwXY_misc.cc:116
static nfu_status ptwXY_createFromFunctionBisect(ptwXYPoints *ptwXY, double x1, double y1, double x2, double y2, ptwXY_createFromFunction_callback func, void *argList, int level, int checkForRoots, double eps)
Definition: ptwXY_misc.cc:97
ptwXYPoint * points
Definition: ptwXY.h:99
static const G4double eps
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
int64_t length
Definition: ptwXY.h:93
nfu_status ptwXY_setValueAtX_overrideIfClose(ptwXYPoints *ptwXY, double x, double y, double eps, int override)
Definition: ptwXY_core.cc:883
#define DBL_EPSILON
Definition: templates.hh:87
const G4int n
double y
Definition: ptwXY.h:62
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: ptwXY_core.cc:29
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529

Here is the call graph for this function:

Here is the caller graph for this function:

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.

90  {
91 
92  return( ptwXY_createFromFunction( (int) xs->length, xs->points, func, argList, accuracy, checkForRoots, biSectionMax, status ) );
93 }
ptwXYPoints * ptwXY_createFromFunction(int n, double *xs, ptwXY_createFromFunction_callback func, void *argList, double accuracy, int checkForRoots, int biSectionMax, nfu_status *status)
Definition: ptwXY_misc.cc:40
int64_t length
Definition: ptwX.h:26
double * points
Definition: ptwX.h:29

Here is the call graph for this function:

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.

567  {
568 
569  int64_t i;
570  ptwXYPoints *gaussian, *sliced;
571  ptwXYPoint *point;
572 
573  if( ( gaussian = ptwXY_createGaussianCenteredSigma1( accuracy, status ) ) == NULL ) return( NULL );
574  for( i = 0, point = gaussian->points; i < gaussian->length; i++, point++ ) {
575  point->x = point->x * sigma + xCenter;
576  point->y *= amplitude;
577  }
578  if( ( gaussian->points[0].x < xMin ) || ( gaussian->points[gaussian->length - 1].x > xMax ) ) {
579  if( ( sliced = ptwXY_xSlice( gaussian, xMin, xMax, 10, 1, status ) ) == NULL ) goto Err;
580  ptwXY_free( gaussian );
581  gaussian = sliced;
582  }
583 
584  return( gaussian );
585 
586 Err:
587  ptwXY_free( gaussian );
588  return( NULL );
589 }
ptwXYPoint * points
Definition: ptwXY.h:99
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
int64_t length
Definition: ptwXY.h:93
ptwXYPoints * ptwXY_createGaussianCenteredSigma1(double accuracy, nfu_status *status)
double y
Definition: ptwXY.h:62
ptwXYPoints * ptwXY_xSlice(ptwXYPoints *ptwXY, double xMin, double xMax, int64_t secondarySize, int fill, nfu_status *status)
Definition: ptwXY_core.cc:274
double x
Definition: ptwXY.h:62

Here is the call graph for this function:

ptwXYPoints* ptwXY_createGaussianCenteredSigma1 ( double  accuracy,
nfu_status status 
)

Definition at line 492 of file ptwXY_convenient.cc.

492  {
493 
494  int64_t i, n;
495  ptwXYPoint *pm, *pp;
496  double x1, y1, x2, y2, accuracy2, yMin = 1e-10;
497  ptwXYPoints *gaussian;
498 
499  if( accuracy < 1e-5 ) accuracy = 1e-5;
500  if( accuracy > 1e-1 ) accuracy = 1e-1;
501  if( ( gaussian = ptwXY_new( ptwXY_interpolationLinLin, NULL, 1., accuracy, 200, 100, status, 0 ) ) == NULL ) return( NULL );
502  accuracy2 = accuracy = gaussian->accuracy;
503  if( accuracy2 > 5e-3 ) accuracy2 = 5e-3;
504 
505  x1 = -std::sqrt( -2. * G4Log( yMin ) );
506  y1 = yMin;
507  x2 = -5.2;
508  y2 = G4Exp( -0.5 * x2 * x2 );
509  if( ( *status = ptwXY_setValueAtX( gaussian, x1, y1 ) ) != nfu_Okay ) goto Err;
510  gaussian->accuracy = 20 * accuracy2;
511  if( ( *status = ptwXY_createGaussianCenteredSigma1_2( gaussian, x1, y1, x2, y2, 1 ) ) != nfu_Okay ) goto Err;
512  x1 = x2;
513  y1 = y2;
514  x2 = -4.;
515  y2 = G4Exp( -0.5 * x2 * x2 );
516  gaussian->accuracy = 5 * accuracy2;
517  if( ( *status = ptwXY_createGaussianCenteredSigma1_2( gaussian, x1, y1, x2, y2, 1 ) ) != nfu_Okay ) goto Err;
518  x1 = x2;
519  y1 = y2;
520  x2 = -1;
521  y2 = G4Exp( -0.5 * x2 * x2 );
522  gaussian->accuracy = accuracy;
523  if( ( *status = ptwXY_createGaussianCenteredSigma1_2( gaussian, x1, y1, x2, y2, 1 ) ) != nfu_Okay ) goto Err;
524  x1 = x2;
525  y1 = y2;
526  x2 = 0;
527  y2 = G4Exp( -0.5 * x2 * x2 );
528  if( ( *status = ptwXY_createGaussianCenteredSigma1_2( gaussian, x1, y1, x2, y2, 1 ) ) != nfu_Okay ) goto Err;
529 
530  n = gaussian->length;
531  if( ( *status = ptwXY_coalescePoints( gaussian, 2 * n + 1, NULL, 0 ) ) != nfu_Okay ) goto Err;
532  if( ( *status = ptwXY_setValueAtX( gaussian, 0., 1. ) ) != nfu_Okay ) goto Err;
533  pp = &(gaussian->points[gaussian->length]);
534  for( i = 0, pm = pp - 2; i < n; i++, pp++, pm-- ) {
535  *pp = *pm;
536  pp->x *= -1;
537  }
538  gaussian->length = 2 * n + 1;
539 
540  return( gaussian );
541 
542 Err:
543  ptwXY_free( gaussian );
544  return( NULL );
545 }
nfu_status ptwXY_coalescePoints(ptwXYPoints *ptwXY, int64_t size, ptwXYPoint *newPoint, int forceSmallerResize)
Definition: ptwXY_core.cc:469
ptwXYPoint * points
Definition: ptwXY.h:99
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
Definition: ptwXY_core.cc:876
int64_t length
Definition: ptwXY.h:93
const G4int n
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
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: ptwXY_core.cc:29
double x
Definition: ptwXY.h:62
static nfu_status ptwXY_createGaussianCenteredSigma1_2(ptwXYPoints *ptwXY, double x1, double y1, double x2, double y2, int addX1Point)
double accuracy
Definition: ptwXY.h:91

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 660 of file ptwXY_core.cc.

660  {
661 
662  int64_t n = ptwXY->length - ( i2 - i1 );
663 
664  if( ( ptwXY->status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( ptwXY->status );
665  if( ( i1 < 0 ) || ( i1 > i2 ) || ( i2 > ptwXY->length ) ) return( nfu_badIndex );
666  if( i1 != i2 ) {
667  for( ; i2 < ptwXY->length; i1++, i2++ ) ptwXY->points[i1] = ptwXY->points[i2];
668  ptwXY->length = n;
669  }
670  return( ptwXY->status );
671 }
ptwXYPoint * points
Definition: ptwXY.h:99
int64_t length
Definition: ptwXY.h:93
const G4int n
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
nfu_status status
Definition: ptwXY.h:85

Here is the call graph for this function:

nfu_status ptwXY_div_doubleFrom ( ptwXYPoints ptwXY,
double  value 
)

Definition at line 44 of file ptwXY_binaryOperators.cc.

44  {
45 
46  if( value == 0. ) {
47  ptwXY->status = nfu_divByZero; }
48  else {
49  ptwXY_slopeOffset( ptwXY, 1. / value, 0. );
50  }
51  return( ptwXY->status );
52 }
const XML_Char int const XML_Char * value
Definition: expat.h:331
nfu_status ptwXY_slopeOffset(ptwXYPoints *ptwXY, double slope, double offset)
nfu_status status
Definition: ptwXY.h:85

Here is the call graph for this function:

nfu_status ptwXY_div_fromDouble ( ptwXYPoints ptwXY,
double  value 
)

Definition at line 53 of file ptwXY_binaryOperators.cc.

53  {
54 /*
55 * This does not do any infilling and it should?????????
56 */
57 
58  int64_t i, nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
59  ptwXYPoint *p;
60  ptwXYOverflowPoint *o, *overflowHeader = &(ptwXY->overflowHeader);
61 
62  if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
64 
65  for( i = 0, p = ptwXY->points; i < nonOverflowLength; i++, p++ ) if( p->y == 0. ) ptwXY->status = nfu_divByZero;
66  for( o = overflowHeader->next; o != overflowHeader; o = o->next ) if( o->point.y == 0. ) ptwXY->status = nfu_divByZero;
67  if( ptwXY->status != nfu_divByZero ) {
68  for( i = 0, p = ptwXY->points; i < nonOverflowLength; i++, p++ ) p->y = value / p->y;
69  for( o = overflowHeader->next; o != overflowHeader; o = o->next ) o->point.y = value / o->point.y;
70  }
71  return( ptwXY->status );
72 }
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
ptwXYPoint point
Definition: ptwXY.h:80
ptwXYPoint * points
Definition: ptwXY.h:99
const char * p
Definition: xmltok.h:285
int64_t ptwXY_getNonOverflowLength(ptwXYPoints const *ptwXY)
Definition: ptwXY_core.cc:590
const XML_Char int const XML_Char * value
Definition: expat.h:331
double y
Definition: ptwXY.h:62
struct ptwXYOverflowPoint_s * next
Definition: ptwXY.h:78
nfu_status status
Definition: ptwXY.h:85
ptwXYOverflowPoint overflowHeader
Definition: ptwXY.h:98

Here is the call graph for this function:

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

Definition at line 288 of file ptwXY_binaryOperators.cc.

288  {
289 
290  int isNAN1, isNAN2;
291  int64_t i, j, k, zeros = 0, length, iYs;
292  double x1, x2, y1, y2, u1, u2, v1, v2, y, xz, nan = nfu_getNAN( ), s1, s2;
293  ptwXYPoints *n = NULL;
294  ptwXYPoint *p;
295 
296  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
297  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY2 ) ) != nfu_Okay ) return( NULL );
298  *status = nfu_otherInterpolation;
299  if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
300  if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( NULL );
301  if( ( ptwXY1->interpolation == ptwXY_interpolationFlat ) || ( ptwXY1->interpolation == ptwXY_interpolationFlat ) )
302  return( ptwXY_div_ptwXY_forFlats( ptwXY1, ptwXY2, status, safeDivide ) );
303 
304  if( ( *status = ptwXY_areDomainsMutual( ptwXY1, ptwXY2 ) ) != nfu_Okay ) return( NULL );
305  if( ( n = ptwXY_union( ptwXY1, ptwXY2, status, ptwXY_union_fill | ptwXY_union_mergeClosePoints ) ) != NULL ) {
306  for( i = 0, p = n->points; i < n->length; i++, p++ ) {
307  if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, p->x, &y ) ) != nfu_Okay ) goto Err;
308  if( y == 0. ) {
309  if( p->y == 0. ) {
310  iYs = 0;
311  y1 = 0.;
312  y2 = 0.;
313  if( i > 0 ) {
314  if( ( *status = ptwXY_getSlopeAtX( ptwXY1, p->x, '-', &s1 ) ) != nfu_Okay ) {
315  if( *status != nfu_XOutsideDomain ) goto Err;
316  s1 = 0.;
317  }
318  if( ( *status = ptwXY_getSlopeAtX( ptwXY2, p->x, '-', &s2 ) ) != nfu_Okay ) goto Err;
319  if( s2 == 0. ) {
320  y1 = nan; }
321  else {
322  y1 = s1 / s2;
323  }
324  iYs++;
325  }
326  if( i < ( n->length - 1 ) ) {
327  if( ( *status = ptwXY_getSlopeAtX( ptwXY1, p->x, '+', &s1 ) ) != nfu_Okay ) {
328  if( *status != nfu_XOutsideDomain ) goto Err;
329  s1 = 0.;
330  }
331  if( ( *status = ptwXY_getSlopeAtX( ptwXY2, p->x, '+', &s2 ) ) != nfu_Okay ) goto Err;
332  if( s2 == 0. ) {
333  y2 = nan; }
334  else {
335  y2 = s1 / s2;
336  }
337  iYs++;
338  }
339  p->y = ( y1 + y2 ) / iYs;
340  if( nfu_isNAN( p->y ) ) zeros++; }
341  else {
342  if( !safeDivide ) {
343  *status = nfu_divByZero;
344  goto Err;
345  }
346  zeros++;
347  p->y = nan;
348  } }
349  else {
350  p->y /= y;
351  }
352  }
353  length = n->length - 1;
354  if( length > 0 ) {
355  x2 = n->points[length].x;
356  for( i = length - 1; i >= 0; i-- ) { /* Find and add y zeros and NAN not currently in n's. */
357  x1 = n->points[i].x;
358  if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x1, &u1 ) ) != nfu_Okay ) goto Err;
359  if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x2, &u2 ) ) != nfu_Okay ) goto Err;
360  if( ( *status = ptwXY_getValueAtX( ptwXY2, x1, &v1 ) ) != nfu_Okay ) goto Err;
361  if( ( *status = ptwXY_getValueAtX( ptwXY2, x2, &v2 ) ) != nfu_Okay ) goto Err;
362  if( u1 * u2 < 0 ) {
363  xz = ( u1 * x2 - u2 * x1 ) / ( u1 - u2 );
364  if( ( *status = ptwXY_setValueAtX( n, xz, 0. ) ) != nfu_Okay ) goto Err;
365  }
366  if( v1 * v2 < 0 ) {
367  if( !safeDivide ) {
368  *status = nfu_divByZero;
369  goto Err;
370  }
371  zeros++;
372  xz = ( v1 * x2 - v2 * x1 ) / ( v1 - v2 );
373  if( ( *status = ptwXY_setValueAtX( n, xz, nan ) ) != nfu_Okay ) goto Err;
374  }
375  x2 = x1;
376  }
377  if( ( *status = ptwXY_simpleCoalescePoints( n ) ) != nfu_Okay ) goto Err;
378  length = n->length;
379  x2 = n->points[n->length-1].x;
380  y2 = n->points[n->length-1].y;
381  isNAN2 = nfu_isNAN( y2 );
382  for( i = n->length - 2; i >= 0; i-- ) { /* Make interpolation fit accuracy. Work backwards so new points will not mess up loop. */
383  x1 = n->points[i].x;
384  y1 = n->points[i].y;
385  isNAN1 = nfu_isNAN( y1 );
386  if( !isNAN1 || !isNAN2 ) {
387  if( ( *status = ptwXY_div_s_ptwXY( n, ptwXY1, ptwXY2, x1, y1, x2, y2, 0, isNAN1, isNAN2 ) ) != nfu_Okay ) goto Err;
388  }
389  x2 = x1;
390  y2 = y1;
391  isNAN2 = isNAN1;
392  }
393  ptwXY_update_biSectionMax( n, (double) length );
394  if( zeros ) {
395  if( ( *status = ptwXY_simpleCoalescePoints( n ) ) != nfu_Okay ) goto Err;
396  for( i = 0; i < n->length; i++ ) if( !nfu_isNAN( n->points[i].y ) ) break;
397  if( nfu_isNAN( n->points[0].y ) ) { /* Special case for first point. */
398  if( i == n->length ) { /* They are all nan's, what now? */
399  zeros = 0;
400  for( i = 0; i < n->length; i++ ) n->points[i].y = 0.; }
401  else {
402  n->points[0].y = 2. * n->points[i].y;
403  zeros--;
404  }
405  }
406  for( i = n->length - 1; i > 0; i-- ) if( !nfu_isNAN( n->points[i].y ) ) break;
407  if( nfu_isNAN( n->points[n->length - 1].y ) ) { /* Special case for last point. */
408  n->points[n->length - 1].y = 2. * n->points[i].y;
409  zeros--;
410  }
411  if( zeros ) {
412  for( i = 0; i < n->length; i++ ) if( nfu_isNAN( n->points[i].y ) ) break;
413  for( k = i + 1, j = i; k < n->length; k++ ) {
414  if( nfu_isNAN( n->points[k].y ) ) continue;
415  n->points[j] = n->points[k];
416  j++;
417  }
418  n->length = j;
419  }
420  }
421  }
422  }
423  return( n );
424 
425 Err:
426  if( n ) ptwXY_free( n );
427  return( NULL );
428 }
nfu_status ptwXY_getValueAtX(ptwXYPoints *ptwXY, double x, double *y)
Definition: ptwXY_core.cc:844
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
ptwXYPoint * points
Definition: ptwXY.h:99
int nfu_isNAN(double d)
Definition: nf_utilities.cc:61
static ptwXYPoints * ptwXY_div_ptwXY_forFlats(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int safeDivide)
const char * p
Definition: xmltok.h:285
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
Definition: ptwXY_core.cc:876
int64_t length
Definition: ptwXY.h:93
#define ptwXY_union_fill
Definition: ptwXY.h:31
ptwXYPoints * ptwXY_union(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int unionOptions)
double nfu_getNAN(void)
Definition: nf_utilities.cc:54
static nfu_status ptwXY_getValueAtX_ignore_XOutsideDomainError(ptwXYPoints *ptwXY1, double x, double *y)
const G4int n
static nfu_status ptwXY_div_s_ptwXY(ptwXYPoints *n, ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, double x1, double y1, double x2, double y2, int level, int isNAN1, int isNAN2)
nfu_status ptwXY_areDomainsMutual(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2)
#define ptwXY_union_mergeClosePoints
Definition: ptwXY.h:33
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
void ptwXY_update_biSectionMax(ptwXYPoints *ptwXY1, double oldLength)
Definition: ptwXY_misc.cc:31
nfu_status ptwXY_getSlopeAtX(ptwXYPoints *ptwXY, double x, const char side, double *slope)
Definition: ptwXY_core.cc:1139

Here is the call graph for this function:

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

Definition at line 42 of file ptwXY_convenient.cc.

42  {
43 
44 #define minEps 5e-16
45 
46  nfu_status status;
47  double xm, xp, dx, y, x1, y1, x2, y2, sign;
48  ptwXYPoint *p;
49 
50 /* This routine can only be used for linear interpolation for the y-axes since for log interpolation, y cannot be 0.
51 This needs to be fixed and documented. */
52  if( ( status = ptwXY->status ) != nfu_Okay ) return( status );
55 
56  if( ptwXY->length < 2 ) return( nfu_Okay );
57 
58  if( lowerEps != 0. ) {
59  if( std::fabs( lowerEps ) < minEps ) {
60  sign = 1;
61  if( lowerEps < 0. ) sign = -1;
62  lowerEps = sign * minEps;
63  }
64 
65  p = ptwXY_getPointAtIndex_Unsafely( ptwXY, 0 );
66  x1 = p->x;
67  y1 = p->y;
68  p = ptwXY_getPointAtIndex_Unsafely( ptwXY, 1 );
69  x2 = p->x;
70  y2 = p->y;
71 
72  if( y1 != 0. ) {
73  dx = std::fabs( x1 * lowerEps );
74  if( x1 == 0 ) dx = std::fabs( lowerEps );
75  xm = x1 - dx;
76  xp = x1 + dx;
77  if( ( xp + dx ) < x2 ) {
78  if( ( status = ptwXY_getValueAtX( ptwXY, xp, &y ) ) != nfu_Okay ) return( status );
79  if( ( status = ptwXY_setValueAtX( ptwXY, xp, y ) ) != nfu_Okay ) return( status ); }
80  else {
81  xp = x2;
82  y = y2;
83  }
84  if( lowerEps > 0 ) {
85  if( ( status = ptwXY_setValueAtX( ptwXY, x1, 0. ) ) != nfu_Okay ) return( status ); }
86  else {
87  if( ( xm < 0. ) && ( x1 >= 0. ) && positiveXOnly ) {
88  if( ( status = ptwXY_setValueAtX( ptwXY, x1, 0. ) ) != nfu_Okay ) return( status ); }
89  else {
90  if( ( status = ptwXY_setValueAtX( ptwXY, xm, 0. ) ) != nfu_Okay ) return( status );
91  if( ( status = ptwXY_interpolatePoint( ptwXY->interpolation, x1, &y, xm, 0., xp, y ) ) != nfu_Okay ) return( status );
92  if( ( status = ptwXY_setValueAtX( ptwXY, x1, y ) ) != nfu_Okay ) return( status );
93  }
94  }
95  }
96  }
97 
98  if( upperEps != 0. ) {
99  if( std::fabs( upperEps ) < minEps ) {
100  sign = 1;
101  if( upperEps < 0. ) sign = -1;
102  upperEps = sign * minEps;
103  }
104 
105  p = ptwXY_getPointAtIndex_Unsafely( ptwXY, ptwXY->length - 2 );
106  x1 = p->x;
107  y1 = p->y;
108  p = ptwXY_getPointAtIndex_Unsafely( ptwXY, ptwXY->length - 1 );
109  x2 = p->x;
110  y2 = p->y;
111 
112  if( y2 != 0. ) {
113  dx = std::fabs( x2 * upperEps );
114  if( x2 == 0 ) dx = std::fabs( upperEps );
115  xm = x2 - dx;
116  xp = x2 + dx;
117  if( ( xm - dx ) > x1 ) {
118  if( ( status = ptwXY_getValueAtX( ptwXY, xm, &y ) ) != nfu_Okay ) return( status );
119  if( ( status = ptwXY_setValueAtX( ptwXY, xm, y ) ) != nfu_Okay ) return( status ); }
120  else {
121  xm = x1;
122  y = y1;
123  }
124  if( upperEps < 0 ) {
125  if( ( status = ptwXY_setValueAtX( ptwXY, x2, 0. ) ) != nfu_Okay ) return( status ); }
126  else {
127  if( ( status = ptwXY_setValueAtX( ptwXY, xp, 0. ) ) != nfu_Okay ) return( status );
128  if( ( status = ptwXY_interpolatePoint( ptwXY->interpolation, x2, &y, xm, y, xp, 0. ) ) != nfu_Okay ) return( status );
129  if( ( status = ptwXY_setValueAtX( ptwXY, x2, y ) ) != nfu_Okay ) return( status );
130  }
131  }
132  }
133 
134  return( ptwXY->status );
135 
136 #undef minEps
137 }
nfu_status ptwXY_getValueAtX(ptwXYPoints *ptwXY, double x, double *y)
Definition: ptwXY_core.cc:844
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
nfu_status ptwXY_interpolatePoint(ptwXY_interpolation interpolation, double x, double *y, double x1, double y1, double x2, double y2)
const char * p
Definition: xmltok.h:285
ptwXYPoint * ptwXY_getPointAtIndex_Unsafely(ptwXYPoints *ptwXY, int64_t index)
Definition: ptwXY_core.cc:684
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
Definition: ptwXY_core.cc:876
int64_t length
Definition: ptwXY.h:93
#define minEps
enum nfu_status_e nfu_status
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
nfu_status status
Definition: ptwXY.h:85
G4int sign(const T t)

Here is the call graph for this function:

Here is the caller graph for this function:

nfu_status ptwXY_exp ( ptwXYPoints ptwXY,
double  a 
)

Definition at line 43 of file ptwXY_functions.cc.

43  {
44 
45  int64_t i, length;
46  nfu_status status;
47  double x1, y1, z1, x2, y2, z2;
48 
49  length = ptwXY->length;
50  if( length < 1 ) return( ptwXY->status );
53 
54  if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status );
55  x2 = ptwXY->points[length-1].x;
56  y2 = a * ptwXY->points[length-1].y;
57  z2 = ptwXY->points[length-1].y = G4Exp( y2 );
58  for( i = length - 2; i >= 0; i-- ) {
59  x1 = ptwXY->points[i].x;
60  y1 = a * ptwXY->points[i].y;
61  z1 = ptwXY->points[i].y = G4Exp( y1 );
62  if( ( status = ptwXY_exp_s( ptwXY, x1, y1, z1, x2, y2, z2, 0 ) ) != nfu_Okay ) return( status );
63  x2 = x1;
64  y2 = y1;
65  }
66  return( status );
67 }
static nfu_status ptwXY_exp_s(ptwXYPoints *ptwXY, double x1, double y1, double z1, double x2, double y2, double z2, int level)
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
ptwXYPoint * points
Definition: ptwXY.h:99
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
int64_t length
Definition: ptwXY.h:93
enum nfu_status_e nfu_status
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
nfu_status status
Definition: ptwXY.h:85

Here is the call graph for this function:

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.

34  {
35 
36  nfu_status status = nfu_Okay;
37  double r;
38 
39  *value = 0.;
40  switch( interpolation ) {
41  case ptwXY_interpolationLinLin : /* x linear, y linear */
42  *value = 0.5 * ( y1 + y2 ) * ( x2 - x1 );
43  break;
44  case ptwXY_interpolationLinLog : /* x linear, y log */
45  if( ( y1 <= 0. ) || ( y2 <= 0. ) ) {
46  status = nfu_badIntegrationInput; }
47  else {
48  r = y2 / y1;
49  if( std::fabs( r - 1. ) < 1e-4 ) {
50  r = r - 1.;
51  *value = y1 * ( x2 - x1 ) / ( 1. + r * ( -0.5 + r * ( 1. / 3. + r * ( -0.25 + .2 * r ) ) ) ); }
52  else {
53  *value = ( y2 - y1 ) * ( x2 - x1 ) / G4Log( r );
54  }
55  }
56  break;
57  case ptwXY_interpolationLogLin : /* x log, y linear */
58  if( ( x1 <= 0. ) || ( x2 <= 0. ) ) {
59  status = nfu_badIntegrationInput; }
60  else {
61  r = x2 / x1;
62  if( std::fabs( r - 1. ) < 1e-4 ) {
63  r = r - 1.;
64  r = r * ( -0.5 + r * ( 1. / 3. + r * ( -0.25 + .2 * r ) ) );
65  *value = x1 * ( y2 - y1 ) * r / ( 1. + r ) + y2 * ( x2 - x1 ); }
66  else {
67  *value = ( y1 - y2 ) * ( x2 - x1 ) / G4Log( r ) + x2 * y2 - x1 * y1;
68  }
69  }
70  break;
71  case ptwXY_interpolationLogLog : /* x log, y log */
72  if( ( x1 <= 0. ) || ( x2 <= 0. ) || ( y1 <= 0. ) || ( y2 <= 0. ) ) {
73  status = nfu_badIntegrationInput; }
74  else {
75  int i, n;
76  double a, z, lx, ly, s, f;
77 
78  r = y2 / y1;
79  if( std::fabs( r - 1. ) < 1e-4 ) {
80  ly = ( y2 - y1 ) / y1;
81  ly = ly * ( 1. + ly * ( -0.5 + ly * ( 1. / 3. - 0.25 * ly ) ) ); }
82  else {
83  ly = G4Log( r );
84  }
85  r = x2 / x1;
86  if( std::fabs( r - 1. ) < 1e-4 ) {
87  lx = ( x2 - x1 ) / x1;
88  lx = lx * ( 1 + lx * ( -0.5 + lx * ( 1. / 3. - 0.25 * lx ) ) ); }
89  else {
90  lx = G4Log( r );
91  }
92  a = ly / lx;
93  if( std::fabs( r - 1. ) < 1e-3 ) {
94  z = ( x2 - x1 ) / x1;
95  n = (int) a;
96  if( n > 10 ) n = 12;
97  if( n < 4 ) n = 6;
98  a = a - n + 1;
99  f = n + 1.;
100  for( i = 0, s = 0.; i < n; i++, a++, f-- ) s = ( 1. + s ) * a * z / f;
101  *value = y1 * ( x2 - x1 ) * ( 1. + s ); }
102  else {
103  *value = y1 * x1 * ( G4Pow::GetInstance()->powA( r, a + 1. ) - 1. ) / ( a + 1. );
104  }
105  }
106  break;
107  case ptwXY_interpolationFlat : /* x ?, y flat */
108  *value = y1 * ( x2 - x1 );
109  break;
111  status = nfu_otherInterpolation;
112  }
113  return( status );
114 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
const XML_Char * s
Definition: expat.h:262
const XML_Char int const XML_Char * value
Definition: expat.h:331
enum nfu_status_e nfu_status
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
const G4int n
G4double G4Log(G4double x)
Definition: G4Log.hh:230
tuple z
Definition: test.py:28

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 74 of file ptwXY_interpolation.cc.

74  {
75 
76  int64_t i, length;
77  double x;
78  ptwXYPoints *n;
79  ptwXYPoint *p1 = NULL, *p2 = NULL, *p3;
80 
81 #define minEps 5e-16
82 
83  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( NULL );
84  *status = nfu_invalidInterpolation;
85  if( ptwXY->interpolation != ptwXY_interpolationFlat ) return( NULL );
86  *status = nfu_badInput;
87  if( ( lowerEps < 0 ) || ( upperEps < 0 ) || ( ( lowerEps == 0 ) && ( upperEps == 0 ) ) ) return( NULL );
88  if( ( lowerEps != 0 ) && ( lowerEps < minEps ) ) lowerEps = minEps;
89  if( ( upperEps != 0 ) && ( upperEps < minEps ) ) upperEps = minEps;
90 
91  length = ptwXY->length * ( 1 + ( lowerEps == 0 ? 0 : 1 ) + ( lowerEps == 0 ? 0 : 1 ) );
92  if( ( n = ptwXY_new( ptwXY_interpolationLinLin, NULL, ptwXY->biSectionMax, ptwXY->accuracy, length, ptwXY->overflowLength, status, ptwXY->userFlag ) ) == NULL ) return( NULL );
93 
94  p3 = ptwXY->points;
95  if( ptwXY->length > 0 ) ptwXY_setValueAtX( n, p3->x, p3->y );
96  for( i = 0; i < ptwXY->length; i++, p3++ ) {
97  if( i > 1 ) {
98  if( lowerEps > 0 ) {
99  x = ptwXY_flatInterpolationToLinear_eps( p2->x, -lowerEps );
100  if( x > p1->x ) {
101  if( ( *status = ptwXY_setValueAtX( n, x, p1->y ) ) != nfu_Okay ) goto Err;
102  }
103  }
104  if( lowerEps == 0 ) if( ( *status = ptwXY_setValueAtX( n, p2->x, p1->y ) ) != nfu_Okay ) goto Err;
105  if( upperEps == 0 ) if( ( *status = ptwXY_setValueAtX( n, p2->x, p2->y ) ) != nfu_Okay ) goto Err;
106  if( upperEps > 0 ) {
107  x = ptwXY_flatInterpolationToLinear_eps( p2->x, upperEps );
108  if( x < p3->x ) {
109  if( ( *status = ptwXY_setValueAtX( n, x, p2->y ) ) != nfu_Okay ) goto Err;
110  }
111  }
112  }
113  p1 = p2;
114  p2 = p3;
115  }
116  if( ptwXY->length > 1 ) {
117  if( ( lowerEps != 0 ) && ( p1->y != p2->y ) ) {
118  x = ptwXY_flatInterpolationToLinear_eps( p2->x, -lowerEps );
119  if( x > p1->x ) {
120  if( ( *status = ptwXY_setValueAtX( n, x, p1->y ) ) != nfu_Okay ) goto Err;
121  }
122  }
123  if( ( *status = ptwXY_setValueAtX( n, p2->x, p2->y ) ) != nfu_Okay ) goto Err;
124  }
125 
126  return( n );
127 
128 Err:
129  ptwXY_free( n );
130  return( NULL );
131 
132 #undef minEps
133 }
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
int userFlag
Definition: ptwXY.h:89
ptwXYPoint * points
Definition: ptwXY.h:99
tuple x
Definition: test.py:50
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
Definition: ptwXY_core.cc:876
double biSectionMax
Definition: ptwXY.h:90
static double ptwXY_flatInterpolationToLinear_eps(double px, double eps)
int64_t length
Definition: ptwXY.h:93
const G4int n
#define minEps
double y
Definition: ptwXY.h:62
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: ptwXY_core.cc:29
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
int64_t overflowLength
Definition: ptwXY.h:95
double accuracy
Definition: ptwXY.h:91

Here is the call graph for this function:

ptwXYPoints* ptwXY_free ( ptwXYPoints ptwXY)

Definition at line 574 of file ptwXY_core.cc.

574  {
575 
576  if( ptwXY != NULL ) ptwXY_release( ptwXY );
577  nfu_free( ptwXY );
578  return( (ptwXYPoints *) NULL );
579 }
nfu_status ptwXY_release(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:549
void * nfu_free(void *p)

Here is the call graph for this function:

Here is the caller graph for this function:

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.

237  {
238 
239  int64_t numberConverted;
240  double *doublePtr;
241  ptwXYPoints *ptwXY = NULL;
242 
243  if( ( *status = nfu_stringToListOfDoubles( str, &numberConverted, &doublePtr, endCharacter ) ) != nfu_Okay ) return( NULL );
244  *status = nfu_oddNumberOfValues;
245  if( ( numberConverted % 2 ) == 0 )
246  ptwXY = ptwXY_create( interpolation, interpolationOtherInfo, biSectionMax, accuracy, numberConverted, 10, numberConverted / 2, doublePtr, status, 0 );
247  nfu_free( doublePtr );
248  return( ptwXY );
249 }
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: ptwXY_core.cc:108
void * nfu_free(void *p)
nfu_status nfu_stringToListOfDoubles(char const *str, int64_t *numberConverted, double **doublePtr, char **endCharacter)

Here is the call graph for this function:

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

Definition at line 331 of file ptwXY_interpolation.cc.

331  {
332 
333  int64_t i, length;
334  ptwXYPoints *n;
335  ptwXYPoint *p, *p2;
336  double dx, inverseDx, xLast = 0.;
337 
338  *status = nfu_tooFewPoints;
339  if( ptwXY->length < 2 ) return( NULL );
340  if( ( n = ptwXY_clone( ptwXY, status ) ) == NULL ) return( NULL );
341 
342  dx = xMax - xMin;
343  inverseDx = 1. / dx;
344  length = n->length;
345  for( i = 0, p2 = p = n->points; i < length; ++i, ++p ) {
346  p2->x = p->x * dx + xMin;
347  if( i > 0 ) {
348  if( std::fabs( p2->x - xLast ) <= 10. * DBL_EPSILON * ( std::fabs( p2->x ) + std::fabs( xLast ) ) ) {
349  --(n->length);
350  continue;
351  }
352  }
353  p2->y = p->y * inverseDx;
354  xLast = p2->x;
355  ++p2;
356  }
357  n->points[n->length-1].x = xMax; /* Make sure last point is realy xMax. */
358  return( n );
359 }
ptwXYPoint * points
Definition: ptwXY.h:99
const char * p
Definition: xmltok.h:285
ptwXYPoints * ptwXY_clone(ptwXYPoints *ptwXY, nfu_status *status)
Definition: ptwXY_core.cc:208
int64_t length
Definition: ptwXY.h:93
#define DBL_EPSILON
Definition: templates.hh:87
const G4int n
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62

Here is the call graph for this function:

Here is the caller graph for this function:

double ptwXY_getAccuracy ( ptwXYPoints ptwXY)

Definition at line 372 of file ptwXY_core.cc.

372  {
373 
374  return( ptwXY->accuracy );
375 }
double accuracy
Definition: ptwXY.h:91
double ptwXY_getBiSectionMax ( ptwXYPoints ptwXY)

Definition at line 390 of file ptwXY_core.cc.

390  {
391 
392  return( ptwXY->biSectionMax );
393 }
double biSectionMax
Definition: ptwXY.h:90
ptwXY_interpolation ptwXY_getInterpolation ( ptwXYPoints ptwXY)

Definition at line 337 of file ptwXY_core.cc.

337  {
338 
339  return( ptwXY->interpolation );
340 }
ptwXY_interpolation interpolation
Definition: ptwXY.h:87

Here is the caller graph for this function:

char const* ptwXY_getInterpolationString ( ptwXYPoints ptwXY)

Definition at line 344 of file ptwXY_core.cc.

344  {
345 
347 }
ptwXY_interpolationOtherInfo interpolationOtherInfo
Definition: ptwXY.h:88
char const * interpolationString
Definition: ptwXY.h:70
int64_t ptwXY_getNonOverflowLength ( ptwXYPoints const *  ptwXY)

Definition at line 590 of file ptwXY_core.cc.

590  {
591 
592  return( ptwXY->length - ptwXY->overflowLength );
593 }

Here is the caller graph for this function:

ptwXYPoint* ptwXY_getPointAtIndex ( ptwXYPoints ptwXY,
int64_t  index 
)

Definition at line 675 of file ptwXY_core.cc.

675  {
676 
677  if( ptwXY->status != nfu_Okay ) return( NULL );
678  if( ( index < 0 ) || ( index >= ptwXY->length ) ) return( NULL );
679  return( ptwXY_getPointAtIndex_Unsafely( ptwXY, index ) );
680 }
ptwXYPoint * ptwXY_getPointAtIndex_Unsafely(ptwXYPoints *ptwXY, int64_t index)
Definition: ptwXY_core.cc:684
int64_t length
Definition: ptwXY.h:93
nfu_status status
Definition: ptwXY.h:85

Here is the call graph for this function:

Here is the caller graph for this function:

ptwXYPoint* ptwXY_getPointAtIndex_Unsafely ( ptwXYPoints ptwXY,
int64_t  index 
)

Definition at line 684 of file ptwXY_core.cc.

684  {
685 
686  int64_t i;
687  ptwXYOverflowPoint *overflowPoint;
688 
689  for( overflowPoint = ptwXY->overflowHeader.next, i = 0; overflowPoint != &(ptwXY->overflowHeader); overflowPoint = overflowPoint->next, i++ ) {
690  if( overflowPoint->index == index ) return( &(overflowPoint->point) );
691  if( overflowPoint->index > index ) break;
692  }
693  return( &(ptwXY->points[index - i]) );
694 }
ptwXYPoint point
Definition: ptwXY.h:80
ptwXYPoint * points
Definition: ptwXY.h:99
struct ptwXYOverflowPoint_s * next
Definition: ptwXY.h:78
ptwXYOverflowPoint overflowHeader
Definition: ptwXY.h:98
int64_t index
Definition: ptwXY.h:79

Here is the caller graph for this function:

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

Definition at line 710 of file ptwXY_core.cc.

710  {
711 
712  int closeIsEqual;
713  ptwXYPoint *closePoint;
714 
715  return( ptwXY_getPointsAroundX_closeIsEqual( ptwXY, x, lessThanEqualXPoint, greaterThanXPoint, 0, &closeIsEqual, &closePoint ) );
716 }
tuple x
Definition: test.py:50
ptwXY_lessEqualGreaterX ptwXY_getPointsAroundX_closeIsEqual(ptwXYPoints *ptwXY, double x, ptwXYOverflowPoint *lessThanEqualXPoint, ptwXYOverflowPoint *greaterThanXPoint, double eps, int *closeIsEqual, ptwXYPoint **closePoint)
Definition: ptwXY_core.cc:720

Here is the call graph for this function:

Here is the caller graph for this function:

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.

721  {
722 
723  int64_t overflowIndex, nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
724  int64_t indexMin, indexMid, indexMax;
725  ptwXY_dataFrom xMinFrom, xMaxFrom;
726  double xMin = ptwXY_getXMinAndFrom( ptwXY, &xMinFrom ), xMax = ptwXY_getXMaxAndFrom( ptwXY, &xMaxFrom );
727  ptwXYOverflowPoint *overflowPoint, *overflowHeader = &(ptwXY->overflowHeader);
729  ptwXYPoint *lowerPoint = NULL, *upperPoint = NULL;
730 
731  ptwXY_initialOverflowPoint( lessThanEqualXPoint, overflowHeader, NULL );
732  ptwXY_initialOverflowPoint( greaterThanXPoint, overflowHeader, NULL );
733  if( ptwXY->length != 0 ) {
734  if( x < xMin ) {
736  if( xMinFrom == ptwXY_dataFrom_Points ) {
737  greaterThanXPoint->prior = overflowHeader;
738  greaterThanXPoint->index = 0;
739  greaterThanXPoint->point = ptwXY->points[0];
740  *closePoint = &(ptwXY->points[0]); }
741  else {
742  *greaterThanXPoint = *(overflowHeader->next);
743  *closePoint = &(overflowHeader->next->point);
744  } }
745  else if( x > xMax ) {
747  if( xMaxFrom == ptwXY_dataFrom_Points ) {
748  lessThanEqualXPoint->prior = overflowHeader->prior;
749  lessThanEqualXPoint->index = nonOverflowLength - 1;
750  lessThanEqualXPoint->point = ptwXY->points[lessThanEqualXPoint->index];
751  *closePoint = &(ptwXY->points[lessThanEqualXPoint->index]); }
752  else {
753  *lessThanEqualXPoint = *(overflowHeader->prior);
754  *closePoint = &(overflowHeader->prior->point);
755  } }
756  else { /* xMin <= x <= xMax */
757  status = ptwXY_lessEqualGreaterX_between; /* Default for this condition, can only be between or equal. */
758  for( overflowPoint = overflowHeader->next, overflowIndex = 0; overflowPoint != overflowHeader;
759  overflowPoint = overflowPoint->next, overflowIndex++ ) if( overflowPoint->point.x > x ) break;
760  overflowPoint = overflowPoint->prior;
761  if( ( overflowPoint != overflowHeader ) && ( overflowPoint->point.x == x ) ) {
763  *lessThanEqualXPoint = *overflowPoint; }
764  else if( ptwXY->length == 1 ) { /* If here and length = 1, then ptwXY->points[0].x == x. */
766  lessThanEqualXPoint->index = 0;
767  lessThanEqualXPoint->point = ptwXY->points[0]; }
768  else { /* ptwXY->length > 1 */
769  indexMin = 0;
770  indexMax = nonOverflowLength - 1;
771  indexMid = ( indexMin + indexMax ) >> 1;
772  while( ( indexMin != indexMid ) && ( indexMid != indexMax ) ) {
773  if( ptwXY->points[indexMid].x > x ) {
774  indexMax = indexMid; }
775  else {
776  indexMin = indexMid;
777  }
778  indexMid = ( indexMin + indexMax ) >> 1;
779  } // Loop checking, 11.06.2015, T. Koi
780  if( ptwXY->points[indexMin].x == x ) {
782  lessThanEqualXPoint->index = indexMin;
783  lessThanEqualXPoint->point = ptwXY->points[indexMin]; }
784  else if( ptwXY->points[indexMax].x == x ) {
786  lessThanEqualXPoint->index = indexMax;
787  lessThanEqualXPoint->point = ptwXY->points[indexMax]; }
788  else {
789  if( ptwXY->points[indexMin].x > x ) indexMax = 0;
790  if( ptwXY->points[indexMax].x < x ) indexMin = indexMax;
791  if( ( overflowPoint == overflowHeader ) || /* x < xMin of overflow points. */
792  ( ( ptwXY->points[indexMin].x > overflowPoint->point.x ) && ( ptwXY->points[indexMin].x < x ) ) ) {
793  if( overflowPoint != overflowHeader ) lessThanEqualXPoint->prior = overflowPoint;
794  lowerPoint = &(ptwXY->points[indexMin]);
795  lessThanEqualXPoint->index = indexMin;
796  lessThanEqualXPoint->point = ptwXY->points[indexMin]; }
797  else {
798  lowerPoint = &(overflowPoint->point);
799  *lessThanEqualXPoint = *overflowPoint;
800  }
801  if( ( overflowPoint->next == overflowHeader ) || /* x > xMax of overflow points. */
802  ( ( ptwXY->points[indexMax].x < overflowPoint->next->point.x ) && ( ptwXY->points[indexMax].x > x ) ) ) {
803  upperPoint = &(ptwXY->points[indexMax]);
804  greaterThanXPoint->index = indexMax;
805  greaterThanXPoint->point = ptwXY->points[indexMax]; }
806  else {
807  upperPoint = &(overflowPoint->next->point);
808  *greaterThanXPoint = *(overflowPoint->next);
809  }
810  }
811  }
812  }
813  }
814 
815  *closeIsEqual = 0;
816  if( eps > 0 ) {
817  double absX = std::fabs( x );
818 
819  if( status == ptwXY_lessEqualGreaterX_lessThan ) {
820  if( absX < std::fabs( greaterThanXPoint->point.x ) ) absX = std::fabs( greaterThanXPoint->point.x );
821  if( ( greaterThanXPoint->point.x - x ) < eps * absX ) *closeIsEqual = 1; }
822  else if( status == ptwXY_lessEqualGreaterX_greater ) {
823  if( absX < std::fabs( lessThanEqualXPoint->point.x ) ) absX = std::fabs( lessThanEqualXPoint->point.x );
824  if( ( x - lessThanEqualXPoint->point.x ) < eps * absX ) *closeIsEqual = -1; }
825  else if( status == ptwXY_lessEqualGreaterX_between ) {
826  if( ( x - lessThanEqualXPoint->point.x ) < ( greaterThanXPoint->point.x - x ) ) { /* x is closer to lower point. */
827  *closePoint = lowerPoint;
828  if( absX < std::fabs( lessThanEqualXPoint->point.x ) ) absX = std::fabs( lessThanEqualXPoint->point.x );
829  if( ( x - lessThanEqualXPoint->point.x ) < eps * absX ) *closeIsEqual = -1; }
830  else { /* x is closer to upper point. */
831  *closePoint = upperPoint;
832  if( absX < std::fabs( greaterThanXPoint->point.x ) ) absX = std::fabs( greaterThanXPoint->point.x );
833  if( ( greaterThanXPoint->point.x - x ) < eps * absX ) *closeIsEqual = 1;
834  } }
835  else if( status == ptwXY_lessEqualGreaterX_equal ) {
836  *closeIsEqual = 1;
837  }
838  }
839  return( status );
840 }
ptwXYPoint point
Definition: ptwXY.h:80
ptwXYPoint * points
Definition: ptwXY.h:99
double ptwXY_getXMaxAndFrom(ptwXYPoints *ptwXY, ptwXY_dataFrom *dataFrom)
Definition: ptwXY_core.cc:1215
static const G4double eps
enum ptwXY_lessEqualGreaterX_e ptwXY_lessEqualGreaterX
tuple x
Definition: test.py:50
struct ptwXYOverflowPoint_s * prior
Definition: ptwXY.h:77
static void ptwXY_initialOverflowPoint(ptwXYOverflowPoint *overflowPoint, ptwXYOverflowPoint *prior, ptwXYOverflowPoint *next)
Definition: ptwXY_core.cc:1290
int64_t length
Definition: ptwXY.h:93
int64_t ptwXY_getNonOverflowLength(ptwXYPoints const *ptwXY)
Definition: ptwXY_core.cc:590
enum ptwXY_dataFrom_e ptwXY_dataFrom
struct ptwXYOverflowPoint_s * next
Definition: ptwXY.h:78
double x
Definition: ptwXY.h:62
double ptwXY_getXMinAndFrom(ptwXYPoints *ptwXY, ptwXY_dataFrom *dataFrom)
Definition: ptwXY_core.cc:1182
ptwXYOverflowPoint overflowHeader
Definition: ptwXY.h:98
int64_t index
Definition: ptwXY.h:79

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 1139 of file ptwXY_core.cc.

1139  {
1140 
1141  nfu_status status = nfu_Okay;
1142  ptwXYOverflowPoint lessThanEqualXPoint, greaterThanXPoint;
1143  ptwXY_lessEqualGreaterX legx = ptwXY_getPointsAroundX( ptwXY, x, &lessThanEqualXPoint, &greaterThanXPoint );
1144  ptwXYPoint *point;
1145 
1146  *slope = 0.;
1147  if( ( side != '-' ) && ( side != '+' ) ) return( nfu_badInput );
1148 
1149  switch( legx ) {
1153  status = nfu_XOutsideDomain;
1154  break;
1156  *slope = ( greaterThanXPoint.point.y - lessThanEqualXPoint.point.y ) /
1157  ( greaterThanXPoint.point.x - lessThanEqualXPoint.point.x );
1158  break;
1160  if( side == '-' ) {
1161  if( lessThanEqualXPoint.index == 0 ) {
1162  status = nfu_XOutsideDomain; }
1163  else {
1164  point = ptwXY_getPointAtIndex_Unsafely( ptwXY, lessThanEqualXPoint.index - 1 );
1165  *slope = ( lessThanEqualXPoint.point.y - point->y ) / ( lessThanEqualXPoint.point.x - point->x );
1166  } }
1167  else {
1168  if( lessThanEqualXPoint.index == ( ptwXY->length - 1 ) ) {
1169  status = nfu_XOutsideDomain; }
1170  else {
1171  point = ptwXY_getPointAtIndex_Unsafely( ptwXY, lessThanEqualXPoint.index + 1 );
1172  *slope = ( point->y - lessThanEqualXPoint.point.y ) / ( point->x - lessThanEqualXPoint.point.x );
1173  }
1174  }
1175  }
1176 
1177  return( status );
1178 }
ptwXYPoint point
Definition: ptwXY.h:80
enum ptwXY_lessEqualGreaterX_e ptwXY_lessEqualGreaterX
tuple x
Definition: test.py:50
ptwXYPoint * ptwXY_getPointAtIndex_Unsafely(ptwXYPoints *ptwXY, int64_t index)
Definition: ptwXY_core.cc:684
int64_t length
Definition: ptwXY.h:93
enum nfu_status_e nfu_status
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
int64_t index
Definition: ptwXY.h:79
ptwXY_lessEqualGreaterX ptwXY_getPointsAroundX(ptwXYPoints *ptwXY, double x, ptwXYOverflowPoint *lessThanEqualXPoint, ptwXYOverflowPoint *greaterThanXPoint)
Definition: ptwXY_core.cc:710

Here is the call graph for this function:

Here is the caller graph for this function:

nfu_status ptwXY_getStatus ( ptwXYPoints ptwXY)

Definition at line 351 of file ptwXY_core.cc.

351  {
352 
353  return( ptwXY->status );
354 }
nfu_status status
Definition: ptwXY.h:85

Here is the caller graph for this function:

int ptwXY_getUserFlag ( ptwXYPoints ptwXY)

Definition at line 358 of file ptwXY_core.cc.

358  {
359 
360  return( ptwXY->userFlag );
361 }
int userFlag
Definition: ptwXY.h:89
nfu_status ptwXY_getValueAtX ( ptwXYPoints ptwXY,
double  x,
double *  y 
)

Definition at line 844 of file ptwXY_core.cc.

844  {
845 
847  ptwXYOverflowPoint lessThanEqualXPoint, greaterThanXPoint;
848  ptwXY_lessEqualGreaterX legx = ptwXY_getPointsAroundX( ptwXY, x, &lessThanEqualXPoint, &greaterThanXPoint );
849 
850  *y = 0.;
851  if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
852  switch( legx ) {
856  break;
858  status = nfu_Okay;
859  *y = lessThanEqualXPoint.point.y;
860  break;
862  if( ptwXY->interpolationOtherInfo.getValueFunc != NULL ) {
864  lessThanEqualXPoint.point.x, lessThanEqualXPoint.point.y, greaterThanXPoint.point.x, greaterThanXPoint.point.y ); }
865  else {
866  status = ptwXY_interpolatePoint( ptwXY->interpolation, x, y, lessThanEqualXPoint.point.x, lessThanEqualXPoint.point.y,
867  greaterThanXPoint.point.x, greaterThanXPoint.point.y );
868  }
869  break;
870  }
871  return( status );
872 }
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
ptwXYPoint point
Definition: ptwXY.h:80
nfu_status ptwXY_interpolatePoint(ptwXY_interpolation interpolation, double x, double *y, double x1, double y1, double x2, double y2)
ptwXY_getValue_callback getValueFunc
Definition: ptwXY.h:71
enum ptwXY_lessEqualGreaterX_e ptwXY_lessEqualGreaterX
tuple x
Definition: test.py:50
enum nfu_status_e nfu_status
ptwXY_interpolationOtherInfo interpolationOtherInfo
Definition: ptwXY.h:88
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
nfu_status status
Definition: ptwXY.h:85
ptwXY_lessEqualGreaterX ptwXY_getPointsAroundX(ptwXYPoints *ptwXY, double x, ptwXYOverflowPoint *lessThanEqualXPoint, ptwXYOverflowPoint *greaterThanXPoint)
Definition: ptwXY_core.cc:710

Here is the call graph for this function:

Here is the caller graph for this function:

ptwXPoints* ptwXY_getXArray ( ptwXYPoints ptwXY,
nfu_status status 
)

Definition at line 24 of file ptwXY_convenient.cc.

24  {
25 
26  int64_t i, n;
27  ptwXPoints *xArray;
28 
29  if( ( *status = ptwXY->status ) != nfu_Okay ) return( NULL );
30  n = ptwXY->length;
31 
32  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( NULL );
33  if( ( xArray = ptwX_new( n, status ) ) == NULL ) return( NULL );
34  for( i = 0; i < n; i++ ) xArray->points[i] = ptwXY->points[i].x;
35  xArray->length = n;
36 
37  return( xArray );
38 }
ptwXYPoint * points
Definition: ptwXY.h:99
int64_t length
Definition: ptwX.h:26
int64_t length
Definition: ptwXY.h:93
ptwXPoints * ptwX_new(int64_t size, nfu_status *status)
Definition: ptwX_core.cc:24
const G4int n
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
nfu_status status
Definition: ptwXY.h:85
double * points
Definition: ptwX.h:29

Here is the call graph for this function:

double ptwXY_getXMax ( ptwXYPoints ptwXY)

Definition at line 1239 of file ptwXY_core.cc.

1239  {
1240 
1241  ptwXY_dataFrom dataFrom;
1242 
1243  return( ptwXY_getXMaxAndFrom( ptwXY, &dataFrom ) );
1244 }
double ptwXY_getXMaxAndFrom(ptwXYPoints *ptwXY, ptwXY_dataFrom *dataFrom)
Definition: ptwXY_core.cc:1215
enum ptwXY_dataFrom_e ptwXY_dataFrom

Here is the call graph for this function:

Here is the caller graph for this function:

double ptwXY_getXMaxAndFrom ( ptwXYPoints ptwXY,
ptwXY_dataFrom dataFrom 
)

Definition at line 1215 of file ptwXY_core.cc.

1215  {
1216 
1217  int64_t nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
1218  double xMax = nfu_getNAN( );
1219 
1220  *dataFrom = ptwXY_dataFrom_Unknown;
1221  if( ptwXY->overflowLength > 0 ) {
1222  *dataFrom = ptwXY_dataFrom_Overflow;
1223  xMax = ptwXY->overflowHeader.prior->point.x;
1224  if( ( nonOverflowLength > 0 ) ) {
1225  if( xMax < ptwXY->points[nonOverflowLength-1].x ) {
1226  *dataFrom = ptwXY_dataFrom_Points;
1227  xMax = ptwXY->points[nonOverflowLength-1].x;
1228  }
1229  } }
1230  else if( ptwXY->length > 0 ) {
1231  *dataFrom = ptwXY_dataFrom_Points;
1232  xMax = ptwXY->points[nonOverflowLength-1].x;
1233  }
1234  return( xMax );
1235 }
ptwXYPoint point
Definition: ptwXY.h:80
ptwXYPoint * points
Definition: ptwXY.h:99
tuple x
Definition: test.py:50
struct ptwXYOverflowPoint_s * prior
Definition: ptwXY.h:77
int64_t length
Definition: ptwXY.h:93
int64_t ptwXY_getNonOverflowLength(ptwXYPoints const *ptwXY)
Definition: ptwXY_core.cc:590
double nfu_getNAN(void)
Definition: nf_utilities.cc:54
double x
Definition: ptwXY.h:62
int64_t overflowLength
Definition: ptwXY.h:95
ptwXYOverflowPoint overflowHeader
Definition: ptwXY.h:98

Here is the call graph for this function:

Here is the caller graph for this function:

double ptwXY_getXMin ( ptwXYPoints ptwXY)

Definition at line 1206 of file ptwXY_core.cc.

1206  {
1207 
1208  ptwXY_dataFrom dataFrom;
1209 
1210  return( ptwXY_getXMinAndFrom( ptwXY, &dataFrom ) );
1211 }
enum ptwXY_dataFrom_e ptwXY_dataFrom
double ptwXY_getXMinAndFrom(ptwXYPoints *ptwXY, ptwXY_dataFrom *dataFrom)
Definition: ptwXY_core.cc:1182

Here is the call graph for this function:

Here is the caller graph for this function:

double ptwXY_getXMinAndFrom ( ptwXYPoints ptwXY,
ptwXY_dataFrom dataFrom 
)

Definition at line 1182 of file ptwXY_core.cc.

1182  {
1183 
1184  int64_t nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
1185  double xMin = nfu_getNAN( );
1186 
1187  *dataFrom = ptwXY_dataFrom_Unknown;
1188  if( ptwXY->overflowLength > 0 ) {
1189  *dataFrom = ptwXY_dataFrom_Overflow;
1190  xMin = ptwXY->overflowHeader.next->point.x;
1191  if( nonOverflowLength >= 0 ) {
1192  if( xMin > ptwXY->points[0].x ) {
1193  *dataFrom = ptwXY_dataFrom_Points;
1194  xMin = ptwXY->points[0].x;
1195  }
1196  } }
1197  else if( nonOverflowLength > 0 ) {
1198  *dataFrom = ptwXY_dataFrom_Points;
1199  xMin = ptwXY->points[0].x;
1200  }
1201  return( xMin );
1202 }
ptwXYPoint point
Definition: ptwXY.h:80
ptwXYPoint * points
Definition: ptwXY.h:99
int64_t ptwXY_getNonOverflowLength(ptwXYPoints const *ptwXY)
Definition: ptwXY_core.cc:590
double nfu_getNAN(void)
Definition: nf_utilities.cc:54
struct ptwXYOverflowPoint_s * next
Definition: ptwXY.h:78
double x
Definition: ptwXY.h:62
int64_t overflowLength
Definition: ptwXY.h:95
ptwXYOverflowPoint overflowHeader
Definition: ptwXY.h:98

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 698 of file ptwXY_core.cc.

698  {
699 
700  ptwXYPoint *p = ptwXY_getPointAtIndex( ptwXY, index );
701 
702  if( p == NULL ) return( nfu_badIndex );
703  *x = p->x;
704  *y = p->y;
705  return( nfu_Okay );
706 }
ptwXYPoint * ptwXY_getPointAtIndex(ptwXYPoints *ptwXY, int64_t index)
Definition: ptwXY_core.cc:675
const char * p
Definition: xmltok.h:285
tuple x
Definition: test.py:50
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62

Here is the call graph for this function:

Here is the caller graph for this function:

double ptwXY_getYMax ( ptwXYPoints ptwXY)

Definition at line 1269 of file ptwXY_core.cc.

1269  {
1270 
1271  int64_t i, n = ptwXY_getNonOverflowLength( ptwXY );
1272  ptwXYPoint *p = ptwXY->points;
1273  ptwXYOverflowPoint *overflowPoint = ptwXY->overflowHeader.next;
1274  double yMax;
1275 
1276  if( ptwXY->length == 0 ) return( 0. );
1277  if( n > 0 ) {
1278  yMax = p->y;
1279  for( i = 1, p++; i < n; i++, p++ ) yMax = ( ( yMax > p->y ) ? yMax : p->y ); }
1280  else {
1281  yMax = overflowPoint->point.y;
1282  }
1283  for( ; overflowPoint != &(ptwXY->overflowHeader); overflowPoint = overflowPoint->next )
1284  yMax = ( ( yMax > overflowPoint->point.y ) ? yMax : overflowPoint->point.y );
1285  return( yMax );
1286 }
ptwXYPoint point
Definition: ptwXY.h:80
ptwXYPoint * points
Definition: ptwXY.h:99
const char * p
Definition: xmltok.h:285
int64_t length
Definition: ptwXY.h:93
int64_t ptwXY_getNonOverflowLength(ptwXYPoints const *ptwXY)
Definition: ptwXY_core.cc:590
const G4int n
double y
Definition: ptwXY.h:62
struct ptwXYOverflowPoint_s * next
Definition: ptwXY.h:78
ptwXYOverflowPoint overflowHeader
Definition: ptwXY.h:98

Here is the call graph for this function:

Here is the caller graph for this function:

double ptwXY_getYMin ( ptwXYPoints ptwXY)

Definition at line 1248 of file ptwXY_core.cc.

1248  {
1249 
1250  int64_t i, n = ptwXY_getNonOverflowLength( ptwXY );
1251  ptwXYPoint *p = ptwXY->points;
1252  ptwXYOverflowPoint *overflowPoint = ptwXY->overflowHeader.next;
1253  double yMin;
1254 
1255  if( ptwXY->length == 0 ) return( 0. );
1256  if( n > 0 ) {
1257  yMin = p->y;
1258  for( i = 1, p++; i < n; i++, p++ ) yMin = ( ( yMin < p->y ) ? yMin : p->y ); }
1259  else {
1260  yMin = overflowPoint->point.y;
1261  }
1262  for( ; overflowPoint != &(ptwXY->overflowHeader); overflowPoint = overflowPoint->next )
1263  yMin = ( ( yMin < overflowPoint->point.y ) ? yMin : overflowPoint->point.y );
1264  return( yMin );
1265 }
ptwXYPoint point
Definition: ptwXY.h:80
ptwXYPoint * points
Definition: ptwXY.h:99
const char * p
Definition: xmltok.h:285
int64_t length
Definition: ptwXY.h:93
int64_t ptwXY_getNonOverflowLength(ptwXYPoints const *ptwXY)
Definition: ptwXY_core.cc:590
const G4int n
double y
Definition: ptwXY.h:62
struct ptwXYOverflowPoint_s * next
Definition: ptwXY.h:78
ptwXYOverflowPoint overflowHeader
Definition: ptwXY.h:98

Here is the call graph for this function:

Here is the caller graph for this function:

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.

363  {
364 
365  int64_t i, igs, ngs;
366  double x1, y1, x2, y2, y2p, xg1, xg2, sum;
367  ptwXYPoints *f;
368  ptwXPoints *groupedData = NULL;
369 
370  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( NULL );
371  if( ( *status = groupBoundaries->status ) != nfu_Okay ) return( NULL );
372  *status = nfu_otherInterpolation;
373  if( ptwXY->interpolation == ptwXY_interpolationOther ) return( NULL );
374 
375  ngs = ptwX_length( groupBoundaries ) - 1;
376  if( normType == ptwXY_group_normType_norm ) {
377  if( ptwX_norm == NULL ) {
378  *status = nfu_badNorm;
379  return( NULL );
380  }
381  *status = ptwX_norm->status;
382  if( ptwX_norm->status != nfu_Okay ) return( NULL );
383  if( ptwX_length( ptwX_norm ) != ngs ) {
384  *status = nfu_badNorm;
385  return( NULL );
386  }
387  }
388 
389  if( ( f = ptwXY_intersectionWith_ptwX( ptwXY, groupBoundaries, status ) ) == NULL ) return( NULL );
390  if( f->length == 0 ) return( ptwX_createLine( ngs, ngs, 0, 0, status ) );
391 
392  if( ( groupedData = ptwX_new( ngs, status ) ) == NULL ) goto err;
393  xg1 = groupBoundaries->points[0];
394  x1 = f->points[0].x;
395  y1 = f->points[0].y;
396  for( igs = 0, i = 1; igs < ngs; igs++ ) {
397  xg2 = groupBoundaries->points[igs+1];
398  sum = 0;
399  if( xg2 > x1 ) {
400  for( ; i < f->length; i++, x1 = x2, y1 = y2 ) {
401  x2 = f->points[i].x;
402  if( x2 > xg2 ) break;
403  y2p = y2 = f->points[i].y;
404  if( f->interpolation == ptwXY_interpolationFlat ) y2p = y1;
405  sum += ( y1 + y2p ) * ( x2 - x1 );
406  }
407  }
408  if( sum != 0. ) {
409  if( normType == ptwXY_group_normType_dx ) {
410  sum /= ( xg2 - xg1 ); }
411  else if( normType == ptwXY_group_normType_norm ) {
412  if( ptwX_norm->points[igs] == 0. ) {
413  *status = nfu_divByZero;
414  goto err;
415  }
416  sum /= ptwX_norm->points[igs];
417  }
418  }
419  groupedData->points[igs] = 0.5 * sum;
420  groupedData->length++;
421  xg1 = xg2;
422  }
423 
424  ptwXY_free( f );
425  return( groupedData );
426 
427 err:
428  ptwXY_free( f );
429  if( groupedData != NULL ) ptwX_free( groupedData );
430  return( NULL );
431 }
ptwXYPoints * ptwXY_intersectionWith_ptwX(ptwXYPoints *ptwXY, ptwXPoints *ptwX, nfu_status *status)
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
ptwXPoints * ptwX_createLine(int64_t size, int64_t length, double slope, double offset, nfu_status *status)
Definition: ptwX_core.cc:62
ptwXYPoint * points
Definition: ptwXY.h:99
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
int64_t length
Definition: ptwX.h:26
int64_t length
Definition: ptwXY.h:93
ptwXPoints * ptwX_new(int64_t size, nfu_status *status)
Definition: ptwX_core.cc:24
ptwXPoints * ptwX_free(ptwXPoints *ptwX)
Definition: ptwX_core.cc:158
double y
Definition: ptwXY.h:62
int64_t ptwX_length(ptwXPoints *ptwX)
Definition: ptwX_core.cc:166
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
nfu_status status
Definition: ptwX.h:25
double * points
Definition: ptwX.h:29

Here is the call graph for this function:

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.

529  {
530 
531  int64_t i, igs, ngs;
532  double x1, fy1, gy1, hy1, x2, fy2, gy2, hy2, fy2p, gy2p, hy2p, xg1, xg2, sum;
533  ptwXYPoints *f = NULL, *ff, *fff = NULL, *g = NULL, *gg = NULL, *h = NULL, *hh = NULL;
534  ptwXPoints *groupedData = NULL;
535 
536  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
537  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY2 ) ) != nfu_Okay ) return( NULL );
538  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY3 ) ) != nfu_Okay ) return( NULL );
539  if( ( *status = groupBoundaries->status ) != nfu_Okay ) return( NULL );
540  *status = nfu_otherInterpolation;
541  if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
542  if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( NULL );
543  if( ptwXY3->interpolation == ptwXY_interpolationOther ) return( NULL );
544 
545  ngs = ptwX_length( groupBoundaries ) - 1;
546  if( normType == ptwXY_group_normType_norm ) {
547  if( ptwX_norm == NULL ) {
548  *status = nfu_badNorm;
549  return( NULL );
550  }
551  if( ( *status = ptwX_norm->status ) != nfu_Okay ) return( NULL );
552  if( ptwX_length( ptwX_norm ) != ngs ) {
553  *status = nfu_badNorm;
554  return( NULL );
555  }
556  }
557 
558  if( ( ff = ptwXY_intersectionWith_ptwX( ptwXY1, groupBoundaries, status ) ) == NULL ) return( NULL );
559  if( ( gg = ptwXY_intersectionWith_ptwX( ptwXY2, groupBoundaries, status ) ) == NULL ) goto err;
560  if( ( hh = ptwXY_intersectionWith_ptwX( ptwXY3, groupBoundaries, status ) ) == NULL ) goto err;
561  if( ( ff->length == 0 ) || ( gg->length == 0 ) || ( hh->length == 0 ) ) return( ptwX_createLine( ngs, ngs, 0, 0, status ) );
562 
563  if( ( *status = ptwXY_tweakDomainsToMutualify( ff, gg, 4, 0 ) ) != nfu_Okay ) goto err;
564  if( ( *status = ptwXY_tweakDomainsToMutualify( ff, hh, 4, 0 ) ) != nfu_Okay ) goto err;
565  if( ( *status = ptwXY_tweakDomainsToMutualify( gg, hh, 4, 0 ) ) != nfu_Okay ) goto err;
566  if( ( fff = ptwXY_union( ff, gg, status, ptwXY_union_fill ) ) == NULL ) goto err;
567  if( ( h = ptwXY_union( hh, fff, status, ptwXY_union_fill ) ) == NULL ) goto err;
568  if( ( f = ptwXY_union( fff, h, status, ptwXY_union_fill ) ) == NULL ) goto err;
569  if( ( g = ptwXY_union( gg, h, status, ptwXY_union_fill ) ) == NULL ) goto err;
570 
571  if( ( groupedData = ptwX_new( ngs, status ) ) == NULL ) goto err;
572  xg1 = groupBoundaries->points[0];
573  x1 = f->points[0].x;
574  fy1 = f->points[0].y;
575  gy1 = g->points[0].y;
576  hy1 = h->points[0].y;
577  for( igs = 0, i = 1; igs < ngs; igs++ ) {
578  xg2 = groupBoundaries->points[igs+1];
579  sum = 0;
580  if( xg2 > x1 ) {
581  for( ; i < f->length; i++, x1 = x2, fy1 = fy2, gy1 = gy2, hy1 = hy2 ) {
582  x2 = f->points[i].x;
583  if( x2 > xg2 ) break;
584  fy2p = fy2 = f->points[i].y;
585  if( f->interpolation == ptwXY_interpolationFlat ) fy2p = fy1;
586  gy2p = gy2 = g->points[i].y;
587  if( g->interpolation == ptwXY_interpolationFlat ) gy2p = gy1;
588  hy2p = hy2 = h->points[i].y;
589  if( h->interpolation == ptwXY_interpolationFlat ) hy2p = hy1;
590  sum += ( ( fy1 + fy2p ) * ( gy1 + gy2p ) * ( hy1 + hy2p ) + 2 * fy1 * gy1 * hy1 + 2 * fy2p * gy2p * hy2p ) * ( x2 - x1 );
591  }
592  }
593  if( sum != 0. ) {
594  if( normType == ptwXY_group_normType_dx ) {
595  sum /= ( xg2 - xg1 ); }
596  else if( normType == ptwXY_group_normType_norm ) {
597  if( ptwX_norm->points[igs] == 0. ) {
598  *status = nfu_divByZero;
599  goto err;
600  }
601  sum /= ptwX_norm->points[igs];
602  }
603  }
604  groupedData->points[igs] = sum / 12.;
605  groupedData->length++;
606  xg1 = xg2;
607  }
608 
609  ptwXY_free( f );
610  ptwXY_free( g );
611  ptwXY_free( h );
612  ptwXY_free( ff );
613  ptwXY_free( gg );
614  ptwXY_free( hh );
615  ptwXY_free( fff );
616  return( groupedData );
617 
618 err:
619  ptwXY_free( ff );
620  if( fff != NULL ) ptwXY_free( fff );
621  if( gg != NULL ) ptwXY_free( gg );
622  if( hh != NULL ) ptwXY_free( hh );
623  if( f != NULL ) ptwXY_free( f );
624  if( g != NULL ) ptwXY_free( g );
625  if( h != NULL ) ptwXY_free( h );
626  if( groupedData != NULL ) ptwX_free( groupedData );
627  return( NULL );
628 }
ptwXYPoints * ptwXY_intersectionWith_ptwX(ptwXYPoints *ptwXY, ptwXPoints *ptwX, nfu_status *status)
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
ptwXPoints * ptwX_createLine(int64_t size, int64_t length, double slope, double offset, nfu_status *status)
Definition: ptwX_core.cc:62
ptwXYPoint * points
Definition: ptwXY.h:99
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
int64_t length
Definition: ptwX.h:26
int64_t length
Definition: ptwXY.h:93
#define ptwXY_union_fill
Definition: ptwXY.h:31
ptwXPoints * ptwX_new(int64_t size, nfu_status *status)
Definition: ptwX_core.cc:24
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
ptwXYPoints * ptwXY_union(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int unionOptions)
ptwXPoints * ptwX_free(ptwXPoints *ptwX)
Definition: ptwX_core.cc:158
double y
Definition: ptwXY.h:62
int64_t ptwX_length(ptwXPoints *ptwX)
Definition: ptwX_core.cc:166
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
nfu_status ptwXY_tweakDomainsToMutualify(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, int epsilonFactor, double epsilon)
nfu_status status
Definition: ptwX.h:25
double * points
Definition: ptwX.h:29

Here is the call graph for this function:

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.

436  {
437 
438  int64_t i, igs, ngs;
439  double x1, fy1, gy1, x2, fy2, gy2, fy2p, gy2p, xg1, xg2, sum;
440  ptwXYPoints *f = NULL, *ff, *g = NULL, *gg = NULL;
441  ptwXPoints *groupedData = NULL;
442 
443  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
444  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY2 ) ) != nfu_Okay ) return( NULL );
445  if( ( *status = groupBoundaries->status ) != nfu_Okay ) return( NULL );
446  *status = nfu_otherInterpolation;
447  if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
448  if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( NULL );
449 
450  ngs = ptwX_length( groupBoundaries ) - 1;
451  if( normType == ptwXY_group_normType_norm ) {
452  if( ptwX_norm == NULL ) {
453  *status = nfu_badNorm;
454  return( NULL );
455  }
456  if( ( *status = ptwX_norm->status ) != nfu_Okay ) return( NULL );
457  if( ptwX_length( ptwX_norm ) != ngs ) {
458  *status = nfu_badNorm;
459  return( NULL );
460  }
461  }
462 
463  if( ( ff = ptwXY_intersectionWith_ptwX( ptwXY1, groupBoundaries, status ) ) == NULL ) return( NULL );
464  if( ( gg = ptwXY_intersectionWith_ptwX( ptwXY2, groupBoundaries, status ) ) == NULL ) goto err;
465  if( ( ff->length == 0 ) || ( gg->length == 0 ) ) {
466  ptwXY_free( ff );
467  ptwXY_free( gg );
468  return( ptwX_createLine( ngs, ngs, 0, 0, status ) );
469  }
470 
471  if( ( *status = ptwXY_tweakDomainsToMutualify( ff, gg, 4, 0 ) ) != nfu_Okay ) goto err;
472  if( ( f = ptwXY_union( ff, gg, status, ptwXY_union_fill ) ) == NULL ) goto err;
473  if( ( g = ptwXY_union( gg, f, status, ptwXY_union_fill ) ) == NULL ) goto err;
474 
475  if( ( groupedData = ptwX_new( ngs, status ) ) == NULL ) goto err;
476  xg1 = groupBoundaries->points[0];
477  x1 = f->points[0].x;
478  fy1 = f->points[0].y;
479  gy1 = g->points[0].y;
480  for( igs = 0, i = 1; igs < ngs; igs++ ) {
481  xg2 = groupBoundaries->points[igs+1];
482  sum = 0;
483  if( xg2 > x1 ) {
484  for( ; i < f->length; i++, x1 = x2, fy1 = fy2, gy1 = gy2 ) {
485  x2 = f->points[i].x;
486  if( x2 > xg2 ) break;
487  fy2p = fy2 = f->points[i].y;
488  if( f->interpolation == ptwXY_interpolationFlat ) fy2p = fy1;
489  gy2p = gy2 = g->points[i].y;
490  if( g->interpolation == ptwXY_interpolationFlat ) gy2p = gy1;
491  sum += ( ( fy1 + fy2p ) * ( gy1 + gy2p ) + fy1 * gy1 + fy2p * gy2p ) * ( x2 - x1 );
492  }
493  }
494  if( sum != 0. ) {
495  if( normType == ptwXY_group_normType_dx ) {
496  sum /= ( xg2 - xg1 ); }
497  else if( normType == ptwXY_group_normType_norm ) {
498  if( ptwX_norm->points[igs] == 0. ) {
499  *status = nfu_divByZero;
500  goto err;
501  }
502  sum /= ptwX_norm->points[igs];
503  }
504  }
505  groupedData->points[igs] = sum / 6.;
506  groupedData->length++;
507  xg1 = xg2;
508  }
509 
510  ptwXY_free( f );
511  ptwXY_free( g );
512  ptwXY_free( ff );
513  ptwXY_free( gg );
514  return( groupedData );
515 
516 err:
517  ptwXY_free( ff );
518  if( gg != NULL ) ptwXY_free( gg );
519  // Coverity #63063
520  if( f != NULL ) ptwXY_free( f );
521  if( g != NULL ) ptwXY_free( g );
522  if( groupedData != NULL ) ptwX_free( groupedData );
523  return( NULL );
524 }
ptwXYPoints * ptwXY_intersectionWith_ptwX(ptwXYPoints *ptwXY, ptwXPoints *ptwX, nfu_status *status)
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
ptwXPoints * ptwX_createLine(int64_t size, int64_t length, double slope, double offset, nfu_status *status)
Definition: ptwX_core.cc:62
ptwXYPoint * points
Definition: ptwXY.h:99
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
int64_t length
Definition: ptwX.h:26
int64_t length
Definition: ptwXY.h:93
#define ptwXY_union_fill
Definition: ptwXY.h:31
ptwXPoints * ptwX_new(int64_t size, nfu_status *status)
Definition: ptwX_core.cc:24
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
ptwXYPoints * ptwXY_union(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int unionOptions)
ptwXPoints * ptwX_free(ptwXPoints *ptwX)
Definition: ptwX_core.cc:158
double y
Definition: ptwXY.h:62
int64_t ptwX_length(ptwXPoints *ptwX)
Definition: ptwX_core.cc:166
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
nfu_status ptwXY_tweakDomainsToMutualify(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, int epsilonFactor, double epsilon)
nfu_status status
Definition: ptwX.h:25
double * points
Definition: ptwX.h:29

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 118 of file ptwXY_integration.cc.

118  {
119 
120  int64_t i, n = ptwXY->length;
121  double sum = 0., dSum, x, y, x1, x2, y1, y2, _sign = 1.;
122  ptwXYPoint *point;
123 
124  if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
125  *status = nfu_otherInterpolation;
126  if( ptwXY->interpolation == ptwXY_interpolationOther ) return( 0. );
127 
128  if( xMax < xMin ) {
129  x = xMin;
130  xMin = xMax;
131  xMax = x;
132  _sign = -1.;
133  }
134  if( n < 2 ) return( 0. );
135 
136  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( 0. );
137  for( i = 0, point = ptwXY->points; i < n; i++, point++ ) {
138  if( point->x >= xMin ) break;
139  }
140  if( i == n ) return( 0. );
141  x2 = point->x;
142  y2 = point->y;
143  if( i > 0 ) {
144  if( x2 > xMin ) {
145  x1 = point[-1].x;
146  y1 = point[-1].y;
147  if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMin, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( 0. );
148  if( x2 > xMax ) {
149  double yMax;
150 
151  if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMax, &yMax, x1, y1, x2, y2 ) ) != nfu_Okay ) return( 0. );
152  if( ( *status = ptwXY_f_integrate( ptwXY->interpolation, xMin, y, xMax, yMax, &sum ) ) != nfu_Okay ) return( 0. );
153  return( sum ); }
154  else {
155  if( ( *status = ptwXY_f_integrate( ptwXY->interpolation, xMin, y, x2, y2, &sum ) ) != nfu_Okay ) return( 0. );
156  }
157  }
158  }
159  i++;
160  point++;
161  for( ; i < n; i++, point++ ) {
162  x1 = x2;
163  y1 = y2;
164  x2 = point->x;
165  y2 = point->y;
166  if( x2 > xMax ) {
167  if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMax, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( 0. );
168  if( ( *status = ptwXY_f_integrate( ptwXY->interpolation, x1, y1, xMax, y, &dSum ) ) != nfu_Okay ) return( 0. );
169  sum += dSum;
170  break;
171  }
172  if( ( *status = ptwXY_f_integrate( ptwXY->interpolation, x1, y1, x2, y2, &dSum ) ) != nfu_Okay ) return( 0. );
173  sum += dSum;
174  }
175 
176  return( _sign * sum );
177 }
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
nfu_status ptwXY_interpolatePoint(ptwXY_interpolation interpolation, double x, double *y, double x1, double y1, double x2, double y2)
ptwXYPoint * points
Definition: ptwXY.h:99
tuple x
Definition: test.py:50
int64_t length
Definition: ptwXY.h:93
nfu_status ptwXY_f_integrate(ptwXY_interpolation interpolation, double x1, double y1, double x2, double y2, double *value)
const G4int n
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
nfu_status status
Definition: ptwXY.h:85

Here is the call graph for this function:

Here is the caller graph for this function:

double ptwXY_integrateDomain ( ptwXYPoints ptwXY,
nfu_status status 
)

Definition at line 181 of file ptwXY_integration.cc.

181  {
182 
183  if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
184  if( ptwXY->length > 0 ) return( ptwXY_integrate( ptwXY, ptwXY_getXMin( ptwXY ), ptwXY_getXMax( ptwXY ), status ) );
185  return( 0. );
186 }
int64_t length
Definition: ptwXY.h:93
double ptwXY_getXMin(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:1206
double ptwXY_integrate(ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status)
nfu_status status
Definition: ptwXY.h:85
double ptwXY_getXMax(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:1239

Here is the call graph for this function:

Here is the caller graph for this function:

double ptwXY_integrateDomainWithWeight_sqrt_x ( ptwXYPoints ptwXY,
nfu_status status 
)

Definition at line 284 of file ptwXY_integration.cc.

284  {
285 
286  if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
287  if( ptwXY->length < 2 ) return( 0. );
288  return( ptwXY_integrateWithWeight_sqrt_x( ptwXY, ptwXY_getXMin( ptwXY ), ptwXY_getXMax( ptwXY ), status ) );
289 }
double ptwXY_integrateWithWeight_sqrt_x(ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status)
int64_t length
Definition: ptwXY.h:93
double ptwXY_getXMin(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:1206
nfu_status status
Definition: ptwXY.h:85
double ptwXY_getXMax(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:1239

Here is the call graph for this function:

double ptwXY_integrateDomainWithWeight_x ( ptwXYPoints ptwXY,
nfu_status status 
)

Definition at line 210 of file ptwXY_integration.cc.

210  {
211 
212  if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
213  if( ptwXY->length < 2 ) return( 0. );
214  return( ptwXY_integrateWithWeight_x( ptwXY, ptwXY_getXMin( ptwXY ), ptwXY_getXMax( ptwXY ), status ) );
215 }
int64_t length
Definition: ptwXY.h:93
double ptwXY_getXMin(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:1206
nfu_status status
Definition: ptwXY.h:85
double ptwXY_integrateWithWeight_x(ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status)
double ptwXY_getXMax(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:1239

Here is the call graph for this function:

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.

658  {
659 
660  int64_t i1, i2, n1 = ptwXY->length;
661  long evaluations;
662  double integral = 0., integral_, sign = -1., xa, xb;
663  ptwXY_integrateWithFunctionInfo integrateWithFunctionInfo;
664  ptwXYPoint *point;
665 
666  if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
667 
668  if( xMin == xMax ) return( 0. );
669  if( n1 < 2 ) return( 0. );
670 
672 
673  if( xMin > xMax ) {
674  sign = xMin;
675  xMin = xMax;
676  xMax = sign;
677  sign = -1.;
678  }
679  if( xMin >= ptwXY->points[n1-1].x ) return( 0. );
680  if( xMax <= ptwXY->points[0].x ) return( 0. );
681 
682  for( i1 = 0; i1 < ( n1 - 1 ); i1++ ) {
683  if( ptwXY->points[i1+1].x > xMin ) break;
684  }
685  for( i2 = n1 - 1; i2 > i1; i2-- ) {
686  if( ptwXY->points[i2-1].x < xMax ) break;
687  }
688  point = &(ptwXY->points[i1]);
689 
690  integrateWithFunctionInfo.degree = degree;
691  integrateWithFunctionInfo.func = func;
692  integrateWithFunctionInfo.argList = argList;
693  integrateWithFunctionInfo.interpolation = ptwXY->interpolation;
694  integrateWithFunctionInfo.x2 = point->x;
695  integrateWithFunctionInfo.y2 = point->y;
696 
697  xa = xMin;
698  for( ; i1 < i2; i1++ ) {
699  integrateWithFunctionInfo.x1 = integrateWithFunctionInfo.x2;
700  integrateWithFunctionInfo.y1 = integrateWithFunctionInfo.y2;
701  ++point;
702  integrateWithFunctionInfo.x2 = point->x;
703  integrateWithFunctionInfo.y2 = point->y;
704  xb = point->x;
705  if( xb > xMax ) xb = xMax;
707  xa, xb, recursionLimit, tolerance, &integral_, &evaluations );
708  if( *status != nfu_Okay ) return( 0. );
709  integral += integral_;
710  xa = xb;
711  }
712 
713  return( integral );
714 }
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
ptwXYPoint * points
Definition: ptwXY.h:99
tuple x
Definition: test.py:50
int64_t length
Definition: ptwXY.h:93
static constexpr double degree
Definition: G4SIunits.hh:144
static nfu_status ptwXY_integrateWithFunction3(double x, double *y, void *argList)
static nfu_status ptwXY_integrateWithFunction2(nf_Legendre_GaussianQuadrature_callback integrandFunction, void *argList, double x1, double x2, double *integral)
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
nfu_status status
Definition: ptwXY.h:85
ptwXY_createFromFunction_callback func
G4int sign(const T t)
nfu_status nf_GnG_adaptiveQuadrature(nf_GnG_adaptiveQuadrature_callback quadratureFunction, nf_Legendre_GaussianQuadrature_callback integrandFunction, void *argList, double x1, double x2, int maxDepth, double tolerance, double *integral, long *evaluations)

Here is the call graph for this function:

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

Definition at line 293 of file ptwXY_integration.cc.

293  {
294 
295  int64_t i, n = ptwXY->length;
296  double sum = 0., x, y, x1, x2, y1, y2, _sign = 1., sqrt_x1, sqrt_x2, inv_apb, c;
297  ptwXYPoint *point;
298 
299  if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
301  if( ( ptwXY->interpolation != ptwXY_interpolationLinLin ) &&
302  ( ptwXY->interpolation != ptwXY_interpolationFlat ) ) return( 0. );
303 
304  if( n < 2 ) return( 0. );
305  if( xMax < xMin ) {
306  x = xMin;
307  xMin = xMax;
308  xMax = x;
309  _sign = -1.;
310  }
311 
312  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( 0. );
313  for( i = 0, point = ptwXY->points; i < n; ++i, ++point ) {
314  if( point->x >= xMin ) break;
315  }
316  if( i == n ) return( 0. );
317  x2 = point->x;
318  y2 = point->y;
319  if( i > 0 ) {
320  if( x2 > xMin ) {
321  if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMin, &y, point[-1].x, point[-1].y, x2, y2 ) ) != nfu_Okay ) return( 0. );
322  x2 = xMin;
323  y2 = y;
324  --i;
325  --point;
326  }
327  }
328  ++i;
329  ++point;
330  sqrt_x2 = std::sqrt( x2 );
331  for( ; i < n; ++i, ++point ) {
332  x1 = x2;
333  y1 = y2;
334  sqrt_x1 = sqrt_x2;
335  x2 = point->x;
336  y2 = point->y;
337  if( x2 > xMax ) {
338  if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMax, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( 0. );
339  x2 = xMax;
340  y2 = y;
341  }
342  sqrt_x2 = std::sqrt( x2 );
343  inv_apb = sqrt_x1 + sqrt_x2;
344  c = 2. * ( sqrt_x1 * sqrt_x2 + x1 + x2 );
345  switch( ptwXY->interpolation ) {
347  sum += ( sqrt_x2 - sqrt_x1 ) * y1 * 2.5 * c;
348  break;
350  sum += ( sqrt_x2 - sqrt_x1 ) * ( y1 * ( c + x1 * ( 1. + sqrt_x2 / inv_apb ) ) + y2 * ( c + x2 * ( 1. + sqrt_x1 / inv_apb ) ) );
351  break;
352  default : /* Only to stop compilers from complaining. */
353  break;
354  }
355  if( x2 == xMax ) break;
356  }
357 
358  return( 2. / 15. * _sign * sum );
359 }
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
nfu_status ptwXY_interpolatePoint(ptwXY_interpolation interpolation, double x, double *y, double x1, double y1, double x2, double y2)
ptwXYPoint * points
Definition: ptwXY.h:99
tuple x
Definition: test.py:50
int64_t length
Definition: ptwXY.h:93
const G4int n
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
nfu_status status
Definition: ptwXY.h:85
tuple c
Definition: test.py:13

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 219 of file ptwXY_integration.cc.

219  {
220 
221  int64_t i, n = ptwXY->length;
222  double sum = 0., x, y, x1, x2, y1, y2, _sign = 1.;
223  ptwXYPoint *point;
224 
225  if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
227  if( ( ptwXY->interpolation != ptwXY_interpolationLinLin ) &&
228  ( ptwXY->interpolation != ptwXY_interpolationFlat ) ) return( 0. );
229 
230  if( n < 2 ) return( 0. );
231  if( xMax < xMin ) {
232  x = xMin;
233  xMin = xMax;
234  xMax = x;
235  _sign = -1.;
236  }
237 
238  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( 0. );
239  for( i = 0, point = ptwXY->points; i < n; ++i, ++point ) {
240  if( point->x >= xMin ) break;
241  }
242  if( i == n ) return( 0. );
243  x2 = point->x;
244  y2 = point->y;
245  if( i > 0 ) {
246  if( x2 > xMin ) {
247  if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMin, &y, point[-1].x, point[-1].y, x2, y2 ) ) != nfu_Okay ) return( 0. );
248  x2 = xMin;
249  y2 = y;
250  --i;
251  --point;
252  }
253  }
254  ++i;
255  ++point;
256  for( ; i < n; ++i, ++point ) {
257  x1 = x2;
258  y1 = y2;
259  x2 = point->x;
260  y2 = point->y;
261  if( x2 > xMax ) {
262  if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMax, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( 0. );
263  x2 = xMax;
264  y2 = y;
265  }
266  switch( ptwXY->interpolation ) {
268  sum += ( x2 - x1 ) * y1 * 3 * ( x1 + x2 );
269  break;
271  sum += ( x2 - x1 ) * ( y1 * ( 2 * x1 + x2 ) + y2 * ( x1 + 2 * x2 ) );
272  break;
273  default : /* Only to stop compilers from complaining. */
274  break;
275  }
276  if( x2 == xMax ) break;
277  }
278 
279  return( _sign * sum / 6 );
280 }
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
nfu_status ptwXY_interpolatePoint(ptwXY_interpolation interpolation, double x, double *y, double x1, double y1, double x2, double y2)
ptwXYPoint * points
Definition: ptwXY.h:99
tuple x
Definition: test.py:50
int64_t length
Definition: ptwXY.h:93
const G4int n
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
nfu_status status
Definition: ptwXY.h:85

Here is the call graph for this function:

Here is the caller graph for this function:

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.

30  {
31 
32  nfu_status status = nfu_Okay;
33 
34  if( interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
35  if( ( x1 > x2 ) || ( x < x1 ) || ( x > x2 ) ) return( nfu_invalidInterpolation );
36  if( y1 == y2 ) {
37  *y = y1; }
38  else if( x1 == x2 ) {
39  *y = 0.5 * ( y1 + y2 ); }
40  else if( x == x1 ) {
41  *y = y1; }
42  else if( x == x2 ) {
43  *y = y2; }
44  else {
45  switch( interpolation ) {
47  *y = ( y1 * ( x2 - x ) + y2 * ( x - x1 ) ) / ( x2 - x1 );
48  break;
50  if( ( x <= 0. ) || ( x1 <= 0. ) || ( x2 <= 0. ) ) return( nfu_invalidInterpolation );
51  *y = ( y1 * G4Log( x2 / x ) + y2 * G4Log( x / x1 ) ) / G4Log( x2 / x1 );
52  break;
54  if( ( y1 <= 0. ) || ( y2 <= 0. ) ) return( nfu_invalidInterpolation );
55  *y = G4Exp( ( G4Log( y1 ) * ( x2 - x ) + G4Log( y2 ) * ( x - x1 ) ) / ( x2 - x1 ) );
56  break;
58  if( ( x <= 0. ) || ( x1 <= 0. ) || ( x2 <= 0. ) ) return( nfu_invalidInterpolation );
59  if( ( y1 <= 0. ) || ( y2 <= 0. ) ) return( nfu_invalidInterpolation );
60  *y = G4Exp( ( G4Log( y1 ) * G4Log( x2 / x ) + G4Log( y2 ) * G4Log( x / x1 ) ) / G4Log( x2 / x1 ) );
61  break;
63  *y = y1;
64  break;
65  default :
66  status = nfu_invalidInterpolation;
67  }
68  }
69  return( status );
70 }
tuple x
Definition: test.py:50
enum nfu_status_e nfu_status
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 194 of file ptwXY_convenient.cc.

194  {
195 
196  int64_t i, i1, i2, lengthX = ptwX_length( ptwX );
197  double x, y, xMin, xMax;
198  ptwXYPoints *n = NULL;
199 
200  if( ( *status = ptwXY->status ) != nfu_Okay ) return( NULL );
201  if( ( *status = ptwX->status ) != nfu_Okay ) return( NULL );
202  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) goto Err;
203  *status = nfu_otherInterpolation;
204  if( ptwXY->interpolation == ptwXY_interpolationOther ) return( NULL );
205  if( ( n = ptwXY_clone( ptwXY, status ) ) == NULL ) return( NULL );
206  if( ptwXY->length == 0 ) return( n );
207  xMin = ptwXY->points[0].x;
208  xMax = ptwXY->points[ptwXY->length - 1].x;
209 
210  if( ( xMin >= ptwX->points[lengthX-1] ) || ( xMax <= ptwX->points[0] ) ) { /* No overlap. */
211  n->length = 0;
212  return( n );
213  }
214 
215  for( i = 0; i < lengthX; i++ ) { /* Fill in ptwXY at x-points in ptwX. */
216  x = ptwX->points[i];
217  if( x <= xMin ) continue;
218  if( x >= xMax ) break;
219  if( ( *status = ptwXY_getValueAtX( ptwXY, x, &y ) ) != nfu_Okay ) goto Err;
220  if( ( *status = ptwXY_setValueAtX( n, x, y ) ) != nfu_Okay ) goto Err;
221  }
222  if( ( *status = ptwXY_simpleCoalescePoints( n ) ) != nfu_Okay ) goto Err;
223 
224  i1 = 0;
225  i2 = n->length - 1;
226  if( lengthX > 0 ) {
227  x = ptwX->points[0];
228  if( x > n->points[i1].x ) {
229  for( ; i1 < n->length; i1++ ) {
230  if( n->points[i1].x == x ) break;
231  }
232  }
233 
234  x = ptwX->points[lengthX - 1];
235  if( x < n->points[i2].x ) {
236  for( ; i2 > i1; i2-- ) {
237  if( n->points[i2].x == x ) break;
238  }
239  }
240  }
241  i2++;
242 
243  if( i1 != 0 ) {
244  for( i = i1; i < i2; i++ ) n->points[i - i1] = n->points[i];
245  }
246  n->length = i2 - i1;
247 
248  return( n );
249 
250 Err:
251  ptwXY_free( n );
252  return( NULL );
253 }
nfu_status ptwXY_getValueAtX(ptwXYPoints *ptwXY, double x, double *y)
Definition: ptwXY_core.cc:844
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
ptwXYPoint * points
Definition: ptwXY.h:99
ptwXYPoints * ptwXY_clone(ptwXYPoints *ptwXY, nfu_status *status)
Definition: ptwXY_core.cc:208
tuple x
Definition: test.py:50
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
Definition: ptwXY_core.cc:876
int64_t length
Definition: ptwXY.h:93
const G4int n
int64_t ptwX_length(ptwXPoints *ptwX)
Definition: ptwX_core.cc:166
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
nfu_status status
Definition: ptwXY.h:85
nfu_status status
Definition: ptwX.h:25
double * points
Definition: ptwX.h:29

Here is the call graph for this function:

Here is the caller graph for this function:

int64_t ptwXY_length ( ptwXYPoints ptwXY)

Definition at line 583 of file ptwXY_core.cc.

583  {
584 
585  return( ptwXY->length );
586 }
int64_t length
Definition: ptwXY.h:93

Here is the caller graph for this function:

nfu_status ptwXY_mergeClosePoints ( ptwXYPoints ptwXY,
double  epsilon 
)

Definition at line 141 of file ptwXY_convenient.cc.

141  {
142 
143  int64_t i, i1, j, k, n = ptwXY->length;
144  double x, y;
145  ptwXYPoint *p1, *p2;
146 
147  if( n < 2 ) return( ptwXY->status );
148  if( epsilon < 4 * DBL_EPSILON ) epsilon = 4 * DBL_EPSILON;
149  if( ptwXY_simpleCoalescePoints( ptwXY ) != nfu_Okay ) return( ptwXY->status );
150 
151  p2 = ptwXY->points;
152  x = p2->x;
153  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. */
154  if( ( p2->x - x ) > 0.5 * epsilon * ( std::fabs( x ) + std::fabs( p2->x ) ) ) break;
155  }
156  if( i1 != 1 ) {
157  for( i = i1, p1 = &(ptwXY->points[1]); i < n; i++, p1++, p2++ ) *p1 = *p2;
158  n = ptwXY->length = ptwXY->length - i1 + 1;
159  }
160 
161  p1 = &(ptwXY->points[n-1]);
162  x = p1->x;
163  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. */
164  if( x - p1->x > 0.5 * epsilon * ( std::fabs( x ) + std::fabs( p1->x ) ) ) break;
165  }
166  if( i1 != ( n - 2 ) ) {
167  ptwXY->points[i1 + 1] = ptwXY->points[n - 1];
168  n = ptwXY->length = i1 + 2;
169  }
170 
171  for( i = 1; i < n - 1; i++ ) {
172  p1 = &(ptwXY->points[i]);
173  x = p1->x;
174  y = p1->y;
175  for( j = i + 1, p2 = &(ptwXY->points[i+1]); j < n - 1; j++, p2++ ) {
176  if( ( p2->x - p1->x ) > 0.5 * epsilon * ( std::fabs( p2->x ) + std::fabs( p1->x ) ) ) break;
177  x += p2->x;
178  y += p2->y;
179  }
180  if( ( k = ( j - i ) ) > 1 ) {
181  p1->x = x / k;
182  p1->y = y / k;
183  for( p1 = &(ptwXY->points[i+1]); j < n; j++, p1++, p2++ ) *p1 = *p2;
184  n -= ( k - 1 );
185  }
186  }
187  ptwXY->length = n;
188 
189  return( ptwXY->status );
190 }
ptwXYPoint * points
Definition: ptwXY.h:99
tuple x
Definition: test.py:50
int64_t length
Definition: ptwXY.h:93
#define DBL_EPSILON
Definition: templates.hh:87
const G4int n
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
nfu_status status
Definition: ptwXY.h:85
double epsilon(double density, double temperature)

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 966 of file ptwXY_core.cc.

966  {
967 
968  return( ptwXY_mergeFrom( ptwXY, 1, length, xs, ys ) );
969 }
static nfu_status ptwXY_mergeFrom(ptwXYPoints *ptwXY, int incY, int length, double *xs, double *ys)
Definition: ptwXY_core.cc:991

Here is the call graph for this function:

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

Definition at line 973 of file ptwXY_core.cc.

973  {
974 
975  int i;
976  double *xs, *p1, *p2;
977  nfu_status status;
978 
979  if( length < 0 ) return( nfu_badInput );
980  if( length == 0 ) return( nfu_Okay );
981  if( ( xs = (double *) nfu_malloc( length * sizeof( double ) ) ) == NULL ) return( nfu_mallocError );
982  for( i = 0, p1 = xs, p2 = xys; i < length; i++, p1++, p2 += 2 ) *p1 = *p2;
983  status = ptwXY_mergeFrom( ptwXY, 2, length, xs, xys );
984  nfu_free( xs );
985 
986  return( status );
987 }
enum nfu_status_e nfu_status
void * nfu_free(void *p)
static nfu_status ptwXY_mergeFrom(ptwXYPoints *ptwXY, int incY, int length, double *xs, double *ys)
Definition: ptwXY_core.cc:991
void * nfu_malloc(size_t size)

Here is the call graph for this function:

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

Definition at line 76 of file ptwXY_binaryOperators.cc.

76  {
77 
78  int64_t i, nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
79  ptwXYPoint *p;
80  ptwXYOverflowPoint *o, *overflowHeader = &(ptwXY->overflowHeader);
81 
82  if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
83  if( m == 0 ) return( ptwXY->status = nfu_divByZero );
84 
85  for( i = 0, p = ptwXY->points; i < nonOverflowLength; i++, p++ ) p->y = ptwXY_mod2( p->y, m, pythonMod );
86  for( o = overflowHeader->next; o != overflowHeader; o = o->next ) o->point.y = ptwXY_mod2( o->point.y, m, pythonMod );
87  return( ptwXY->status );
88 }
ptwXYPoint point
Definition: ptwXY.h:80
ptwXYPoint * points
Definition: ptwXY.h:99
const char * p
Definition: xmltok.h:285
static double ptwXY_mod2(double v, double m, int pythonMod)
int64_t ptwXY_getNonOverflowLength(ptwXYPoints const *ptwXY)
Definition: ptwXY_core.cc:590
static constexpr double m
Definition: G4SIunits.hh:129
double y
Definition: ptwXY.h:62
struct ptwXYOverflowPoint_s * next
Definition: ptwXY.h:78
nfu_status status
Definition: ptwXY.h:85
ptwXYOverflowPoint overflowHeader
Definition: ptwXY.h:98

Here is the call graph for this function:

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

Definition at line 187 of file ptwXY_binaryOperators.cc.

187  {
188 
189  int64_t i, length;
190  ptwXYPoints *n = NULL;
191  int found;
192  double x1, y1, x2, y2, u1, u2, v1, v2, xz1 = 0, xz2 = 0, x;
193 
194  *status = nfu_otherInterpolation;
195  if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
196  if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( NULL );
197  if( ( n = ptwXY_mul_ptwXY( ptwXY1, ptwXY2, status ) ) == NULL ) return( n );
198  if( ptwXY1->interpolation == ptwXY_interpolationFlat ) return( n );
199  if( ptwXY2->interpolation == ptwXY_interpolationFlat ) return( n );
200  length = n->length - 1;
201  if( length > 0 ) {
202  x2 = n->points[length].x;
203  for( i = length - 1; i >= 0; i-- ) { /* Find and add y zeros not currently in n's. */
204  x1 = n->points[i].x;
205  if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x1, &u1 ) ) != nfu_Okay ) goto Err;
206  if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x2, &u2 ) ) != nfu_Okay ) goto Err;
207  if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, x1, &v1 ) ) != nfu_Okay ) goto Err;
208  if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, x2, &v2 ) ) != nfu_Okay ) goto Err;
209  found = 0;
210  if( u1 * u2 < 0 ) {
211  xz1 = ( u1 * x2 - u2 * x1 ) / ( u1 - u2 );
212  if( ( *status = ptwXY_setValueAtX( n, xz1, 0. ) ) != nfu_Okay ) goto Err;
213  found = 1;
214  }
215  if( v1 * v2 < 0 ) {
216  xz2 = ( v1 * x2 - v2 * x1 ) / ( v1 - v2 );
217  if( ( *status = ptwXY_setValueAtX( n, xz2, 0. ) ) != nfu_Okay ) goto Err;
218  found += 1;
219  }
220  if( found > 1 ) {
221  x = 0.5 * ( xz1 + xz2 );
222  if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x, &u1 ) ) != nfu_Okay ) goto Err;
223  if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, x, &v1 ) ) != nfu_Okay ) goto Err;
224  if( ( *status = ptwXY_setValueAtX( n, x, u1 * v1 ) ) != nfu_Okay ) goto Err;
225  }
226  x2 = x1;
227  }
228 
229  if( ( *status = ptwXY_simpleCoalescePoints( n ) ) != nfu_Okay ) goto Err;
230  length = n->length;
231  x2 = n->points[n->length-1].x;
232  y2 = n->points[n->length-1].y;
233  for( i = n->length - 2; i >= 0; i-- ) { /* Make interpolation fit accuracy. Work backwards so new points will not mess up loop. */
234  x1 = n->points[i].x;
235  y1 = n->points[i].y;
236  if( ( *status = ptwXY_mul2_s_ptwXY( n, ptwXY1, ptwXY2, x1, y1, x2, y2, 0 ) ) != nfu_Okay ) goto Err;
237  x2 = x1;
238  y2 = y1;
239  }
240  ptwXY_update_biSectionMax( n, (double) length );
241  }
242  return( n );
243 
244 Err:
245  if( n ) ptwXY_free( n );
246  return( NULL );
247 }
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
ptwXYPoint * points
Definition: ptwXY.h:99
tuple x
Definition: test.py:50
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
Definition: ptwXY_core.cc:876
int64_t length
Definition: ptwXY.h:93
static nfu_status ptwXY_getValueAtX_ignore_XOutsideDomainError(ptwXYPoints *ptwXY1, double x, double *y)
const G4int n
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
void ptwXY_update_biSectionMax(ptwXYPoints *ptwXY1, double oldLength)
Definition: ptwXY_misc.cc:31
ptwXYPoints * ptwXY_mul_ptwXY(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status)
static nfu_status ptwXY_mul2_s_ptwXY(ptwXYPoints *n, ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, double x1, double y1, double x2, double y2, int level)

Here is the call graph for this function:

nfu_status ptwXY_mul_double ( ptwXYPoints ptwXY,
double  value 
)

Definition at line 43 of file ptwXY_binaryOperators.cc.

43 { return( ptwXY_slopeOffset( ptwXY, value, 0. ) ); }
const XML_Char int const XML_Char * value
Definition: expat.h:331
nfu_status ptwXY_slopeOffset(ptwXYPoints *ptwXY, double slope, double offset)

Here is the call graph for this function:

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

Definition at line 171 of file ptwXY_binaryOperators.cc.

171  {
172 
173  ptwXYPoints *mul;
174 
175  if( ptwXY1->length == 0 ) {
176  mul = ptwXY_clone( ptwXY1, status ); }
177  else if( ptwXY2->length == 0 ) {
178  mul = ptwXY_clone( ptwXY2, status ); }
179  else {
180  mul = ptwXY_binary_ptwXY( ptwXY1, ptwXY2, 0., 0., 1., status );
181  }
182  return( mul );
183 }
ptwXYPoints * ptwXY_clone(ptwXYPoints *ptwXY, nfu_status *status)
Definition: ptwXY_core.cc:208
int64_t length
Definition: ptwXY.h:93
ptwXYPoints * ptwXY_binary_ptwXY(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, double v1, double v2, double v1v2, nfu_status *status)

Here is the call graph for this function:

Here is the caller graph for this function:

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.

369  {
370 
371  nfu_status status;
372  int64_t n1 = ptwXY1->length, n2 = ptwXY2->length;
373  ptwXYPoint *xy1, *xy2;
374 
375  switch( status = ptwXY_areDomainsMutual( ptwXY1, ptwXY2 ) ) {
376  case nfu_Okay :
377  case nfu_empty :
378  return( nfu_Okay );
379  case nfu_domainsNotMutual :
380  break;
381  default :
382  return( status );
383  }
384 
389 
390  xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, 0 );
391  xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, 0 );
392  if( xy1->x < xy2->x ) {
393  lowerEps1 = 0.;
394  if( xy2->y == 0. ) lowerEps2 = 0.; }
395  else if( xy1->x > xy2->x ) {
396  lowerEps2 = 0.;
397  if( xy1->y == 0. ) lowerEps1 = 0.; }
398  else {
399  lowerEps1 = lowerEps2 = 0.;
400  }
401 
402  xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, n1 - 1 );
403  xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, n2 - 1 );
404  if( xy1->x < xy2->x ) {
405  upperEps2 = 0.;
406  if( xy1->y == 0. ) upperEps1 = 0.; }
407  else if( xy1->x > xy2->x ) {
408  upperEps1 = 0.;
409  if( xy2->y == 0. ) upperEps2 = 0.; }
410  else {
411  upperEps1 = upperEps2 = 0.;
412  }
413 
414  if( ( lowerEps1 != 0. ) || ( upperEps1 != 0. ) )
415  if( ( status = ptwXY_dullEdges( ptwXY1, lowerEps1, upperEps1, positiveXOnly1 ) ) != nfu_Okay ) return( status );
416  if( ( lowerEps2 != 0. ) || ( upperEps2 != 0. ) )
417  if( ( status = ptwXY_dullEdges( ptwXY2, lowerEps2, upperEps2, positiveXOnly2 ) ) != nfu_Okay ) return( status );
418 
419  return( status );
420 }
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
ptwXYPoint * ptwXY_getPointAtIndex_Unsafely(ptwXYPoints *ptwXY, int64_t index)
Definition: ptwXY_core.cc:684
int64_t length
Definition: ptwXY.h:93
enum nfu_status_e nfu_status
nfu_status ptwXY_dullEdges(ptwXYPoints *ptwXY, double lowerEps, double upperEps, int positiveXOnly)
nfu_status ptwXY_areDomainsMutual(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2)
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62

Here is the call graph for this function:

nfu_status ptwXY_neg ( ptwXYPoints ptwXY)

Definition at line 34 of file ptwXY_unitaryOperators.cc.

34  {
35 
36  int64_t i, nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
37  ptwXYPoint *p;
38  ptwXYOverflowPoint *o, *overflowHeader = &(ptwXY->overflowHeader);
39 
40  if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
41 
42  for( i = 0, p = ptwXY->points; i < nonOverflowLength; i++, p++ ) p->y = -p->y;
43  for( o = overflowHeader->next; o != overflowHeader; o = o->next ) o->point.y = -o->point.y;
44  return( ptwXY->status );
45 }
ptwXYPoint point
Definition: ptwXY.h:80
ptwXYPoint * points
Definition: ptwXY.h:99
const char * p
Definition: xmltok.h:285
int64_t ptwXY_getNonOverflowLength(ptwXYPoints const *ptwXY)
Definition: ptwXY_core.cc:590
double y
Definition: ptwXY.h:62
struct ptwXYOverflowPoint_s * next
Definition: ptwXY.h:78
nfu_status status
Definition: ptwXY.h:85
ptwXYOverflowPoint overflowHeader
Definition: ptwXY.h:98

Here is the call graph for this function:

Here is the caller graph for this function:

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.

30  {
31 
32  ptwXYPoints *ptwXY = (ptwXYPoints *) nfu_calloc( sizeof( ptwXYPoints ), 1 );
33 
34  *status = nfu_mallocError;
35  if( ptwXY == NULL ) return( NULL );
36  ptwXY_setup( ptwXY, interpolation, interpolationOtherInfo, biSectionMax, accuracy, primarySize,
37  secondarySize, userFlag );
38  if( ( *status = ptwXY->status ) != nfu_Okay ) {
39  ptwXY = (ptwXYPoints *) nfu_free( ptwXY );
40  }
41  return( ptwXY );
42 }
void * nfu_calloc(size_t size, size_t n)
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: ptwXY_core.cc:46
void * nfu_free(void *p)
nfu_status status
Definition: ptwXY.h:85

Here is the call graph for this function:

Here is the caller graph for this function:

nfu_status ptwXY_normalize ( ptwXYPoints ptwXY1)

Definition at line 190 of file ptwXY_integration.cc.

190  {
191 /*
192 * This function assumes ptwXY_integrateDomain checks status and coalesces the points.
193 */
194 
195  int64_t i;
196  nfu_status status;
197  double sum = ptwXY_integrateDomain( ptwXY, &status );
198 
199  if( status != nfu_Okay ) return( status );
200  if( sum == 0. ) {
201  status = nfu_badNorm; }
202  else {
203  for( i = 0; i < ptwXY->length; i++ ) ptwXY->points[i].y /= sum;
204  }
205  return( status );
206 }
double ptwXY_integrateDomain(ptwXYPoints *ptwXY, nfu_status *status)
enum nfu_status_e nfu_status

Here is the call graph for this function:

Here is the caller graph for this function:

nfu_status ptwXY_pow ( ptwXYPoints ptwXY,
double  p 
)

Definition at line 24 of file ptwXY_functions.cc.

24  {
25 
26  return( ptwXY_applyFunction( ptwXY, ptwXY_pow_callback, (void *) &v, 0 ) );
27 }
nfu_status ptwXY_applyFunction(ptwXYPoints *ptwXY1, ptwXY_applyFunction_callback func, void *argList, int checkForRoots)
Definition: ptwXY_misc.cc:146
static nfu_status ptwXY_pow_callback(ptwXYPoint *point, void *argList)
tuple v
Definition: test.py:18

Here is the call graph for this function:

nfu_status ptwXY_reallocateOverflowPoints ( ptwXYPoints ptwXY,
int64_t  size 
)

Definition at line 439 of file ptwXY_core.cc.

439  {
440 /*
441 * This is for allocating/reallocating the secondary data memory.
442 */
443  nfu_status status = nfu_Okay;
444 
445  if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
446 
447  if( size < ptwXY_minimumOverflowSize ) size = ptwXY_minimumOverflowSize; /* ptwXY_minimumOverflowSize must be > 0. */
448  if( size < ptwXY->overflowLength ) status = ptwXY_coalescePoints( ptwXY, ptwXY->length + ptwXY->overflowAllocatedSize, NULL, 0 );
449  if( status == nfu_Okay ) {
450  if( size != ptwXY->overflowAllocatedSize ) {
451  ptwXY->overflowPoints = (ptwXYOverflowPoint *) nfu_realloc( (size_t) size * sizeof( ptwXYOverflowPoint ), ptwXY->overflowPoints );
452  if( ptwXY->overflowPoints == NULL ) {
453  ptwXY->length = 0;
454  ptwXY->overflowLength = 0;
455  ptwXY->mallocFailedSize = size;
456  size = 0;
457  ptwXY->status = nfu_mallocError;
458  }
459  }
460  ptwXY->overflowAllocatedSize = size; }
461  else {
462  ptwXY->status = status;
463  }
464  return( ptwXY->status );
465 }
#define ptwXY_minimumOverflowSize
Definition: ptwXY.h:21
int64_t mallocFailedSize
Definition: ptwXY.h:97
nfu_status ptwXY_coalescePoints(ptwXYPoints *ptwXY, int64_t size, ptwXYPoint *newPoint, int forceSmallerResize)
Definition: ptwXY_core.cc:469
int64_t length
Definition: ptwXY.h:93
enum nfu_status_e nfu_status
int64_t overflowAllocatedSize
Definition: ptwXY.h:96
int64_t overflowLength
Definition: ptwXY.h:95
nfu_status status
Definition: ptwXY.h:85
void * nfu_realloc(size_t size, void *old)
struct ptwXYOverflowPoint_s ptwXYOverflowPoint
ptwXYOverflowPoint * overflowPoints
Definition: ptwXY.h:100

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 410 of file ptwXY_core.cc.

410  {
411 /*
412 * This is for allocating/reallocating the primary data memory.
413 */
414  if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
415 
416  if( size < ptwXY_minimumSize ) size = ptwXY_minimumSize; /* ptwXY_minimumSize must be > 0. */
417  if( size < ptwXY->length ) size = ptwXY->length;
418  if( size != ptwXY->allocatedSize ) {
419  if( size > ptwXY->allocatedSize ) { /* Increase size of allocated points. */
420  ptwXY->points = (ptwXYPoint *) nfu_realloc( (size_t) size * sizeof( ptwXYPoint ), ptwXY->points ); }
421  else if( ( ptwXY->allocatedSize > 2 * size ) || forceSmallerResize ) { /* Decrease size, if at least 1/2 size reduction or if forced to. */
422  ptwXY->points = (ptwXYPoint *) nfu_realloc( (size_t) size * sizeof( ptwXYPoint ), ptwXY->points ); }
423  else {
424  size = ptwXY->allocatedSize; /* Size is < ptwXY->allocatedSize, but realloc not called. */
425  }
426  if( ptwXY->points == NULL ) {
427  ptwXY->length = 0;
428  ptwXY->mallocFailedSize = size;
429  size = 0;
430  ptwXY->status = nfu_mallocError;
431  }
432  ptwXY->allocatedSize = size;
433  }
434  return( ptwXY->status );
435 }
int64_t mallocFailedSize
Definition: ptwXY.h:97
ptwXYPoint * points
Definition: ptwXY.h:99
int64_t length
Definition: ptwXY.h:93
int64_t allocatedSize
Definition: ptwXY.h:94
nfu_status status
Definition: ptwXY.h:85
void * nfu_realloc(size_t size, void *old)
#define ptwXY_minimumSize
Definition: ptwXY.h:20
struct ptwXYPoint_s ptwXYPoint

Here is the call graph for this function:

Here is the caller graph for this function:

nfu_status ptwXY_release ( ptwXYPoints ptwXY)

Definition at line 549 of file ptwXY_core.cc.

549  {
550 /*
551 * Note, this routine does not free ptwXY (i.e., it does not undo all of ptwXY_new).
552 */
553 
554  if( ptwXY->interpolation == ptwXY_interpolationOther ) {
555  if( ptwXY->interpolationOtherInfo.interpolationString != NULL )
557  }
559  ptwXY->interpolationOtherInfo.getValueFunc = NULL;
560  ptwXY->interpolationOtherInfo.argList = NULL;
561  ptwXY->length = 0;
562  ptwXY->allocatedSize = 0;
563  ptwXY->points = (ptwXYPoint *) nfu_free( ptwXY->points );
564 
565  ptwXY->overflowLength = 0;
566  ptwXY->overflowAllocatedSize = 0;
568 
569  return( nfu_Okay );
570 }
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
ptwXYPoint * points
Definition: ptwXY.h:99
ptwXY_getValue_callback getValueFunc
Definition: ptwXY.h:71
int64_t length
Definition: ptwXY.h:93
void * nfu_free(void *p)
ptwXY_interpolationOtherInfo interpolationOtherInfo
Definition: ptwXY.h:88
int64_t allocatedSize
Definition: ptwXY.h:94
int64_t overflowAllocatedSize
Definition: ptwXY.h:96
int64_t overflowLength
Definition: ptwXY.h:95
char const * interpolationString
Definition: ptwXY.h:70
ptwXYOverflowPoint * overflowPoints
Definition: ptwXY.h:100

Here is the call graph for this function:

Here is the caller graph for this function:

ptwXPoints* ptwXY_runningIntegral ( ptwXYPoints ptwXY,
nfu_status status 
)

Definition at line 632 of file ptwXY_integration.cc.

632  {
633 
634  int i;
635  ptwXPoints *runningIntegral = NULL;
636  double integral = 0., sum;
637 
638  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( NULL );
639  if( ( runningIntegral = ptwX_new( ptwXY->length, status ) ) == NULL ) goto err;
640 
641  if( ( *status = ptwX_setPointAtIndex( runningIntegral, 0, 0. ) ) != nfu_Okay ) goto err;
642  for( i = 1; i < ptwXY->length; i++ ) {
643  if( ( *status = ptwXY_f_integrate( ptwXY->interpolation, ptwXY->points[i-1].x, ptwXY->points[i-1].y,
644  ptwXY->points[i].x, ptwXY->points[i].y, &sum ) ) != nfu_Okay ) goto err;
645  integral += sum;
646  if( ( *status = ptwX_setPointAtIndex( runningIntegral, i, integral ) ) != nfu_Okay ) goto err;
647  }
648  return( runningIntegral );
649 
650 err:
651  if( runningIntegral != NULL ) ptwX_free( runningIntegral );
652  return( NULL );
653 }
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
ptwXYPoint * points
Definition: ptwXY.h:99
int64_t length
Definition: ptwXY.h:93
ptwXPoints * ptwX_new(int64_t size, nfu_status *status)
Definition: ptwX_core.cc:24
nfu_status ptwXY_f_integrate(ptwXY_interpolation interpolation, double x1, double y1, double x2, double y2, double *value)
ptwXPoints * ptwX_free(ptwXPoints *ptwX)
Definition: ptwX_core.cc:158
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
nfu_status ptwX_setPointAtIndex(ptwXPoints *ptwX, int64_t index, double x)
Definition: ptwX_core.cc:222

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 481 of file ptwXY_methods.cc.

481  {
482 
483  int64_t i1, length = ptwXY->length;
484  ptwXYPoint *p1;
485  nfu_status status;
486 
487  if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
488  if( xScale == 0 ) return( nfu_XNotAscending );
489 
490  if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status );
491 
492  for( i1 = 0, p1 = ptwXY->points; i1 < length; i1++, p1++ ) {
493  p1->x = xScale * p1->x + xOffset;
494  p1->y = yScale * p1->y + yOffset;
495  }
496 
497  if( xScale < 0 ) {
498  int64_t length_2 = length / 2;
499  ptwXYPoint tmp, *p2;
500 
501  for( i1 = 0, p1 = ptwXY->points, p2 = &(ptwXY->points[length-1]); i1 < length_2; i1++ ) {
502  tmp = *p1;
503  *p1 = *p2;
504  *p2 = tmp;
505  }
506  }
507 
508  return( ptwXY->status );
509 }
ptwXYPoint * points
Definition: ptwXY.h:99
int64_t length
Definition: ptwXY.h:93
enum nfu_status_e nfu_status
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
nfu_status status
Definition: ptwXY.h:85

Here is the call graph for this function:

Here is the caller graph for this function:

double ptwXY_setAccuracy ( ptwXYPoints ptwXY,
double  accuracy 
)

Definition at line 379 of file ptwXY_core.cc.

379  {
380 
381  if( accuracy < ptwXY_minAccuracy ) accuracy = ptwXY_minAccuracy;
382  if( accuracy < ptwXY->accuracy ) accuracy = ptwXY->accuracy;
383  if( accuracy > 1 ) accuracy = 1.;
384  ptwXY->accuracy = accuracy;
385  return( ptwXY->accuracy );
386 }
double accuracy
Definition: ptwXY.h:91
#define ptwXY_minAccuracy
Definition: ptwXY.h:23

Here is the caller graph for this function:

double ptwXY_setBiSectionMax ( ptwXYPoints ptwXY,
double  biSectionMax 
)

Definition at line 397 of file ptwXY_core.cc.

397  {
398 
399  if( biSectionMax < 0 ) {
400  biSectionMax = 0; }
401  else if( biSectionMax > ptwXY_maxBiSectionMax ) {
402  biSectionMax = ptwXY_maxBiSectionMax;
403  }
404  ptwXY->biSectionMax = biSectionMax;
405  return( ptwXY->biSectionMax );
406 }
double biSectionMax
Definition: ptwXY.h:90
#define ptwXY_maxBiSectionMax
Definition: ptwXY.h:22

Here is the caller graph for this function:

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.

47  {
48 
49  ptwXY->status = nfu_Okay;
50  ptwXY->typeX = ptwXY_sigma_none;
51  ptwXY->typeY = ptwXY_sigma_none;
52  ptwXY->interpolation = interpolation;
55  ptwXY->interpolationOtherInfo.argList = NULL;
56  switch( interpolation ) {
67  case ptwXY_interpolationOther : /* For ptwXY_interpolationOther, interpolationOtherInfo and interpolationString must be defined. */
68  if( interpolationOtherInfo == NULL ) {
69  ptwXY->status = nfu_otherInterpolation; }
70  else {
71  if( interpolationOtherInfo->interpolationString == NULL ) {
72  ptwXY->status = nfu_otherInterpolation; }
73  else {
74  if( ( ptwXY->interpolationOtherInfo.interpolationString = strdup( interpolationOtherInfo->interpolationString ) ) == NULL ) {
75  ptwXY->status = nfu_mallocError;
76  }
77  }
78  ptwXY->interpolationOtherInfo.getValueFunc = interpolationOtherInfo->getValueFunc;
79  ptwXY->interpolationOtherInfo.argList = interpolationOtherInfo->argList;
80  }
81  }
82  ptwXY->userFlag = 0;
83  ptwXY_setUserFlag( ptwXY, userFlag );
85  ptwXY_setBiSectionMax( ptwXY, biSectionMax );
86  ptwXY->accuracy = ptwXY_minAccuracy;
87  ptwXY_setAccuracy( ptwXY, accuracy );
88 
89  ptwXY->length = 0;
90  ptwXY->allocatedSize = 0;
91  ptwXY->overflowLength = 0;
92  ptwXY->overflowAllocatedSize = 0;
93  ptwXY->mallocFailedSize = 0;
94 
96 
97  ptwXY->points = NULL;
98  ptwXY->overflowPoints = NULL;
99 
100  ptwXY_reallocatePoints( ptwXY, primarySize, 0 );
101  ptwXY_reallocateOverflowPoints( ptwXY, secondarySize );
102  if( ptwXY->status != nfu_Okay ) ptwXY_release( ptwXY );
103  return( ptwXY->status );
104 }
static char const logLinInterpolationString[]
Definition: ptwXY_core.cc:19
int64_t mallocFailedSize
Definition: ptwXY.h:97
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
int userFlag
Definition: ptwXY.h:89
ptwXYPoint * points
Definition: ptwXY.h:99
double ptwXY_setAccuracy(ptwXYPoints *ptwXY, double accuracy)
Definition: ptwXY_core.cc:379
static char const logLogInterpolationString[]
Definition: ptwXY_core.cc:20
nfu_status ptwXY_reallocateOverflowPoints(ptwXYPoints *ptwXY, int64_t size)
Definition: ptwXY_core.cc:439
ptwXY_getValue_callback getValueFunc
Definition: ptwXY.h:71
static void ptwXY_initialOverflowPoint(ptwXYOverflowPoint *overflowPoint, ptwXYOverflowPoint *prior, ptwXYOverflowPoint *next)
Definition: ptwXY_core.cc:1290
double biSectionMax
Definition: ptwXY.h:90
int64_t length
Definition: ptwXY.h:93
void ptwXY_setUserFlag(ptwXYPoints *ptwXY, int userFlag)
Definition: ptwXY_core.cc:365
nfu_status ptwXY_release(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:549
double ptwXY_setBiSectionMax(ptwXYPoints *ptwXY, double biSectionMax)
Definition: ptwXY_core.cc:397
static char const linLinInterpolationString[]
Definition: ptwXY_core.cc:17
ptwXY_interpolationOtherInfo interpolationOtherInfo
Definition: ptwXY.h:88
int64_t allocatedSize
Definition: ptwXY.h:94
ptwXY_sigma typeY
Definition: ptwXY.h:86
int64_t overflowAllocatedSize
Definition: ptwXY.h:96
int64_t overflowLength
Definition: ptwXY.h:95
nfu_status status
Definition: ptwXY.h:85
#define ptwXY_maxBiSectionMax
Definition: ptwXY.h:22
char const * interpolationString
Definition: ptwXY.h:70
ptwXYOverflowPoint overflowHeader
Definition: ptwXY.h:98
double accuracy
Definition: ptwXY.h:91
static char const flatInterpolationString[]
Definition: ptwXY_core.cc:21
nfu_status ptwXY_reallocatePoints(ptwXYPoints *ptwXY, int64_t size, int forceSmallerResize)
Definition: ptwXY_core.cc:410
#define ptwXY_minAccuracy
Definition: ptwXY.h:23
ptwXY_sigma typeX
Definition: ptwXY.h:86
ptwXYOverflowPoint * overflowPoints
Definition: ptwXY.h:100
static char const linLogInterpolationString[]
Definition: ptwXY_core.cc:18

Here is the call graph for this function:

Here is the caller graph for this function:

void ptwXY_setUserFlag ( ptwXYPoints ptwXY,
int  userFlag 
)

Definition at line 365 of file ptwXY_core.cc.

365  {
366 
367  ptwXY->userFlag = userFlag;
368 }
int userFlag
Definition: ptwXY.h:89

Here is the caller graph for this function:

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

Definition at line 876 of file ptwXY_core.cc.

876  {
877 
878  return( ptwXY_setValueAtX_overrideIfClose( ptwXY, x, y, 0., 0 ) );
879 }
tuple x
Definition: test.py:50
nfu_status ptwXY_setValueAtX_overrideIfClose(ptwXYPoints *ptwXY, double x, double y, double eps, int override)
Definition: ptwXY_core.cc:883

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 883 of file ptwXY_core.cc.

883  {
884 
885  int closeIsEqual;
886  int64_t nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY ), i;
887  nfu_status status = nfu_Okay;
889  ptwXYPoint *point = NULL, newPoint = { x, y };
890  ptwXYOverflowPoint *overflowPoint, *p, *overflowHeader = &(ptwXY->overflowHeader);
891  ptwXYOverflowPoint lessThanEqualXPoint, greaterThanXPoint;
892  ptwXYPoint *closePoint = NULL;
893 
894  if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
895 
896  legx = ptwXY_getPointsAroundX_closeIsEqual( ptwXY, x, &lessThanEqualXPoint, &greaterThanXPoint, eps, &closeIsEqual, &closePoint );
897  switch( legx ) {
901  if( closeIsEqual ) {
902  if( !override ) return( status );
903  point = closePoint;
905  x = point->x; }
906  else {
907  if( ( legx == ptwXY_lessEqualGreaterX_greater ) && ( nonOverflowLength < ptwXY->allocatedSize ) ) {
908  point = &(ptwXY->points[nonOverflowLength]); }
909  else {
910  if( ptwXY->overflowLength == ptwXY->overflowAllocatedSize )
911  return( ptwXY_coalescePoints( ptwXY, ptwXY->length + ptwXY->overflowAllocatedSize, &newPoint, 0 ) );
912  overflowPoint = &(ptwXY->overflowPoints[ptwXY->overflowLength]);
913  if( legx == ptwXY_lessEqualGreaterX_lessThan ) {
914  overflowPoint->prior = greaterThanXPoint.prior;
915  overflowPoint->index = 0; }
916  else { /* Between or greater and must go in overflow area. */
917  if( legx == ptwXY_lessEqualGreaterX_greater ) {
918  overflowPoint->prior = overflowHeader->prior;
919  overflowPoint->index = ptwXY->length; }
920  else {
921  overflowPoint->prior = lessThanEqualXPoint.prior;
922  if( lessThanEqualXPoint.next != NULL ) {
923  if( lessThanEqualXPoint.point.x < x )
924  overflowPoint->prior = lessThanEqualXPoint.prior->next;
925  i = 1; }
926  else {
927  for( p = overflowHeader->next, i = 1; p != overflowHeader; p = p->next, i++ )
928  if( p->point.x > x ) break;
929  }
930  overflowPoint->index = lessThanEqualXPoint.index + i;
931  }
932  }
933  overflowPoint->next = overflowPoint->prior->next;
934  overflowPoint->prior->next = overflowPoint;
935  overflowPoint->next->prior = overflowPoint;
936  point = &(overflowPoint->point);
937  for( overflowPoint = overflowPoint->next; overflowPoint != overflowHeader; overflowPoint = overflowPoint->next ) {
938  overflowPoint->index++;
939  }
940  ptwXY->overflowLength++;
941  }
942  }
943  break;
945  point = ptwXY->points; /* ptwXY_minimumSize must be > 0 so there is always space here. */
946  break;
948  if( closeIsEqual && !override ) return( status );
949  if( lessThanEqualXPoint.next == NULL ) {
950  point = &(ptwXY->points[lessThanEqualXPoint.index]); }
951  else {
952  point = &(lessThanEqualXPoint.prior->next->point);
953  }
954  break;
955  }
956  if( status == nfu_Okay ) {
957  point->x = x;
958  point->y = y;
959  if( legx != ptwXY_lessEqualGreaterX_equal ) ptwXY->length++;
960  }
961  return( status );
962 }
nfu_status ptwXY_coalescePoints(ptwXYPoints *ptwXY, int64_t size, ptwXYPoint *newPoint, int forceSmallerResize)
Definition: ptwXY_core.cc:469
ptwXYPoint point
Definition: ptwXY.h:80
ptwXYPoint * points
Definition: ptwXY.h:99
const char * p
Definition: xmltok.h:285
static const G4double eps
enum ptwXY_lessEqualGreaterX_e ptwXY_lessEqualGreaterX
tuple x
Definition: test.py:50
struct ptwXYOverflowPoint_s * prior
Definition: ptwXY.h:77
int64_t length
Definition: ptwXY.h:93
int64_t ptwXY_getNonOverflowLength(ptwXYPoints const *ptwXY)
Definition: ptwXY_core.cc:590
ptwXY_lessEqualGreaterX ptwXY_getPointsAroundX_closeIsEqual(ptwXYPoints *ptwXY, double x, ptwXYOverflowPoint *lessThanEqualXPoint, ptwXYOverflowPoint *greaterThanXPoint, double eps, int *closeIsEqual, ptwXYPoint **closePoint)
Definition: ptwXY_core.cc:720
enum nfu_status_e nfu_status
double y
Definition: ptwXY.h:62
struct ptwXYOverflowPoint_s * next
Definition: ptwXY.h:78
double x
Definition: ptwXY.h:62
int64_t overflowAllocatedSize
Definition: ptwXY.h:96
int64_t overflowLength
Definition: ptwXY.h:95
nfu_status status
Definition: ptwXY.h:85
ptwXYOverflowPoint overflowHeader
Definition: ptwXY.h:98
int64_t index
Definition: ptwXY.h:79
ptwXYOverflowPoint * overflowPoints
Definition: ptwXY.h:100

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 597 of file ptwXY_core.cc.

597  {
598 
599  nfu_status status = nfu_Okay;
600  int64_t i;
601  ptwXYPoint *p;
602  double const *d = xy;
603  double xOld = 0.;
604 
605  if( length > ptwXY->allocatedSize ) {
606  status = ptwXY_reallocatePoints( ptwXY, length, 0 );
607  if( status != nfu_Okay ) return( status );
608  }
609  for( i = 0, p = ptwXY->points; i < length; i++, p++ ) {
610  if( i != 0 ) {
611  if( *d <= xOld ) {
612  status = nfu_XNotAscending;
613  ptwXY->status = nfu_XNotAscending;
614  length = 0;
615  break;
616  }
617  }
618  xOld = *d;
619  p->x = *(d++);
620  p->y = *(d++);
621  }
622  ptwXY->overflowHeader.next = &(ptwXY->overflowHeader);
623  ptwXY->overflowHeader.prior = &(ptwXY->overflowHeader);
624  ptwXY->overflowLength = 0;
625  ptwXY->length = length;
626  return( ptwXY->status = status );
627 }
ptwXYPoint * points
Definition: ptwXY.h:99
const char * p
Definition: xmltok.h:285
struct ptwXYOverflowPoint_s * prior
Definition: ptwXY.h:77
int64_t length
Definition: ptwXY.h:93
enum nfu_status_e nfu_status
int64_t allocatedSize
Definition: ptwXY.h:94
double y
Definition: ptwXY.h:62
struct ptwXYOverflowPoint_s * next
Definition: ptwXY.h:78
double x
Definition: ptwXY.h:62
int64_t overflowLength
Definition: ptwXY.h:95
nfu_status status
Definition: ptwXY.h:85
ptwXYOverflowPoint overflowHeader
Definition: ptwXY.h:98
nfu_status ptwXY_reallocatePoints(ptwXYPoints *ptwXY, int64_t size, int forceSmallerResize)
Definition: ptwXY_core.cc:410

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 631 of file ptwXY_core.cc.

631  {
632 
633  nfu_status status;
634  int64_t i;
635  ptwXYPoint *p;
636  double xOld = 0.;
637 
638  if( ( status = ptwXY_clear( ptwXY ) ) != nfu_Okay ) return( status );
639  if( length > ptwXY->allocatedSize ) {
640  if( ( status = ptwXY_reallocatePoints( ptwXY, length, 0 ) ) != nfu_Okay ) return( status );
641  }
642  for( i = 0, p = ptwXY->points; i < length; i++, p++, x++, y++ ) {
643  if( i != 0 ) {
644  if( *x <= xOld ) {
645  status = ptwXY->status = nfu_XNotAscending;
646  length = 0;
647  break;
648  }
649  }
650  xOld = *x;
651  p->x = *x;
652  p->y = *y;
653  }
654  ptwXY->length = length;
655  return( status );
656 }
ptwXYPoint * points
Definition: ptwXY.h:99
const char * p
Definition: xmltok.h:285
tuple x
Definition: test.py:50
int64_t length
Definition: ptwXY.h:93
enum nfu_status_e nfu_status
int64_t allocatedSize
Definition: ptwXY.h:94
nfu_status ptwXY_clear(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:536
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
nfu_status status
Definition: ptwXY.h:85
nfu_status ptwXY_reallocatePoints(ptwXYPoints *ptwXY, int64_t size, int forceSmallerResize)
Definition: ptwXY_core.cc:410

Here is the call graph for this function:

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

Definition at line 1098 of file ptwXY_core.cc.

1098  {
1099 
1100  int64_t i, ip1;
1101  ptwXYOverflowPoint *overflowPoint, *pm1, *pp1;
1102 
1103  if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
1104 
1105  if( ( index < 0 ) || ( index >= ptwXY->length ) ) return( nfu_badIndex );
1106  for( overflowPoint = ptwXY->overflowHeader.next, i = 0; overflowPoint != &(ptwXY->overflowHeader); overflowPoint = overflowPoint->next, i++ ) {
1107  if( overflowPoint->index >= index ) break;
1108  }
1109  ip1 = i;
1110  pm1 = pp1 = overflowPoint;
1111  if( overflowPoint->index == index ) { /* Note, if overflowPoint is header, then its index = -1. */
1112  pp1 = overflowPoint->next;
1113  ip1++;
1114  }
1115  if( ( pp1 != &(ptwXY->overflowHeader) ) && ( pp1->index == ( index + 1 ) ) ) { /* This if and else check that x < element[index+1]'s x values. */
1116  if( pp1->point.x <= x ) return( nfu_badIndexForX ); }
1117  else {
1118  if( ( ( index + 1 ) < ptwXY->length ) && ( ptwXY->points[index + 1 - ip1].x <= x ) ) return( nfu_badIndexForX );
1119  }
1120  if( overflowPoint != &(ptwXY->overflowHeader) ) pm1 = overflowPoint->prior;
1121  if( ( pm1 != &(ptwXY->overflowHeader) ) && ( pm1->index == ( index - 1 ) ) ) { /* This if and else check that x > element[index-1]'s x values. */
1122  if( pm1->point.x >= x ) return( nfu_badIndexForX ); }
1123  else {
1124  if( ( ( index - 1 ) >= 0 ) && ( ptwXY->points[index - 1 - i].x >= x ) ) return( nfu_badIndexForX );
1125  }
1126  if( ( overflowPoint != &(ptwXY->overflowHeader) ) && ( overflowPoint->index == index ) ) {
1127  overflowPoint->point.x = x;
1128  overflowPoint->point.y = y; }
1129  else {
1130  index -= i;
1131  ptwXY->points[index].x = x;
1132  ptwXY->points[index].y = y;
1133  }
1134  return( nfu_Okay );
1135 }
ptwXYPoint point
Definition: ptwXY.h:80
ptwXYPoint * points
Definition: ptwXY.h:99
tuple x
Definition: test.py:50
struct ptwXYOverflowPoint_s * prior
Definition: ptwXY.h:77
int64_t length
Definition: ptwXY.h:93
double y
Definition: ptwXY.h:62
struct ptwXYOverflowPoint_s * next
Definition: ptwXY.h:78
double x
Definition: ptwXY.h:62
nfu_status status
Definition: ptwXY.h:85
ptwXYOverflowPoint overflowHeader
Definition: ptwXY.h:98
int64_t index
Definition: ptwXY.h:79
void ptwXY_showInteralStructure ( ptwXYPoints ptwXY,
FILE *  f,
int  printPointersAsNull 
)

Definition at line 253 of file ptwXY_misc.cc.

253  {
254 
255  int64_t i, n = ptwXY_getNonOverflowLength( ptwXY );
256  ptwXYPoint *point = ptwXY->points;
257  ptwXYOverflowPoint *overflowPoint;
258 
259  fprintf( f, "status = %d interpolation = %d length = %d allocatedSize = %d\n",
260  (int) ptwXY->status, (int) ptwXY->interpolation, (int) ptwXY->length, (int) ptwXY->allocatedSize );
261  fprintf( f, "userFlag = %d biSectionMax = %.8e accuracy = %.2e minFractional_dx = %.6e\n",
262  ptwXY->userFlag, ptwXY->biSectionMax, ptwXY->accuracy, ptwXY->minFractional_dx );
263  fprintf( f, "interpolationString = %s\n", ptwXY->interpolationOtherInfo.interpolationString );
264  fprintf( f, "getValueFunc is NULL = %d. argList is NULL = %d.\n",
265  ptwXY->interpolationOtherInfo.getValueFunc == NULL, ptwXY->interpolationOtherInfo.argList == NULL );
266  fprintf( f, " overflowLength = %d overflowAllocatedSize = %d mallocFailedSize = %d\n",
267  (int) ptwXY->overflowLength, (int) ptwXY->overflowAllocatedSize, (int) ptwXY->mallocFailedSize );
268  fprintf( f, " Points data, points = %20p\n", ( printPointersAsNull ? NULL : (void*)ptwXY->points ) );
269  for( i = 0; i < n; i++, point++ ) fprintf( f, " %14.7e %14.7e\n", point->x, point->y );
270  fprintf( f, " Overflow points data; %20p\n", ( printPointersAsNull ? NULL : (void*)&(ptwXY->overflowHeader) ) );
271  for( overflowPoint = ptwXY->overflowHeader.next; overflowPoint != &(ptwXY->overflowHeader); overflowPoint = overflowPoint->next ) {
272  fprintf( f, " %14.7e %14.7e %8d %20p %20p %20p\n", overflowPoint->point.x, overflowPoint->point.y, (int) overflowPoint->index,
273  (void*) ( printPointersAsNull ? NULL : overflowPoint ), (void*) ( printPointersAsNull ? NULL : overflowPoint->prior ),
274  (void*) ( printPointersAsNull ? NULL : overflowPoint->next ) );
275  }
276  fprintf( f, " Points in order\n" );
277  for( i = 0; i < ptwXY->length; i++ ) {
278  point = ptwXY_getPointAtIndex( ptwXY, i );
279  fprintf( f, " %14.7e %14.7e\n", point->x, point->y );
280  }
281 }
double minFractional_dx
Definition: ptwXY.h:92
ptwXYPoint * ptwXY_getPointAtIndex(ptwXYPoints *ptwXY, int64_t index)
Definition: ptwXY_core.cc:675
int64_t mallocFailedSize
Definition: ptwXY.h:97
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
ptwXYPoint point
Definition: ptwXY.h:80
int userFlag
Definition: ptwXY.h:89
ptwXYPoint * points
Definition: ptwXY.h:99
ptwXY_getValue_callback getValueFunc
Definition: ptwXY.h:71
double biSectionMax
Definition: ptwXY.h:90
int64_t length
Definition: ptwXY.h:93
int64_t ptwXY_getNonOverflowLength(ptwXYPoints const *ptwXY)
Definition: ptwXY_core.cc:590
const G4int n
ptwXY_interpolationOtherInfo interpolationOtherInfo
Definition: ptwXY.h:88
int64_t allocatedSize
Definition: ptwXY.h:94
double y
Definition: ptwXY.h:62
struct ptwXYOverflowPoint_s * next
Definition: ptwXY.h:78
double x
Definition: ptwXY.h:62
int64_t overflowAllocatedSize
Definition: ptwXY.h:96
int64_t overflowLength
Definition: ptwXY.h:95
nfu_status status
Definition: ptwXY.h:85
char const * interpolationString
Definition: ptwXY.h:70
ptwXYOverflowPoint overflowHeader
Definition: ptwXY.h:98
int64_t index
Definition: ptwXY.h:79
double accuracy
Definition: ptwXY.h:91

Here is the call graph for this function:

nfu_status ptwXY_simpleCoalescePoints ( ptwXYPoints ptwXY)

Definition at line 529 of file ptwXY_core.cc.

529  {
530 
531  return( ptwXY_coalescePoints( ptwXY, ptwXY->length, NULL, 0 ) );
532 }
nfu_status ptwXY_coalescePoints(ptwXYPoints *ptwXY, int64_t size, ptwXYPoint *newPoint, int forceSmallerResize)
Definition: ptwXY_core.cc:469
int64_t length
Definition: ptwXY.h:93

Here is the call graph for this function:

Here is the caller graph for this function:

void ptwXY_simplePrint ( ptwXYPoints ptwXY,
char *  format 
)

Definition at line 298 of file ptwXY_misc.cc.

298  {
299 
300  ptwXY_simpleWrite( ptwXY, stdout, format );
301 }
void ptwXY_simpleWrite(ptwXYPoints *ptwXY, FILE *f, char *format)
Definition: ptwXY_misc.cc:285

Here is the call graph for this function:

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

Definition at line 285 of file ptwXY_misc.cc.

285  {
286 
287  int64_t i;
288  ptwXYPoint *point;
289 
290  for( i = 0; i < ptwXY->length; i++ ) {
291  point = ptwXY_getPointAtIndex( ptwXY, i );
292  fprintf( f, format, point->x, point->y );
293  }
294 }
ptwXYPoint * ptwXY_getPointAtIndex(ptwXYPoints *ptwXY, int64_t index)
Definition: ptwXY_core.cc:675
int64_t length
Definition: ptwXY.h:93
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62

Here is the call graph for this function:

Here is the caller graph for this function:

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.

248  {
249 
250  int64_t i, length;
251  ptwXYPoints *n;
252 
253  *status = nfu_badSelf;
254  if( ptwXY->status != nfu_Okay ) return( NULL );
255 
256  *status = nfu_badIndex;
257  if( index2 < index1 ) return( NULL );
258  if( index1 < 0 ) index1 = 0;
259  if( index2 > ptwXY->length ) index2 = ptwXY->length;
260 
261  length = index2 - index1;
262  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( NULL );
263  if( ( n = ptwXY_new( ptwXY->interpolation, &(ptwXY->interpolationOtherInfo), ptwXY->biSectionMax,
264  ptwXY->accuracy, length, secondarySize, status, ptwXY->userFlag ) ) == NULL ) return( NULL );
265 
266  *status = n->status = ptwXY->status;
267  for( i = index1; i < index2; i++ ) n->points[i - index1] = ptwXY->points[i];
268  n->length = length;
269  return( n );
270 }
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
int userFlag
Definition: ptwXY.h:89
ptwXYPoint * points
Definition: ptwXY.h:99
double biSectionMax
Definition: ptwXY.h:90
int64_t length
Definition: ptwXY.h:93
const G4int n
ptwXY_interpolationOtherInfo interpolationOtherInfo
Definition: ptwXY.h:88
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: ptwXY_core.cc:29
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
nfu_status status
Definition: ptwXY.h:85
double accuracy
Definition: ptwXY.h:91

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 25 of file ptwXY_binaryOperators.cc.

25  {
26 
27  int64_t i, nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
28  ptwXYPoint *p;
29  ptwXYOverflowPoint *o, *overflowHeader = &(ptwXY->overflowHeader);
30 
31  if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
32 
33  for( i = 0, p = ptwXY->points; i < nonOverflowLength; i++, p++ ) p->y = slope * p->y + offset;
34  for( o = overflowHeader->next; o != overflowHeader; o = o->next ) o->point.y = slope * o->point.y + offset;
35  return( ptwXY->status );
36 }
ptwXYPoint point
Definition: ptwXY.h:80
ptwXYPoint * points
Definition: ptwXY.h:99
const char * p
Definition: xmltok.h:285
int64_t ptwXY_getNonOverflowLength(ptwXYPoints const *ptwXY)
Definition: ptwXY_core.cc:590
double y
Definition: ptwXY.h:62
struct ptwXYOverflowPoint_s * next
Definition: ptwXY.h:78
nfu_status status
Definition: ptwXY.h:85
ptwXYOverflowPoint overflowHeader
Definition: ptwXY.h:98

Here is the call graph for this function:

Here is the caller graph for this function:

nfu_status ptwXY_sub_doubleFrom ( ptwXYPoints ptwXY,
double  value 
)

Definition at line 41 of file ptwXY_binaryOperators.cc.

41 { return( ptwXY_slopeOffset( ptwXY, 1., -value ) ); }
const XML_Char int const XML_Char * value
Definition: expat.h:331
nfu_status ptwXY_slopeOffset(ptwXYPoints *ptwXY, double slope, double offset)

Here is the call graph for this function:

nfu_status ptwXY_sub_fromDouble ( ptwXYPoints ptwXY,
double  value 
)

Definition at line 42 of file ptwXY_binaryOperators.cc.

42 { return( ptwXY_slopeOffset( ptwXY, -1., value ) ); }
const XML_Char int const XML_Char * value
Definition: expat.h:331
nfu_status ptwXY_slopeOffset(ptwXYPoints *ptwXY, double slope, double offset)

Here is the call graph for this function:

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

Definition at line 154 of file ptwXY_binaryOperators.cc.

154  {
155 
156  ptwXYPoints *diff;
157 
158  if( ptwXY1->length == 0 ) {
159  diff = ptwXY_clone( ptwXY2, status );
160  if( ( *status = ptwXY_neg( diff ) ) != nfu_Okay ) diff = ptwXY_free( diff ); }
161  else if( ptwXY2->length == 0 ) {
162  diff = ptwXY_clone( ptwXY1, status ); }
163  else {
164  diff = ptwXY_binary_ptwXY( ptwXY1, ptwXY2, 1., -1., 0., status );
165  }
166  return( diff );
167 }
nfu_status ptwXY_neg(ptwXYPoints *ptwXY)
ptwXYPoints * ptwXY_clone(ptwXYPoints *ptwXY, nfu_status *status)
Definition: ptwXY_core.cc:208
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
int64_t length
Definition: ptwXY.h:93
ptwXYPoints * ptwXY_binary_ptwXY(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, double v1, double v2, double v1v2, nfu_status *status)

Here is the call graph for this function:

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

Definition at line 144 of file ptwXY_methods.cc.

144  {
145 
146  double x1, x2 = 0., y1, y2 = 0., fx = 1.1, x, dx, dxp, lfx, y; /* fx initialized so compilers want complain. */
147  int64_t i, notFirstPass = 0;
148  int nfx, nDone, doLinear;
149  nfu_status status;
150 
152  if( ( sectionSubdivideMax < 1 ) || ( dxMax < 0. ) || ( fxMax < 1. ) ) return( nfu_badInput );
153  if( sectionSubdivideMax > ptwXY_sectionSubdivideMax ) sectionSubdivideMax = ptwXY_sectionSubdivideMax;
154  if( ( status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( status );
155  for( i = ptwXY1->length - 1; i >= 0; i-- ) {
156  x1 = ptwXY1->points[i].x;
157  y1 = ptwXY1->points[i].y;
158  if( notFirstPass ) {
159  dx = ptwXY_thicken_linear_dx( sectionSubdivideMax, dxMax, x1, x2 );
160 
161  if( x1 == 0. ) {
162  doLinear = 1; }
163  else {
164  fx = x2 / x1;
165  if( fx > 0. ) {
166  lfx = G4Log( fx );
167  if( fxMax == 1. ) {
168  nfx = sectionSubdivideMax; }
169  else {
170  nfx = ( (int) ( lfx / G4Log( fxMax ) ) ) + 1;
171  if( nfx > sectionSubdivideMax ) nfx = sectionSubdivideMax;
172  }
173  if( nfx > 0 ) fx = G4Exp( lfx / nfx );
174  doLinear = 0;
175  if( dx < ( fx - 1 ) * x1 ) doLinear = 1; }
176  else {
177  doLinear = 1;
178  }
179  }
180  x = x1;
181  dxp = dx;
182  nDone = 0;
183  while( 1 ) {
184  if( doLinear ) {
185  x += dx; }
186  else {
187  dx = ptwXY_thicken_linear_dx( sectionSubdivideMax - nDone, dxMax, x, x2 );
188  if( dx <= ( fx - 1 ) * x ) {
189  dxp = dx;
190  doLinear = 1;
191  continue;
192  }
193  dxp = ( fx - 1. ) * x;
194  x *= fx;
195  }
196  if( ( x2 - x ) < 0.05 * std::fabs( dxp ) ) break;
197  if( ( status = ptwXY_interpolatePoint( ptwXY1->interpolation, x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( status );
198  if( ( status = ptwXY_setValueAtX( ptwXY1, x, y ) ) != nfu_Okay ) return( status );
199  nDone++;
200  } // Loop checking, 11.06.2015, T. Koi
201  }
202  notFirstPass = 1;
203  x2 = x1;
204  y2 = y1;
205  }
206  return( status );
207 }
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
nfu_status ptwXY_interpolatePoint(ptwXY_interpolation interpolation, double x, double *y, double x1, double y1, double x2, double y2)
ptwXYPoint * points
Definition: ptwXY.h:99
static double ptwXY_thicken_linear_dx(int sectionSubdivideMax, double dxMax, double x1, double x2)
tuple x
Definition: test.py:50
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
Definition: ptwXY_core.cc:876
int64_t length
Definition: ptwXY.h:93
enum nfu_status_e nfu_status
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
#define ptwXY_sectionSubdivideMax
Definition: ptwXY.h:24

Here is the call graph for this function:

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

Definition at line 211 of file ptwXY_methods.cc.

211  {
212 
213  int64_t i, j, length = ptwXY1->length;
214  ptwXYPoints *thinned = NULL;
215  double y1, y2, y3;
216  char *thin = NULL;
217 
218  if( length < 3 ) return( ptwXY_clone( ptwXY1, status ) ); /* Logic below requires at least 2 points. */
219  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
220  *status = nfu_otherInterpolation;
221  if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
222  if( accuracy < ptwXY1->accuracy ) accuracy = ptwXY1->accuracy;
223  if( ( thinned = ptwXY_new( ptwXY1->interpolation, &(ptwXY1->interpolationOtherInfo),
224  ptwXY1->biSectionMax, accuracy, length, ptwXY1->overflowLength, status, ptwXY1->userFlag ) ) == NULL ) return( NULL );
225 
226  thinned->points[0] = ptwXY1->points[0]; /* This sections removes middle point if surrounding points have the same y-value. */
227  y1 = ptwXY1->points[0].y;
228  y2 = ptwXY1->points[1].y;
229  for( i = 2, j = 1; i < length; i++ ) {
230  y3 = ptwXY1->points[i].y;
231  if( ( y1 != y2 ) || ( y2 != y3 ) ) {
232  thinned->points[j++] = ptwXY1->points[i - 1];
233  y1 = y2;
234  y2 = y3;
235  }
236  }
237  thinned->points[j++] = ptwXY1->points[length - 1];
238 
239  if( ptwXY1->interpolation != ptwXY_interpolationFlat ) { /* Now call ptwXY_thin2 for more thinning. */
240  length = thinned->length = j;
241  if( ( thin = (char *) nfu_calloc( 1, (size_t) length ) ) == NULL ) goto Err;
242  if( ( *status = ptwXY_thin2( thinned, thin, accuracy, 0, length - 1 ) ) != nfu_Okay ) goto Err;
243  for( j = 1; j < length; j++ ) if( thin[j] != 0 ) break;
244  for( i = j + 1; i < length; i++ ) {
245  if( thin[i] == 0 ) {
246  thinned->points[j] = thinned->points[i];
247  j++;
248  }
249  }
250  nfu_free( thin );
251  }
252  thinned->length = j;
253 
254  return( thinned );
255 
256 Err:
257  ptwXY_free( thinned );
258  if( thin != NULL ) nfu_free( thin );
259  return( NULL );
260 }
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
int userFlag
Definition: ptwXY.h:89
ptwXYPoint * points
Definition: ptwXY.h:99
ptwXYPoints * ptwXY_clone(ptwXYPoints *ptwXY, nfu_status *status)
Definition: ptwXY_core.cc:208
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
double biSectionMax
Definition: ptwXY.h:90
int64_t length
Definition: ptwXY.h:93
void * nfu_calloc(size_t size, size_t n)
void * nfu_free(void *p)
ptwXY_interpolationOtherInfo interpolationOtherInfo
Definition: ptwXY.h:88
double y
Definition: ptwXY.h:62
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: ptwXY_core.cc:29
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
int64_t overflowLength
Definition: ptwXY.h:95
static nfu_status ptwXY_thin2(ptwXYPoints *thinned, char *thin, double accuracy, int64_t i1, int64_t i2)
double accuracy
Definition: ptwXY.h:91

Here is the call graph for this function:

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

Definition at line 153 of file ptwXY_interpolation.cc.

153  {
154 /*
155 * This function only works when 'ptwXY->interpolation == interpolationTo' or when interpolationTo is ptwXY_interpolationLinLin.
156 */
157  ptwXYPoints *n1;
158  interpolation_func func = NULL;
159 
160  if( ( *status = ptwXY->status ) != nfu_Okay ) return( NULL );
161  if( ptwXY->interpolation == interpolationTo ) {
162  *status = nfu_Okay;
163  return( ptwXY_clone( ptwXY, status ) ); }
164  else {
165  if( interpolationTo == ptwXY_interpolationLinLin ) {
166  switch( ptwXY->interpolation ) {
168  func = ptwXY_LogLogToLinLin; break;
170  func = ptwXY_LinLogToLinLin; break;
172  func = ptwXY_LogLinToLinLin; break;
174  if( ptwXY->interpolationOtherInfo.getValueFunc != NULL ) func = ptwXY_otherToLinLin;
175  break;
176  case ptwXY_interpolationLinLin : /* Stops compilers from complaining. */
178  break;
179  }
180  }
181  }
183  if( func == NULL ) return( NULL );
184 
185  *status = nfu_Okay;
186  if( ( n1 = ptwXY_cloneToInterpolation( ptwXY, interpolationTo, status ) ) == NULL ) return( NULL );
187  if( accuracy < ptwXY->accuracy ) accuracy = ptwXY->accuracy;
188  n1->accuracy = accuracy;
189 
192  *status = ptwXY_toOtherInterpolation2( n1, ptwXY, func );
194  n1->interpolationOtherInfo.argList = NULL;
195  if( *status != nfu_Okay ) n1 = ptwXY_free( n1 );
196  return( n1 );
197 }
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
ptwXYPoints * ptwXY_clone(ptwXYPoints *ptwXY, nfu_status *status)
Definition: ptwXY_core.cc:208
ptwXYPoints * ptwXY_cloneToInterpolation(ptwXYPoints *ptwXY, ptwXY_interpolation interpolationTo, nfu_status *status)
Definition: ptwXY_core.cc:215
ptwXY_getValue_callback getValueFunc
Definition: ptwXY.h:71
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
static nfu_status ptwXY_LogLinToLinLin(ptwXYPoints *desc, double x1, double y1, double x2, double y2, int depth)
static nfu_status ptwXY_LogLogToLinLin(ptwXYPoints *desc, double x1, double y1, double x2, double y2, int depth)
static nfu_status ptwXY_otherToLinLin(ptwXYPoints *desc, double x1, double y1, double x2, double y2, int depth)
ptwXY_interpolationOtherInfo interpolationOtherInfo
Definition: ptwXY.h:88
static nfu_status ptwXY_LinLogToLinLin(ptwXYPoints *desc, double x1, double y1, double x2, double y2, int depth)
nfu_status status
Definition: ptwXY.h:85
nfu_status(* interpolation_func)(ptwXYPoints *desc, double x1, double y1, double x2, double y2, int depth)
double accuracy
Definition: ptwXY.h:91
static nfu_status ptwXY_toOtherInterpolation2(ptwXYPoints *desc, ptwXYPoints *src, interpolation_func func)

Here is the call graph for this function:

ptwXYPoints* ptwXY_toUnitbase ( ptwXYPoints ptwXY,
nfu_status status 
)

Definition at line 306 of file ptwXY_interpolation.cc.

306  {
307 
308  int64_t i;
309  ptwXYPoints *n;
310  ptwXYPoint *p;
311  double xMin, xMax, dx, inverseDx;
312 
313  *status = nfu_tooFewPoints;
314  if( ptwXY->length < 2 ) return( NULL );
315  if( ( n = ptwXY_clone( ptwXY, status ) ) == NULL ) return( NULL );
316 
317  xMin = n->points[0].x;
318  xMax = n->points[n->length-1].x;
319  dx = xMax - xMin;
320  inverseDx = 1. / dx;
321  for( i = 0, p = n->points; i < n->length; i++, p++ ) {
322  p->x = ( p->x - xMin ) * inverseDx;
323  p->y = p->y * dx;
324  }
325  n->points[n->length-1].x = 1.; /* Make sure last point is realy 1. */
326  return( n );
327 }
ptwXYPoint * points
Definition: ptwXY.h:99
const char * p
Definition: xmltok.h:285
ptwXYPoints * ptwXY_clone(ptwXYPoints *ptwXY, nfu_status *status)
Definition: ptwXY_core.cc:208
int64_t length
Definition: ptwXY.h:93
const G4int n
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62

Here is the call graph for this function:

Here is the caller graph for this function:

nfu_status ptwXY_trim ( ptwXYPoints ptwXY)

Definition at line 315 of file ptwXY_methods.cc.

315  {
316 /*
317 c Remove extra zeros at beginning and end.
318 */
319 
320  int64_t i, i1, i2;
321  nfu_status status;
322 
323  if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
324  if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status );
325  for( i1 = 0; i1 < ptwXY->length; i1++ ) {
326  if( ptwXY->points[i1].y != 0 ) break;
327  }
328  if( i1 > 0 ) i1--;
329  for( i2 = ptwXY->length - 1; i2 >= 0; i2-- ) {
330  if( ptwXY->points[i2].y != 0 ) break;
331  }
332  i2++;
333  if( i2 < ptwXY->length ) i2++;
334  if( i2 > i1 ) {
335  if( i1 > 0 ) {
336  for( i = i1; i < i2; i++ ) ptwXY->points[i - i1] = ptwXY->points[i];
337  }
338  ptwXY->length = i2 - i1; }
339  else if( i2 < i1 ) { /* Remove all zeros between endpoints. */
340  ptwXY->points[1] = ptwXY->points[ptwXY->length - 1];
341  ptwXY->length = 2;
342  }
343 
344  return( nfu_Okay );
345 }
ptwXYPoint * points
Definition: ptwXY.h:99
int64_t length
Definition: ptwXY.h:93
enum nfu_status_e nfu_status
double y
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
nfu_status status
Definition: ptwXY.h:85

Here is the call graph for this function:

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

Definition at line 295 of file ptwXY_convenient.cc.

295  {
296 
297  nfu_status status;
298  int64_t n1 = ptwXY1->length, n2 = ptwXY2->length;
299  double sum, diff;
300  ptwXYPoint *xy1, *xy2;
301 
302  epsilon = std::fabs( epsilon ) + std::fabs( epsilonFactor * DBL_EPSILON );
303 
304  if( ( status = ptwXY1->status ) != nfu_Okay ) return( status );
305  if( ( status = ptwXY2->status ) != nfu_Okay ) return( status );
306  if( n1 == 0 ) return( nfu_empty );
307  if( n2 == 0 ) return( nfu_empty );
308  if( n1 < 2 ) {
309  status = nfu_tooFewPoints; }
310  else if( n2 < 2 ) {
311  status = nfu_tooFewPoints; }
312  else {
313  xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, 0 );
314  xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, 0 );
315  if( xy1->x < xy2->x ) {
316  if( xy2->y != 0. ) {
317  sum = std::fabs( xy1->x ) + std::fabs( xy2->x );
318  diff = std::fabs( xy2->x - xy1->x );
319  if( diff > epsilon * sum ) {
320  status = nfu_domainsNotMutual; }
321  else {
322  xy1->x = xy2->x;
323  }
324  } }
325  else if( xy1->x > xy2->x ) {
326  if( xy1->y != 0. ) {
327  sum = std::fabs( xy1->x ) + std::fabs( xy2->x );
328  diff = std::fabs( xy2->x - xy1->x );
329  if( diff > epsilon * sum ) {
330  status = nfu_domainsNotMutual; }
331  else {
332  xy2->x = xy1->x;
333  }
334  }
335  }
336 
337  if( status == nfu_Okay ) {
338  xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, n1 - 1 );
339  xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, n2 - 1 );
340  if( xy1->x < xy2->x ) {
341  if( xy1->y != 0. ) {
342  sum = std::fabs( xy1->x ) + std::fabs( xy2->x );
343  diff = std::fabs( xy2->x - xy1->x );
344  if( diff > epsilon * sum ) {
345  status = nfu_domainsNotMutual; }
346  else {
347  xy2->x = xy1->x;
348  }
349  } }
350  else if( xy1->x > xy2->x ) {
351  if( xy2->y != 0. ) {
352  sum = std::fabs( xy1->x ) + std::fabs( xy2->x );
353  diff = std::fabs( xy2->x - xy1->x );
354  if( diff > epsilon * sum ) {
355  status = nfu_domainsNotMutual; }
356  else {
357  xy1->x = xy2->x;
358  }
359  }
360  }
361  }
362  }
363  return( status );
364 }
ptwXYPoint * ptwXY_getPointAtIndex_Unsafely(ptwXYPoints *ptwXY, int64_t index)
Definition: ptwXY_core.cc:684
int64_t length
Definition: ptwXY.h:93
enum nfu_status_e nfu_status
#define DBL_EPSILON
Definition: templates.hh:87
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
nfu_status status
Definition: ptwXY.h:85
double epsilon(double density, double temperature)

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 349 of file ptwXY_methods.cc.

349  {
350 
351  int64_t overflowSize, i, i1 = 0, i2 = 0, n1 = ptwXY1->length, n2 = ptwXY2->length, length;
352  int fillWithFirst = unionOptions & ptwXY_union_fill, trim = unionOptions & ptwXY_union_trim;
353  ptwXYPoints *n;
354  double x1 = 0., x2 = 0., y1 = 0., y2 = 0., y, biSectionMax, accuracy;
355 
356  if( ( *status = ptwXY1->status ) != nfu_Okay ) return( NULL );
357  if( ( *status = ptwXY2->status ) != nfu_Okay ) return( NULL );
358  *status = nfu_otherInterpolation;
359  if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
360 /*
361 * Many other routines use the fact that ptwXY_union calls ptwXY_coalescePoints for ptwXY1 and ptwXY2 so do not change it.
362 */
363  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
364  if( ( *status = ptwXY_simpleCoalescePoints( ptwXY2 ) ) != nfu_Okay ) return( NULL );
365 
366  if( ( n1 == 1 ) || ( n2 == 1 ) ) {
367  *status = nfu_tooFewPoints;
368  return( NULL );
369  }
370  if( trim ) {
371  if( n1 > 0 ) {
372  if( n2 > 0 ) {
373  if( ptwXY1->points[0].x < ptwXY2->points[0].x ) {
374  while( i1 < n1 ) { // Loop checking, 11.05.2015, T. Koi
375  if( ptwXY1->points[i1].x >= ptwXY2->points[0].x ) break;
376  if( fillWithFirst ) {
377  if( i1 < ( ptwXY1->length - 1 ) ) {
378  x1 = ptwXY1->points[i1].x;
379  y1 = ptwXY1->points[i1].y;
380  x2 = ptwXY1->points[i1+1].x;
381  y2 = ptwXY1->points[i1+1].y;
382  }
383  }
384  i1++;
385  } }
386  else {
387  while( i2 < n2 ) { // Loop checking, 11.06.2015, T. Koi
388  if( ptwXY2->points[i2].x >= ptwXY1->points[0].x ) break;
389  i2++;
390  }
391  }
392  if( ptwXY1->points[n1-1].x > ptwXY2->points[n2-1].x ) {
393  while( i1 < n1 ) { // Loop checking, 11.06.2015, T. Koi
394  if( ptwXY1->points[n1-1].x <= ptwXY2->points[n2-1].x ) break;
395  n1--;
396  } }
397  else {
398  while( i2 < n2 ) { // Loop checking, 11.06.2015, T. Koi
399  if( ptwXY2->points[n2-1].x <= ptwXY1->points[n1-1].x ) break;
400  n2--;
401  }
402  } }
403  else {
404  n1 = 0;
405  } }
406  else {
407  n2 = 0;
408  }
409  }
410  overflowSize = ptwXY1->overflowAllocatedSize;
411  if( overflowSize < ptwXY2->overflowAllocatedSize ) overflowSize = ptwXY2->overflowAllocatedSize;
412  length = ( n1 - i1 ) + ( n2 - i2 );
413  if( length == 0 ) length = ptwXY_minimumSize;
414  biSectionMax = ptwXY1->biSectionMax;
415  if( biSectionMax < ptwXY2->biSectionMax ) biSectionMax = ptwXY2->biSectionMax;
416  accuracy = ptwXY1->accuracy;
417  if( accuracy < ptwXY2->accuracy ) accuracy = ptwXY2->accuracy;
418  n = ptwXY_new( ptwXY1->interpolation, NULL, biSectionMax, accuracy, length, overflowSize, status, ptwXY1->userFlag );
419  if( n == NULL ) return( NULL );
420 
421  for( i = 0; ( i1 < n1 ) && ( i2 < n2 ); i++ ) {
422  y = 0.;
423  if( ptwXY1->points[i1].x <= ptwXY2->points[i2].x ) {
424  n->points[i].x = ptwXY1->points[i1].x;
425  if( fillWithFirst ) {
426  y = ptwXY1->points[i1].y;
427  if( i1 < ( ptwXY1->length - 1 ) ) {
428  x1 = ptwXY1->points[i1].x;
429  y1 = ptwXY1->points[i1].y;
430  x2 = ptwXY1->points[i1+1].x;
431  y2 = ptwXY1->points[i1+1].y; }
432  else {
433  y1 = 0.;
434  y2 = 0.;
435  }
436  }
437  if( ptwXY1->points[i1].x == ptwXY2->points[i2].x ) i2++;
438  i1++; }
439  else {
440  n->points[i].x = ptwXY2->points[i2].x;
441  if( fillWithFirst && ( ( y1 != 0. ) || ( y2 != 0. ) ) ) {
442  if( ( *status = ptwXY_interpolatePoint( ptwXY1->interpolation, ptwXY2->points[i2].x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) {
443  ptwXY_free( n );
444  return( NULL );
445  }
446  }
447  i2++;
448  }
449  n->points[i].y = y;
450  }
451 
452  y = 0.;
453  for( ; i1 < n1; i1++, i++ ) {
454  n->points[i].x = ptwXY1->points[i1].x;
455  if( fillWithFirst ) y = ptwXY1->points[i1].y;
456  n->points[i].y = y;
457  }
458  for( ; i2 < n2; i2++, i++ ) {
459  n->points[i].x = ptwXY2->points[i2].x;
460  if( fillWithFirst && trim && ( n->points[i].x <= x2 ) ) {
461  if( ( *status = ptwXY_interpolatePoint( ptwXY1->interpolation, n->points[i].x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) {
462  ptwXY_free( n );
463  return( NULL );
464  }
465  }
466  n->points[i].y = y;
467  }
468  n->length = i;
469 
470  if( unionOptions & ptwXY_union_mergeClosePoints ) {
471  if( ( *status = ptwXY_mergeClosePoints( n, 4 * DBL_EPSILON ) ) != nfu_Okay ) {
472  ptwXY_free( n );
473  return( NULL );
474  }
475  }
476  return( n );
477 }
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
nfu_status ptwXY_interpolatePoint(ptwXY_interpolation interpolation, double x, double *y, double x1, double y1, double x2, double y2)
int userFlag
Definition: ptwXY.h:89
ptwXYPoint * points
Definition: ptwXY.h:99
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
double biSectionMax
Definition: ptwXY.h:90
int64_t length
Definition: ptwXY.h:93
#define ptwXY_union_fill
Definition: ptwXY.h:31
#define DBL_EPSILON
Definition: templates.hh:87
const G4int n
#define ptwXY_union_mergeClosePoints
Definition: ptwXY.h:33
double y
Definition: ptwXY.h:62
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: ptwXY_core.cc:29
double x
Definition: ptwXY.h:62
int64_t overflowAllocatedSize
Definition: ptwXY.h:96
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
nfu_status ptwXY_mergeClosePoints(ptwXYPoints *ptwXY, double epsilon)
nfu_status status
Definition: ptwXY.h:85
#define ptwXY_union_trim
Definition: ptwXY.h:32
double accuracy
Definition: ptwXY.h:91
#define ptwXY_minimumSize
Definition: ptwXY.h:20

Here is the call graph for this function:

Here is the caller graph for this function:

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.

363  {
364 /*
365 * Should we not be checking the interpolation members???????
366 */
367  int64_t i;
368  ptwXYPoints *n1 = NULL, *n2 = NULL, *a = NULL, *r = NULL;
369  ptwXYPoint *p;
370  double f, g, xMin, xMax;
371 
372  *status = nfu_XOutsideDomain;
373  if( w <= w1 ) {
374  if( w < w1 ) return( NULL );
375  return( ptwXY_clone( ptwXY1, status ) );
376  }
377  if( w >= w2 ) {
378  if( w > w2 ) return( NULL );
379  return( ptwXY_clone( ptwXY2, status ) );
380  }
381  if( ( n1 = ptwXY_toUnitbase( ptwXY1, status ) ) == NULL ) return( NULL );
382  if( ( n2 = ptwXY_toUnitbase( ptwXY2, status ) ) == NULL ) goto Err;
383  f = ( w - w1 ) / ( w2 - w1 );
384  g = 1. - f;
385  for( i = 0, p = n1->points; i < n1->length; i++, p++ ) p->y *= g;
386  for( i = 0, p = n2->points; i < n2->length; i++, p++ ) p->y *= f;
387  if( ( a = ptwXY_add_ptwXY( n1, n2, status ) ) == NULL ) goto Err;
388 
389  xMin = g * ptwXY1->points[0].x + f * ptwXY2->points[0].x;
390  xMax = g * ptwXY1->points[ptwXY1->length-1].x + f * ptwXY2->points[ptwXY2->length-1].x;
391  if( ( r = ptwXY_fromUnitbase( a, xMin, xMax, status ) ) == NULL ) goto Err;
392  ptwXY_free( n1 );
393  ptwXY_free( n2 );
394  ptwXY_free( a );
395  return( r );
396 
397 Err:
398  if( n1 != NULL ) ptwXY_free( n1 );
399  if( n2 != NULL ) ptwXY_free( n2 );
400  if( a != NULL ) ptwXY_free( a );
401  return( NULL );
402 }
ptwXYPoint * points
Definition: ptwXY.h:99
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
ptwXYPoints * ptwXY_fromUnitbase(ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status)
const char * p
Definition: xmltok.h:285
ptwXYPoints * ptwXY_clone(ptwXYPoints *ptwXY, nfu_status *status)
Definition: ptwXY_core.cc:208
ptwXYPoints * ptwXY_toUnitbase(ptwXYPoints *ptwXY, nfu_status *status)
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
int64_t length
Definition: ptwXY.h:93
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
ptwXYPoints * ptwXY_add_ptwXY(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status)

Here is the call graph for this function:

void ptwXY_update_biSectionMax ( ptwXYPoints ptwXY1,
double  oldLength 
)

Definition at line 31 of file ptwXY_misc.cc.

31  {
32 
33  ptwXY1->biSectionMax = ptwXY1->biSectionMax - 1.442695 * G4Log( ptwXY1->length / oldLength ); /* 1.442695 = 1 / std::log( 2. ) */
34  if( ptwXY1->biSectionMax < 0 ) ptwXY1->biSectionMax = 0;
36 }
double biSectionMax
Definition: ptwXY.h:90
int64_t length
Definition: ptwXY.h:93
G4double G4Log(G4double x)
Definition: G4Log.hh:230
#define ptwXY_maxBiSectionMax
Definition: ptwXY.h:22

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 450 of file ptwXY_convenient.cc.

450  {
451 
452  int64_t i1, length = ptwXY_length( ptwXY );
453  double *xps, *yps;
454  ptwXYPoint *pointFrom;
455  nfu_status status;
456 
457  if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
458  if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status );
459 
460  if( ( *xs = (double *) malloc( length * sizeof( double ) ) ) == NULL ) return( nfu_mallocError );
461  if( ( *ys = (double *) malloc( length * sizeof( double ) ) ) == NULL ) {
462  free( *xs );
463  *xs = NULL;
464  return( nfu_mallocError );
465  }
466 
467  for( i1 = 0, pointFrom = ptwXY->points, xps = *xs, yps = *ys; i1 < length; ++i1, ++pointFrom, ++xps, ++yps ) {
468  *xps = pointFrom->x;
469  *yps = pointFrom->y;
470  }
471 
472  return( nfu_Okay );
473 }
ptwXYPoint * points
Definition: ptwXY.h:99
int64_t ptwXY_length(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:583
enum nfu_status_e nfu_status
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
nfu_status status
Definition: ptwXY.h:85

Here is the call graph for this function:

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

Definition at line 477 of file ptwXY_convenient.cc.

477  {
478 
479  ptwXYPoints *n;
480 
481  *status = nfu_XNotAscending;
482  if( x1 >= x2 ) return( NULL );
483  *status = nfu_Okay;
484  if( ( n = ptwXY_new( ptwXY_interpolationLinLin, NULL, ptwXY_maxBiSectionMax, ptwXY_minAccuracy, 2, 0, status, 0 ) ) == NULL ) return( NULL );
485  ptwXY_setValueAtX( n, x1, y );
486  ptwXY_setValueAtX( n, x2, y );
487  return( n );
488 }
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
Definition: ptwXY_core.cc:876
const G4int n
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: ptwXY_core.cc:29
#define ptwXY_maxBiSectionMax
Definition: ptwXY.h:22
#define ptwXY_minAccuracy
Definition: ptwXY.h:23

Here is the call graph for this function:

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

Definition at line 326 of file ptwXY_core.cc.

326  {
327 
328  double xMin = 0.9 * xMax - 1;
329 
330  if( xMax < 0 ) xMin = 1.1 * xMax - 1;
331  if( ptwXY->length > 0 ) xMin = ptwXY_getXMin( ptwXY );
332  return( ptwXY_xSlice( ptwXY, xMin, xMax, secondarySize, fill, status ) );
333 }
int64_t length
Definition: ptwXY.h:93
double ptwXY_getXMin(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:1206
ptwXYPoints * ptwXY_xSlice(ptwXYPoints *ptwXY, double xMin, double xMax, int64_t secondarySize, int fill, nfu_status *status)
Definition: ptwXY_core.cc:274

Here is the call graph for this function:

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

Definition at line 315 of file ptwXY_core.cc.

315  {
316 
317  double xMax = 1.1 * xMin + 1;
318 
319  if( xMin < 0 ) xMax = 0.9 * xMin + 1;
320  if( ptwXY->length > 0 ) xMax = ptwXY_getXMax( ptwXY );
321  return( ptwXY_xSlice( ptwXY, xMin, xMax, secondarySize, fill, status ) );
322 }
int64_t length
Definition: ptwXY.h:93
ptwXYPoints * ptwXY_xSlice(ptwXYPoints *ptwXY, double xMin, double xMax, int64_t secondarySize, int fill, nfu_status *status)
Definition: ptwXY_core.cc:274
double ptwXY_getXMax(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:1239

Here is the call graph for this function:

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.

274  {
275 
276  int64_t i, i1, i2;
277  double y;
278  ptwXYPoints *n = NULL;
279 
280  if( ( *status = ptwXY->status ) != nfu_Okay ) return( NULL );
281 
282  if( ( ptwXY->length == 0 ) || ( ptwXY_getXMin( ptwXY ) >= xMax ) || ( ptwXY_getXMax( ptwXY ) <= xMin ) ) {
283  n = ptwXY_new( ptwXY->interpolation, &(ptwXY->interpolationOtherInfo), ptwXY->biSectionMax,
284  ptwXY->accuracy, 0, secondarySize, status, ptwXY->userFlag ); }
285  else {
286  if( ( n = ptwXY_clone( ptwXY, status ) ) == NULL ) return( NULL );
287  if( ( n->points[0].x < xMin ) || ( n->points[n->length - 1].x > xMax ) ) {
288  if( fill && ( n->points[n->length - 1].x > xMax ) ) {
289  if( ( *status = ptwXY_getValueAtX( n, xMax, &y ) ) != nfu_Okay ) goto Err;
290  if( ( *status = ptwXY_setValueAtX( n, xMax, y ) ) != nfu_Okay ) goto Err;
291  }
292  if( fill && ( n->points[0].x < xMin ) ) {
293  if( ( *status = ptwXY_getValueAtX( n, xMin, &y ) ) != nfu_Okay ) goto Err;
294  if( ( *status = ptwXY_setValueAtX( n, xMin, y ) ) != nfu_Okay ) goto Err;
295  }
296  ptwXY_coalescePoints( n, n->length + n->overflowAllocatedSize, NULL, 0 );
297  for( i1 = 0; i1 < n->length; i1++ ) if( n->points[i1].x >= xMin ) break;
298  for( i2 = n->length - 1; i2 > 0; i2-- ) if( n->points[i2].x <= xMax ) break;
299  i2++;
300  if( i1 > 0 ) {
301  for( i = i1; i < i2; i++ ) n->points[i- i1] = n->points[i];
302  }
303  n->length = i2 - i1;
304  }
305  }
306  return( n );
307 
308 Err:
309  if( n != NULL ) ptwXY_free( n );
310  return( NULL );
311 }
nfu_status ptwXY_getValueAtX(ptwXYPoints *ptwXY, double x, double *y)
Definition: ptwXY_core.cc:844
ptwXY_interpolation interpolation
Definition: ptwXY.h:87
nfu_status ptwXY_coalescePoints(ptwXYPoints *ptwXY, int64_t size, ptwXYPoint *newPoint, int forceSmallerResize)
Definition: ptwXY_core.cc:469
int userFlag
Definition: ptwXY.h:89
ptwXYPoint * points
Definition: ptwXY.h:99
ptwXYPoints * ptwXY_clone(ptwXYPoints *ptwXY, nfu_status *status)
Definition: ptwXY_core.cc:208
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
Definition: ptwXY_core.cc:876
double biSectionMax
Definition: ptwXY.h:90
int64_t length
Definition: ptwXY.h:93
double ptwXY_getXMin(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:1206
const G4int n
ptwXY_interpolationOtherInfo interpolationOtherInfo
Definition: ptwXY.h:88
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: ptwXY_core.cc:29
double x
Definition: ptwXY.h:62
int64_t overflowAllocatedSize
Definition: ptwXY.h:96
nfu_status status
Definition: ptwXY.h:85
double accuracy
Definition: ptwXY.h:91
double ptwXY_getXMax(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:1239

Here is the call graph for this function:

Here is the caller graph for this function: