Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
ptwXY_convenient.cc File Reference
#include <stdlib.h>
#include <cmath>
#include <float.h>
#include "ptwXY.h"

Go to the source code of this file.

Macros

#define minEps   5e-16
 

Functions

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 *xys)
 
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, nfu_status *status)
 

Macro Definition Documentation

◆ minEps

#define minEps   5e-16

Function Documentation

◆ ptwXY_areDomainsMutual()

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}
@ nfu_domainsNotMutual
Definition: nf_utilities.h:28
@ nfu_Okay
Definition: nf_utilities.h:25
@ nfu_tooFewPoints
Definition: nf_utilities.h:28
@ nfu_empty
Definition: nf_utilities.h:28
enum nfu_status_e nfu_status
ptwXYPoint * ptwXY_getPointAtIndex_Unsafely(ptwXYPoints *ptwXY, int64_t index)
Definition: ptwXY_core.cc:684
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
int64_t length
Definition: ptwXY.h:93
nfu_status status
Definition: ptwXY.h:85

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

◆ ptwXY_copyToC_XY()

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

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}
@ nfu_insufficientMemory
Definition: nf_utilities.h:25
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
ptwXYPoint * points
Definition: ptwXY.h:99

◆ ptwXY_createGaussian()

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

Definition at line 566 of file ptwXY_convenient.cc.

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
586Err:
587 ptwXY_free( gaussian );
588 return( NULL );
589}
ptwXYPoints * ptwXY_xSlice(ptwXYPoints *ptwXY, double xMin, double xMax, int64_t secondarySize, int fill, nfu_status *status)
Definition: ptwXY_core.cc:274
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
ptwXYPoints * ptwXY_createGaussianCenteredSigma1(double accuracy, nfu_status *status)

◆ ptwXY_createGaussianCenteredSigma1()

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
542Err:
543 ptwXY_free( gaussian );
544 return( NULL );
545}
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
G4double G4Log(G4double x)
Definition: G4Log.hh:227
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
Definition: ptwXY_core.cc:876
nfu_status ptwXY_coalescePoints(ptwXYPoints *ptwXY, int64_t size, ptwXYPoint *newPoint, int forceSmallerResize)
Definition: ptwXY_core.cc:469
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
@ ptwXY_interpolationLinLin
Definition: ptwXY.h:35
double accuracy
Definition: ptwXY.h:91

Referenced by ptwXY_createGaussian().

◆ ptwXY_dullEdges()

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.
51This 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}
G4int sign(const T t)
@ nfu_invalidInterpolation
Definition: nf_utilities.h:27
@ nfu_otherInterpolation
Definition: nf_utilities.h:29
nfu_status ptwXY_getValueAtX(ptwXYPoints *ptwXY, double x, double *y)
Definition: ptwXY_core.cc:844
@ ptwXY_interpolationFlat
Definition: ptwXY.h:36
@ ptwXY_interpolationOther
Definition: ptwXY.h:36
nfu_status ptwXY_interpolatePoint(ptwXY_interpolation interpolation, double x, double *y, double x1, double y1, double x2, double y2)
#define minEps
ptwXY_interpolation interpolation
Definition: ptwXY.h:87

Referenced by MCGIDI_reaction_fixDomains(), and ptwXY_mutualifyDomains().

◆ ptwXY_getXArray()

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}
ptwXPoints * ptwX_new(int64_t size, nfu_status *status)
Definition: ptwX_core.cc:24
int64_t length
Definition: ptwX.h:26
double * points
Definition: ptwX.h:29

◆ ptwXY_intersectionWith_ptwX()

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
250Err:
251 ptwXY_free( n );
252 return( NULL );
253}
ptwXYPoints * ptwXY_clone(ptwXYPoints *ptwXY, nfu_status *status)
Definition: ptwXY_core.cc:208
int64_t ptwX_length(ptwXPoints *ptwX)
Definition: ptwX_core.cc:166
nfu_status status
Definition: ptwX.h:25

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

◆ ptwXY_mergeClosePoints()

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}
G4double epsilon(G4double density, G4double temperature)
#define DBL_EPSILON
Definition: templates.hh:66

Referenced by ptwXY_union().

◆ ptwXY_mutualifyDomains()

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

Definition at line 368 of file ptwXY_convenient.cc.

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 );
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}
nfu_status ptwXY_dullEdges(ptwXYPoints *ptwXY, double lowerEps, double upperEps, int positiveXOnly)
nfu_status ptwXY_areDomainsMutual(ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2)

◆ ptwXY_tweakDomainsToMutualify()

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}

Referenced by ptwXY_groupThreeFunctions(), and ptwXY_groupTwoFunctions().

◆ ptwXY_valueTo_ptwXAndY()

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}
@ nfu_mallocError
Definition: nf_utilities.h:25
int64_t ptwXY_length(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:583

◆ ptwXY_valueTo_ptwXY()

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

Definition at line 477 of file ptwXY_convenient.cc.

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_XNotAscending
Definition: nf_utilities.h:26
#define ptwXY_minAccuracy
Definition: ptwXY.h:23
#define ptwXY_maxBiSectionMax
Definition: ptwXY.h:22