12 #if defined __cplusplus
33 if( ( xArray =
ptwX_new( n, status ) ) == NULL )
return( NULL );
58 if( lowerEps != 0. ) {
59 if( std::fabs( lowerEps ) <
minEps ) {
61 if( lowerEps < 0. ) sign = -1;
73 dx = std::fabs( x1 * lowerEps );
74 if( x1 == 0 ) dx = std::fabs( lowerEps );
77 if( ( xp + dx ) < x2 ) {
87 if( ( xm < 0. ) && ( x1 >= 0. ) && positiveXOnly ) {
98 if( upperEps != 0. ) {
99 if( std::fabs( upperEps ) <
minEps ) {
101 if( upperEps < 0. ) sign = -1;
113 dx = std::fabs( x2 * upperEps );
114 if( x2 == 0 ) dx = std::fabs( upperEps );
117 if( ( xm - dx ) > x1 ) {
143 int64_t i, i1, j,
k,
n = ptwXY->
length;
147 if( n < 2 )
return( ptwXY->
status );
153 for( i1 = 1, p2++; i1 < ( n - 1 ); i1++, p2++ ) {
154 if( ( p2->
x - x ) > 0.5 * epsilon * ( std::fabs( x ) + std::fabs( p2->
x ) ) )
break;
157 for( i = i1, p1 = &(ptwXY->
points[1]); i <
n; i++, p1++, p2++ ) *p1 = *p2;
161 p1 = &(ptwXY->
points[n-1]);
163 for( i1 = n - 2, p1--; i1 > 0; i1--, p1-- ) {
164 if( x - p1->
x > 0.5 * epsilon * ( std::fabs( x ) + std::fabs( p1->
x ) ) )
break;
166 if( i1 != ( n - 2 ) ) {
168 n = ptwXY->
length = i1 + 2;
171 for( i = 1; i < n - 1; i++ ) {
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;
180 if( ( k = ( j - i ) ) > 1 ) {
183 for( p1 = &(ptwXY->
points[i+1]); j <
n; j++, p1++, p2++ ) *p1 = *p2;
197 double x,
y, xMin, xMax;
205 if( ( n =
ptwXY_clone( ptwXY, status ) ) == NULL )
return( NULL );
206 if( ptwXY->
length == 0 )
return( n );
210 if( ( xMin >= ptwX->
points[lengthX-1] ) || ( xMax <= ptwX->points[0] ) ) {
215 for( i = 0; i < lengthX; i++ ) {
217 if( x <= xMin )
continue;
218 if( x >= xMax )
break;
229 for( ; i1 < n->
length; i1++ ) {
230 if( n->
points[i1].
x == x )
break;
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;
244 for( i = i1; i < i2; i++ ) n->
points[i - i1] = n->
points[i];
274 if( xy1->
x < xy2->
x ) {
276 else if( xy1->
x > xy2->
x ) {
283 if( xy1->
x < xy2->
x ) {
285 else if( xy1->
x > xy2->
x ) {
302 epsilon = std::fabs( epsilon ) + std::fabs( epsilonFactor *
DBL_EPSILON );
315 if( xy1->
x < xy2->
x ) {
317 sum = std::fabs( xy1->
x ) + std::fabs( xy2->
x );
318 diff = std::fabs( xy2->
x - xy1->
x );
319 if( diff > epsilon * sum ) {
325 else if( xy1->
x > xy2->
x ) {
327 sum = std::fabs( xy1->
x ) + std::fabs( xy2->
x );
328 diff = std::fabs( xy2->
x - xy1->
x );
329 if( diff > epsilon * sum ) {
340 if( xy1->
x < xy2->
x ) {
342 sum = std::fabs( xy1->
x ) + std::fabs( xy2->
x );
343 diff = std::fabs( xy2->
x - xy1->
x );
344 if( diff > epsilon * sum ) {
350 else if( xy1->
x > xy2->
x ) {
352 sum = std::fabs( xy1->
x ) + std::fabs( xy2->
x );
353 diff = std::fabs( xy2->
x - xy1->
x );
354 if( diff > epsilon * sum ) {
369 ptwXYPoints *ptwXY2,
double lowerEps2,
double upperEps2,
int positiveXOnly2 ) {
392 if( xy1->
x < xy2->
x ) {
394 if( xy2->
y == 0. ) lowerEps2 = 0.; }
395 else if( xy1->
x > xy2->
x ) {
397 if( xy1->
y == 0. ) lowerEps1 = 0.; }
399 lowerEps1 = lowerEps2 = 0.;
404 if( xy1->
x < xy2->
x ) {
406 if( xy1->
y == 0. ) upperEps1 = 0.; }
407 else if( xy1->
x > xy2->
x ) {
409 if( xy2->
y == 0. ) upperEps2 = 0.; }
411 upperEps1 = upperEps2 = 0.;
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 );
434 if( index1 < 0 ) index1 = 0;
436 if( index2 < index1 ) index2 =
index1;
437 *numberOfPoints = index2 -
index1;
440 for( i = index1, pointFrom = ptwXY->
points; i < index2; i++, pointFrom++ ) {
441 *(d++) = pointFrom->
x;
442 *(d++) = pointFrom->
y;
460 if( ( *xs = (
double *) malloc( length *
sizeof(
double ) ) ) == NULL )
return(
nfu_mallocError );
461 if( ( *ys = (
double *) malloc( length *
sizeof(
double ) ) ) == NULL ) {
467 for( i1 = 0, pointFrom = ptwXY->
points, xps = *xs, yps = *ys; i1 < length; ++i1, ++pointFrom, ++xps, ++yps ) {
482 if( x1 >= x2 )
return( NULL );
496 double x1,
y1,
x2,
y2, accuracy2, yMin = 1
e-10;
499 if( accuracy < 1
e-5 ) accuracy = 1
e-5;
500 if( accuracy > 1
e-1 ) accuracy = 1
e-1;
502 accuracy2 = accuracy = gaussian->
accuracy;
503 if( accuracy2 > 5
e-3 ) accuracy2 = 5
e-3;
505 x1 = -std::sqrt( -2. *
G4Log( yMin ) );
508 y2 =
G4Exp( -0.5 * x2 * x2 );
510 gaussian->
accuracy = 20 * accuracy2;
515 y2 =
G4Exp( -0.5 * x2 * x2 );
521 y2 =
G4Exp( -0.5 * x2 * x2 );
527 y2 =
G4Exp( -0.5 * x2 * x2 );
534 for( i = 0, pm = pp - 2; i <
n; i++, pp++, pm-- ) {
538 gaussian->
length = 2 * n + 1;
553 double x = 0.5 * ( x1 +
x2 );
554 double y =
G4Exp( -x * x / 2 ), yMin = ( y1 * ( x2 -
x ) + y2 * ( x - x1 ) ) / ( x2 - x1 );
556 if( std::fabs( y - yMin ) > y * ptwXY->
accuracy ) morePoints = 1;
574 for( i = 0, point = gaussian->
points; i < gaussian->
length; i++, point++ ) {
575 point->
x = point->
x * sigma + xCenter;
576 point->
y *= amplitude;
579 if( ( sliced =
ptwXY_xSlice( gaussian, xMin, xMax, 10, 1, status ) ) == NULL )
goto Err;
591 #if defined __cplusplus