Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ptwXY_binaryOperators.cc File Reference
#include <cmath>
#include <float.h>
#include "ptwXY.h"
Include dependency graph for ptwXY_binaryOperators.cc:

Go to the source code of this file.

Functions

static double ptwXY_mod2 (double v, double m, int pythonMod)
 
static nfu_status ptwXY_mul2_s_ptwXY (ptwXYPoints *n, ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, double x1, double y1, double x2, double y2, int level)
 
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)
 
static ptwXYPointsptwXY_div_ptwXY_forFlats (ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int safeDivide)
 
static nfu_status ptwXY_getValueAtX_ignore_XOutsideDomainError (ptwXYPoints *ptwXY1, double x, double *y)
 
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)
 

Function Documentation

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:

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_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:

static ptwXYPoints * ptwXY_div_ptwXY_forFlats ( ptwXYPoints ptwXY1,
ptwXYPoints ptwXY2,
nfu_status status,
int  safeDivide 
)
static

Definition at line 480 of file ptwXY_binaryOperators.cc.

480  {
481 
482  int64_t i;
483  ptwXYPoints *n = NULL;
484  ptwXYPoint *p;
485  double y;
486 
487  *status = nfu_invalidInterpolation;
488  if( ptwXY1->interpolation != ptwXY_interpolationFlat ) return( NULL );
489  if( ptwXY2->interpolation != ptwXY_interpolationFlat ) return( NULL );
490  if( ( n = ptwXY_union( ptwXY1, ptwXY2, status, ptwXY_union_fill | ptwXY_union_mergeClosePoints ) ) != NULL ) {
491  for( i = 0, p = n->points; i < n->length; i++, p++ ) {
492  if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, p->x, &y ) ) != nfu_Okay ) goto Err;
493  if( y == 0. ) {
494  if( ( safeDivide ) && ( p->y == 0 ) ) {
495  *status = nfu_divByZero;
496  goto Err;
497  } }
498  else {
499  p->y /= y;
500  }
501  }
502  }
503  return( n );
504 
505 Err:
506  if( n ) ptwXY_free( n );
507  return( NULL );
508 }
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
#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:

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 
)
static

Definition at line 432 of file ptwXY_binaryOperators.cc.

433  {
434 
435  nfu_status status;
436  double u1, u2, v1, v2, v, x, y, yp, dx, a1, a2;
437 
438  if( ( x2 - x1 ) < ClosestAllowXFactor * DBL_EPSILON * ( std::fabs( x1 ) + std::fabs( x2 ) ) ) return( nfu_Okay );
439  if( level >= n->biSectionMax ) return( nfu_Okay );
440  level++;
441  if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x1, &u1 ) ) != nfu_Okay ) return( status );
442  if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x2, &u2 ) ) != nfu_Okay ) return( status );
443  if( ( status = ptwXY_getValueAtX( ptwXY2, x1, &v1 ) ) != nfu_Okay ) return( status );
444  if( ( status = ptwXY_getValueAtX( ptwXY2, x2, &v2 ) ) != nfu_Okay ) return( status );
445  if( isNAN1 ) {
446  x = 0.5 * ( x1 + x2 );
447  if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x, &u1 ) ) != nfu_Okay ) return( status );
448  if( ( status = ptwXY_getValueAtX( ptwXY2, x, &v1 ) ) != nfu_Okay ) return( status );
449  y = u1 / v1; }
450  else if( isNAN2 ) {
451  x = 0.5 * ( x1 + x2 );
452  if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x, &u2 ) ) != nfu_Okay ) return( status );
453  if( ( status = ptwXY_getValueAtX( ptwXY2, x, &v2 ) ) != nfu_Okay ) return( status );
454  y = u2 / v2; }
455  else {
456  if( ( u1 == u2 ) || ( v1 == v2 ) ) return( nfu_Okay );
457  if( ( y1 == 0. ) || ( y2 == 0. ) ) { /* Handle special case of 0 where accuracy can never be met. */
458  x = 0.5 * ( x1 + x2 ); }
459  else {
460  if( ( u1 * u2 < 0. ) ) return( nfu_Okay ); /* Assume rounding error and no point needed. */
461  a1 = std::sqrt( std::fabs( u1 ) );
462  a2 = std::sqrt( std::fabs( u2 ) );
463  x = ( a2 * x1 + a1 * x2 ) / ( a2 + a1 );
464  }
465  dx = x2 - x1;
466  v = v1 * ( x2 - x ) + v2 * ( x - x1 );
467  if( ( v1 == 0. ) || ( v2 == 0. ) || ( v == 0. ) ) return( nfu_Okay ); /* Probably not correct, but I had to do something. */
468  yp = ( u1 / v1 * ( x2 - x ) + u2 / v2 * ( x - x1 ) ) / dx;
469  y = ( u1 * ( x2 - x ) + u2 * ( x - x1 ) ) / v;
470  if( std::fabs( y - yp ) < std::fabs( y * n->accuracy ) ) return( nfu_Okay );
471  }
472  if( ( status = ptwXY_setValueAtX( n, x, y ) ) != nfu_Okay ) return( status );
473  if( ( status = ptwXY_div_s_ptwXY( n, ptwXY1, ptwXY2, x, y, x2, y2, level, 0, isNAN2 ) ) != nfu_Okay ) return( status );
474  status = ptwXY_div_s_ptwXY( n, ptwXY1, ptwXY2, x1, y1, x, y, level, isNAN1, 0 );
475  return( status );
476 }
#define ClosestAllowXFactor
Definition: ptwXY.h:25
nfu_status ptwXY_getValueAtX(ptwXYPoints *ptwXY, double x, double *y)
Definition: ptwXY_core.cc:844
tuple x
Definition: test.py:50
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
Definition: ptwXY_core.cc:876
double biSectionMax
Definition: ptwXY.h:90
enum nfu_status_e nfu_status
#define DBL_EPSILON
Definition: templates.hh:87
static nfu_status ptwXY_getValueAtX_ignore_XOutsideDomainError(ptwXYPoints *ptwXY1, double x, double *y)
tuple v
Definition: test.py:18
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)
double accuracy
Definition: ptwXY.h:91

Here is the call graph for this function:

Here is the caller graph for this function:

static nfu_status ptwXY_getValueAtX_ignore_XOutsideDomainError ( ptwXYPoints ptwXY1,
double  x,
double *  y 
)
static

Definition at line 512 of file ptwXY_binaryOperators.cc.

512  {
513 
514  nfu_status status = ptwXY_getValueAtX( ptwXY1, x, y );
515 
516  if( status == nfu_XOutsideDomain ) status = nfu_Okay;
517  return( status );
518 }
nfu_status ptwXY_getValueAtX(ptwXYPoints *ptwXY, double x, double *y)
Definition: ptwXY_core.cc:844
tuple x
Definition: test.py:50
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_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:

static double ptwXY_mod2 ( double  v,
double  m,
int  pythonMod 
)
static

Definition at line 92 of file ptwXY_binaryOperators.cc.

92  {
93 
94  double r = std::fmod( std::fabs( v ), std::fabs( m ) );
95 
96  if( pythonMod ) {
97  if( ( v * m ) < 0. ) r = std::fabs( m ) - std::fabs( r );
98  if( m < 0. ) r *= -1.; }
99  else {
100  if( v < 0. ) r *= -1.;
101  }
102 
103  return( r );
104 }
static constexpr double m
Definition: G4SIunits.hh:129
tuple v
Definition: test.py:18

Here is the caller 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:

static nfu_status ptwXY_mul2_s_ptwXY ( ptwXYPoints n,
ptwXYPoints ptwXY1,
ptwXYPoints ptwXY2,
double  x1,
double  y1,
double  x2,
double  y2,
int  level 
)
static

Definition at line 251 of file ptwXY_binaryOperators.cc.

251  {
252 
253  nfu_status status;
254  double u1, u2, v1, v2, x, y, yp, dx, a1, a2;
255 
256  if( ( x2 - x1 ) < ClosestAllowXFactor * DBL_EPSILON * ( std::fabs( x1 ) + std::fabs( x2 ) ) ) return( nfu_Okay );
257  if( level >= n->biSectionMax ) return( nfu_Okay );
258  level++;
259  if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x1, &u1 ) ) != nfu_Okay ) return( status );
260  if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x2, &u2 ) ) != nfu_Okay ) return( status );
261  if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, x1, &v1 ) ) != nfu_Okay ) return( status );
262  if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, x2, &v2 ) ) != nfu_Okay ) return( status );
263  if( ( u1 == u2 ) || ( v1 == v2 ) ) return( nfu_Okay );
264  a1 = u1 * v1;
265  if( y1 == 0 ) a1 = 0.; /* Fix rounding problem. */
266  a2 = u2 * v2;
267  if( y2 == 0 ) a2 = 0.; /* Fix rounding problem. */
268  if( ( a1 == 0. ) || ( a2 == 0. ) ) { /* Handle special case of 0 where accuracy can never be met. */
269  x = 0.5 * ( x1 + x2 ); }
270  else {
271  if( ( a1 * a2 < 0. ) ) return( nfu_Okay ); /* Assume rounding error and no point needed as zero crossings are set in ptwXY_mul2_ptwXY. */
272  a1 = std::sqrt( std::fabs( a1 ) );
273  a2 = std::sqrt( std::fabs( a2 ) );
274  x = ( a2 * x1 + a1 * x2 ) / ( a2 + a1 );
275  }
276  dx = x2 - x1;
277  yp = ( u1 * v1 * ( x2 - x ) + u2 * v2 * ( x - x1 ) ) / dx;
278  y = ( u1 * ( x2 - x ) + u2 * ( x - x1 ) ) * ( v1 * ( x2 - x ) + v2 * ( x - x1 ) ) / ( dx * dx );
279  if( std::fabs( y - yp ) < std::fabs( y * n->accuracy ) ) return( nfu_Okay );
280  if( ( status = ptwXY_setValueAtX( n, x, y ) ) != nfu_Okay ) return( status );
281  if( ( status = ptwXY_mul2_s_ptwXY( n, ptwXY1, ptwXY2, x, y, x2, y2, level ) ) != nfu_Okay ) return( status );
282  status = ptwXY_mul2_s_ptwXY( n, ptwXY1, ptwXY2, x1, y1, x, y, level );
283  return( status );
284 }
#define ClosestAllowXFactor
Definition: ptwXY.h:25
tuple x
Definition: test.py:50
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
Definition: ptwXY_core.cc:876
double biSectionMax
Definition: ptwXY.h:90
enum nfu_status_e nfu_status
#define DBL_EPSILON
Definition: templates.hh:87
static nfu_status ptwXY_getValueAtX_ignore_XOutsideDomainError(ptwXYPoints *ptwXY1, double x, double *y)
double accuracy
Definition: ptwXY.h:91
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:

Here is the caller 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_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: