| Geant4 Cross Reference |
1 /* 1 /*
2 # <<BEGIN-copyright>> 2 # <<BEGIN-copyright>>
3 # <<END-copyright>> 3 # <<END-copyright>>
4 */ 4 */
5 5
6 #include <cmath> 6 #include <cmath>
7 #include <float.h> 7 #include <float.h>
8 8
9 #include "ptwXY.h" 9 #include "ptwXY.h"
10 10
11 #if defined __cplusplus 11 #if defined __cplusplus
12 #include <cmath> 12 #include <cmath>
13 #include "G4Exp.hh" 13 #include "G4Exp.hh"
14 #include "G4Log.hh" 14 #include "G4Log.hh"
15 namespace GIDI { 15 namespace GIDI {
16 using namespace GIDI; 16 using namespace GIDI;
17 #endif 17 #endif
18 18
19 static nfu_status ptwXY_clip2( ptwXYPoints *pt 19 static nfu_status ptwXY_clip2( ptwXYPoints *ptwXY1, double y, double x1, double y1, double x2, double y2 );
20 static double ptwXY_thicken_linear_dx( int sec 20 static double ptwXY_thicken_linear_dx( int sectionSubdivideMax, double dxMax, double x1, double x2 );
21 static nfu_status ptwXY_thin2( ptwXYPoints *th 21 static nfu_status ptwXY_thin2( ptwXYPoints *thinned, char *thin, double accuracy, int64_t i1, int64_t i2 );
22 /* 22 /*
23 ********************************************** 23 ************************************************************
24 */ 24 */
25 nfu_status ptwXY_clip( ptwXYPoints *ptwXY1, do 25 nfu_status ptwXY_clip( ptwXYPoints *ptwXY1, double yMin, double yMax ) {
26 /* 26 /*
27 This function acts oddly for xy = [ [ 1, 0 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 line 28 This function probably only works for linear, linear interpolation (mainly because of ptwXY_clip2).
29 */ 29 */
30 int64_t i, j, n; 30 int64_t i, j, n;
31 double x2, y2; 31 double x2, y2;
32 nfu_status status; 32 nfu_status status;
33 ptwXYPoints *clipped; 33 ptwXYPoints *clipped;
34 ptwXYPoint *points; 34 ptwXYPoint *points;
35 35
36 if( ( status = ptwXY_simpleCoalescePoints( 36 if( ( status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( status );
37 if( ptwXY1->interpolation == ptwXY_interpo 37 if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
38 n = ptwXY1->length; 38 n = ptwXY1->length;
39 if( n > 0 ) { 39 if( n > 0 ) {
40 i = 0; 40 i = 0;
41 if( ptwXY_getYMax( ptwXY1 ) < yMin ) i 41 if( ptwXY_getYMax( ptwXY1 ) < yMin ) i = 1;
42 if( ptwXY_getYMin( ptwXY1 ) > yMax ) i 42 if( ptwXY_getYMin( ptwXY1 ) > yMax ) i = 1;
43 if( i == 1 ) return( ptwXY_clear( ptwX 43 if( i == 1 ) return( ptwXY_clear( ptwXY1 ) );
44 } 44 }
45 if( n == 1 ) { 45 if( n == 1 ) {
46 y2 = ptwXY1->points[0].y; 46 y2 = ptwXY1->points[0].y;
47 if( y2 < yMin ) { 47 if( y2 < yMin ) {
48 ptwXY1->points[0].y = yMin; } 48 ptwXY1->points[0].y = yMin; }
49 else if( y2 > yMax ) { 49 else if( y2 > yMax ) {
50 ptwXY1->points[0].y = yMax; 50 ptwXY1->points[0].y = yMax;
51 } } 51 } }
52 else if( n > 1 ) { 52 else if( n > 1 ) {
53 if( ( clipped = ptwXY_new( ptwXY1->int 53 if( ( clipped = ptwXY_new( ptwXY1->interpolation, &(ptwXY1->interpolationOtherInfo),
54 ptwXY1->biSectionMax, ptwXY1-> 54 ptwXY1->biSectionMax, ptwXY1->accuracy, n, 10, &status, ptwXY1->userFlag ) ) == NULL )
55 return( ptwXY1->status = status ); 55 return( ptwXY1->status = status );
56 for( i = 0; i < n; i++ ) { 56 for( i = 0; i < n; i++ ) {
57 x2 = ptwXY1->points[i].x; 57 x2 = ptwXY1->points[i].x;
58 y2 = ptwXY1->points[i].y; 58 y2 = ptwXY1->points[i].y;
59 if( y2 < yMin ) { 59 if( y2 < yMin ) {
60 if( i > 0 ) { 60 if( i > 0 ) {
61 points = ptwXY_getPointAtI 61 points = ptwXY_getPointAtIndex_Unsafely( clipped, clipped->length - 1 );
62 if( points->y > yMin ) { 62 if( points->y > yMin ) {
63 if( ( status = ptwXY_c 63 if( ( status = ptwXY_clip2( clipped, yMin, points->x, points->y, x2, y2 ) ) != nfu_Okay ) goto Err;
64 } 64 }
65 } 65 }
66 if( ( status = ptwXY_setValueA 66 if( ( status = ptwXY_setValueAtX( clipped, x2, yMin ) ) != nfu_Okay ) goto Err;
67 j = i; 67 j = i;
68 for( i++; i < n; i++ ) if( !( 68 for( i++; i < n; i++ ) if( !( ptwXY1->points[i].y < yMin ) ) break;
69 if( i < n ) { 69 if( i < n ) {
70 x2 = ptwXY1->points[i].x; 70 x2 = ptwXY1->points[i].x;
71 y2 = ptwXY1->points[i].y; 71 y2 = ptwXY1->points[i].y;
72 if( ( status = ptwXY_clip2 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 ) { 73 if( y2 > yMax ) {
74 if( ( status = ptwXY_c 74 if( ( status = ptwXY_clip2( clipped, yMax, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) goto Err;
75 } } 75 } }
76 else if( j != n - 1 ) { 76 else if( j != n - 1 ) {
77 if( ( status = ptwXY_setVa 77 if( ( status = ptwXY_setValueAtX( clipped, ptwXY1->points[n - 1].x, yMin ) ) != nfu_Okay ) goto Err;
78 } 78 }
79 i--; } 79 i--; }
80 else if( y2 > yMax ) { 80 else if( y2 > yMax ) {
81 if( i > 0 ) { 81 if( i > 0 ) {
82 points = ptwXY_getPointAtI 82 points = ptwXY_getPointAtIndex_Unsafely( clipped, clipped->length - 1 );
83 if( points->y < yMax ) { 83 if( points->y < yMax ) {
84 if( ( status = ptwXY_c 84 if( ( status = ptwXY_clip2( clipped, yMax, points->x, points->y, x2, y2 ) ) != nfu_Okay ) goto Err;
85 } 85 }
86 } 86 }
87 if( ( status = ptwXY_setValueA 87 if( ( status = ptwXY_setValueAtX( clipped, x2, yMax ) ) != nfu_Okay ) goto Err;
88 j = i; 88 j = i;
89 for( i++; i < n; i++ ) if( !( 89 for( i++; i < n; i++ ) if( !( ptwXY1->points[i].y > yMax ) ) break;
90 if( i < n ) { 90 if( i < n ) {
91 x2 = ptwXY1->points[i].x; 91 x2 = ptwXY1->points[i].x;
92 y2 = ptwXY1->points[i].y; 92 y2 = ptwXY1->points[i].y;
93 if( ( status = ptwXY_clip2 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 ) { 94 if( y2 < yMin ) {
95 if( ( status = ptwXY_c 95 if( ( status = ptwXY_clip2( clipped, yMin, ptwXY1->points[i-1].x, ptwXY1->points[i-1].y, x2, y2 ) ) != nfu_Okay ) goto Err;
96 } } 96 } }
97 else if( j != n - 1 ) { 97 else if( j != n - 1 ) {
98 if( ( status = ptwXY_setVa 98 if( ( status = ptwXY_setValueAtX( clipped, ptwXY1->points[n - 1].x, yMax ) ) != nfu_Okay ) goto Err;
99 } 99 }
100 i--; } 100 i--; }
101 else { 101 else {
102 if( ( status = ptwXY_setValueA 102 if( ( status = ptwXY_setValueAtX( clipped, x2, y2 ) ) != nfu_Okay ) goto Err;
103 } 103 }
104 } 104 }
105 if( ( status = ptwXY_simpleCoalescePoi 105 if( ( status = ptwXY_simpleCoalescePoints( clipped ) ) != nfu_Okay ) goto Err;
106 ptwXY1->length = clipped->length; /* 106 ptwXY1->length = clipped->length; /* The squeamish may want to skip the next few lines. */
107 clipped->length = n; 107 clipped->length = n;
108 n = ptwXY1->allocatedSize; 108 n = ptwXY1->allocatedSize;
109 ptwXY1->allocatedSize = clipped->alloc 109 ptwXY1->allocatedSize = clipped->allocatedSize;
110 clipped->allocatedSize = n; 110 clipped->allocatedSize = n;
111 points = clipped->points; 111 points = clipped->points;
112 clipped->points = ptwXY1->points; 112 clipped->points = ptwXY1->points;
113 ptwXY1->points = points; 113 ptwXY1->points = points;
114 ptwXY_free( clipped ); 114 ptwXY_free( clipped );
115 } 115 }
116 116
117 return( ptwXY1->status ); 117 return( ptwXY1->status );
118 118
119 Err: 119 Err:
120 ptwXY_free( clipped ); 120 ptwXY_free( clipped );
121 return( ptwXY1->status = status ); 121 return( ptwXY1->status = status );
122 } 122 }
123 /* 123 /*
124 ********************************************** 124 ************************************************************
125 */ 125 */
126 static nfu_status ptwXY_clip2( ptwXYPoints *cl 126 static nfu_status ptwXY_clip2( ptwXYPoints *clipped, double y, double x1, double y1, double x2, double y2 ) {
127 127
128 double x; 128 double x;
129 nfu_status status = nfu_Okay; 129 nfu_status status = nfu_Okay;
130 130
131 x = ( y - y1 ) * ( x2 - x1 ) / ( y2 - y1 ) 131 x = ( y - y1 ) * ( x2 - x1 ) / ( y2 - y1 ) + x1;
132 if( x <= x1 ) { 132 if( x <= x1 ) {
133 x = x1; } 133 x = x1; }
134 else if( x >= x2 ) { 134 else if( x >= x2 ) {
135 x = x1; } 135 x = x1; }
136 else { 136 else {
137 status = ptwXY_setValueAtX( clipped, x 137 status = ptwXY_setValueAtX( clipped, x, y );
138 } 138 }
139 return( status ); 139 return( status );
140 } 140 }
141 /* 141 /*
142 ********************************************** 142 ************************************************************
143 */ 143 */
144 nfu_status ptwXY_thicken( ptwXYPoints *ptwXY1, 144 nfu_status ptwXY_thicken( ptwXYPoints *ptwXY1, int sectionSubdivideMax, double dxMax, double fxMax ) {
145 145
146 double x1, x2 = 0., y1, y2 = 0., fx = 1.1, 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; 147 int64_t i, notFirstPass = 0;
148 int nfx, nDone, doLinear; 148 int nfx, nDone, doLinear;
149 nfu_status status; 149 nfu_status status;
150 150
151 if( ptwXY1->interpolation == ptwXY_interpo 151 if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
152 if( ( sectionSubdivideMax < 1 ) || ( dxMax 152 if( ( sectionSubdivideMax < 1 ) || ( dxMax < 0. ) || ( fxMax < 1. ) ) return( nfu_badInput );
153 if( sectionSubdivideMax > ptwXY_sectionSub 153 if( sectionSubdivideMax > ptwXY_sectionSubdivideMax ) sectionSubdivideMax = ptwXY_sectionSubdivideMax;
154 if( ( status = ptwXY_simpleCoalescePoints( 154 if( ( status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( status );
155 for( i = ptwXY1->length - 1; i >= 0; i-- ) 155 for( i = ptwXY1->length - 1; i >= 0; i-- ) {
156 x1 = ptwXY1->points[i].x; 156 x1 = ptwXY1->points[i].x;
157 y1 = ptwXY1->points[i].y; 157 y1 = ptwXY1->points[i].y;
158 if( notFirstPass ) { 158 if( notFirstPass ) {
159 dx = ptwXY_thicken_linear_dx( sect 159 dx = ptwXY_thicken_linear_dx( sectionSubdivideMax, dxMax, x1, x2 );
160 160
161 if( x1 == 0. ) { 161 if( x1 == 0. ) {
162 doLinear = 1; } 162 doLinear = 1; }
163 else { 163 else {
164 fx = x2 / x1; 164 fx = x2 / x1;
165 if( fx > 0. ) { 165 if( fx > 0. ) {
166 lfx = G4Log( fx ); 166 lfx = G4Log( fx );
167 if( fxMax == 1. ) { 167 if( fxMax == 1. ) {
168 nfx = sectionSubdivide 168 nfx = sectionSubdivideMax; }
169 else { 169 else {
170 nfx = ( (int) ( lfx / 170 nfx = ( (int) ( lfx / G4Log( fxMax ) ) ) + 1;
171 if( nfx > sectionSubdi 171 if( nfx > sectionSubdivideMax ) nfx = sectionSubdivideMax;
172 } 172 }
173 if( nfx > 0 ) fx = G4Exp( 173 if( nfx > 0 ) fx = G4Exp( lfx / nfx );
174 doLinear = 0; 174 doLinear = 0;
175 if( dx < ( fx - 1 ) * x1 ) 175 if( dx < ( fx - 1 ) * x1 ) doLinear = 1; }
176 else { 176 else {
177 doLinear = 1; 177 doLinear = 1;
178 } 178 }
179 } 179 }
180 x = x1; 180 x = x1;
181 dxp = dx; 181 dxp = dx;
182 nDone = 0; 182 nDone = 0;
183 while( 1 ) { 183 while( 1 ) {
184 if( doLinear ) { 184 if( doLinear ) {
185 x += dx; } 185 x += dx; }
186 else { 186 else {
187 dx = ptwXY_thicken_linear_ 187 dx = ptwXY_thicken_linear_dx( sectionSubdivideMax - nDone, dxMax, x, x2 );
188 if( dx <= ( fx - 1 ) * x ) 188 if( dx <= ( fx - 1 ) * x ) {
189 dxp = dx; 189 dxp = dx;
190 doLinear = 1; 190 doLinear = 1;
191 continue; 191 continue;
192 } 192 }
193 dxp = ( fx - 1. ) * x; 193 dxp = ( fx - 1. ) * x;
194 x *= fx; 194 x *= fx;
195 } 195 }
196 if( ( x2 - x ) < 0.05 * std::f 196 if( ( x2 - x ) < 0.05 * std::fabs( dxp ) ) break;
197 if( ( status = ptwXY_interpola 197 if( ( status = ptwXY_interpolatePoint( ptwXY1->interpolation, x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( status );
198 if( ( status = ptwXY_setValueA 198 if( ( status = ptwXY_setValueAtX( ptwXY1, x, y ) ) != nfu_Okay ) return( status );
199 nDone++; 199 nDone++;
200 } // Loop checking, 11.06.2015, T. 200 } // Loop checking, 11.06.2015, T. Koi
201 } 201 }
202 notFirstPass = 1; 202 notFirstPass = 1;
203 x2 = x1; 203 x2 = x1;
204 y2 = y1; 204 y2 = y1;
205 } 205 }
206 return( status ); 206 return( status );
207 } 207 }
208 /* 208 /*
209 ********************************************** 209 ************************************************************
210 */ 210 */
211 ptwXYPoints *ptwXY_thin( ptwXYPoints *ptwXY1, 211 ptwXYPoints *ptwXY_thin( ptwXYPoints *ptwXY1, double accuracy, nfu_status *status ) {
212 212
213 int64_t i, j, length = ptwXY1->length; 213 int64_t i, j, length = ptwXY1->length;
214 ptwXYPoints *thinned = NULL; 214 ptwXYPoints *thinned = NULL;
215 double y1, y2, y3; 215 double y1, y2, y3;
216 char *thin = NULL; 216 char *thin = NULL;
217 217
218 if( length < 3 ) return( ptwXY_clone( ptwX 218 if( length < 3 ) return( ptwXY_clone( ptwXY1, status ) ); /* Logic below requires at least 2 points. */
219 if( ( *status = ptwXY_simpleCoalescePoints 219 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
220 *status = nfu_otherInterpolation; 220 *status = nfu_otherInterpolation;
221 if( ptwXY1->interpolation == ptwXY_interpo 221 if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
222 if( accuracy < ptwXY1->accuracy ) accuracy 222 if( accuracy < ptwXY1->accuracy ) accuracy = ptwXY1->accuracy;
223 if( ( thinned = ptwXY_new( ptwXY1->interpo 223 if( ( thinned = ptwXY_new( ptwXY1->interpolation, &(ptwXY1->interpolationOtherInfo),
224 ptwXY1->biSectionMax, accuracy, length 224 ptwXY1->biSectionMax, accuracy, length, ptwXY1->overflowLength, status, ptwXY1->userFlag ) ) == NULL ) return( NULL );
225 225
226 thinned->points[0] = ptwXY1->points[0]; 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; 227 y1 = ptwXY1->points[0].y;
228 y2 = ptwXY1->points[1].y; 228 y2 = ptwXY1->points[1].y;
229 for( i = 2, j = 1; i < length; i++ ) { 229 for( i = 2, j = 1; i < length; i++ ) {
230 y3 = ptwXY1->points[i].y; 230 y3 = ptwXY1->points[i].y;
231 if( ( y1 != y2 ) || ( y2 != y3 ) ) { 231 if( ( y1 != y2 ) || ( y2 != y3 ) ) {
232 thinned->points[j++] = ptwXY1->poi 232 thinned->points[j++] = ptwXY1->points[i - 1];
233 y1 = y2; 233 y1 = y2;
234 y2 = y3; 234 y2 = y3;
235 } 235 }
236 } 236 }
237 thinned->points[j++] = ptwXY1->points[leng 237 thinned->points[j++] = ptwXY1->points[length - 1];
238 238
239 if( ptwXY1->interpolation != ptwXY_interpo 239 if( ptwXY1->interpolation != ptwXY_interpolationFlat ) { /* Now call ptwXY_thin2 for more thinning. */
240 length = thinned->length = j; 240 length = thinned->length = j;
241 if( ( thin = (char *) nfu_calloc( 1, ( 241 if( ( thin = (char *) nfu_calloc( 1, (size_t) length ) ) == NULL ) goto Err;
242 if( ( *status = ptwXY_thin2( thinned, 242 if( ( *status = ptwXY_thin2( thinned, thin, accuracy, 0, length - 1 ) ) != nfu_Okay ) goto Err;
243 for( j = 1; j < length; j++ ) if( thin 243 for( j = 1; j < length; j++ ) if( thin[j] != 0 ) break;
244 for( i = j + 1; i < length; i++ ) { 244 for( i = j + 1; i < length; i++ ) {
245 if( thin[i] == 0 ) { 245 if( thin[i] == 0 ) {
246 thinned->points[j] = thinned-> 246 thinned->points[j] = thinned->points[i];
247 j++; 247 j++;
248 } 248 }
249 } 249 }
250 nfu_free( thin ); 250 nfu_free( thin );
251 } 251 }
252 thinned->length = j; 252 thinned->length = j;
253 253
254 return( thinned ); 254 return( thinned );
255 255
256 Err: 256 Err:
257 ptwXY_free( thinned ); 257 ptwXY_free( thinned );
258 if( thin != NULL ) nfu_free( thin ); 258 if( thin != NULL ) nfu_free( thin );
259 return( NULL ); 259 return( NULL );
260 } 260 }
261 /* 261 /*
262 ********************************************** 262 ************************************************************
263 */ 263 */
264 static nfu_status ptwXY_thin2( ptwXYPoints *th 264 static nfu_status ptwXY_thin2( ptwXYPoints *thinned, char *thin, double accuracy, int64_t i1, int64_t i2 ) {
265 265
266 int64_t i, iMax = 0; 266 int64_t i, iMax = 0;
267 double y, s, dY, dYMax = 0., dYR, dYRMax = 267 double y, s, dY, dYMax = 0., dYR, dYRMax = 0;
268 double x1 = thinned->points[i1].x, y1 = th 268 double x1 = thinned->points[i1].x, y1 = thinned->points[i1].y, x2 = thinned->points[i2].x, y2 = thinned->points[i2].y;
269 nfu_status status = nfu_Okay; 269 nfu_status status = nfu_Okay;
270 270
271 if( i1 + 1 >= i2 ) return( nfu_Okay ); 271 if( i1 + 1 >= i2 ) return( nfu_Okay );
272 for( i = i1 + 1; i < i2; i++ ) { 272 for( i = i1 + 1; i < i2; i++ ) {
273 if( ( status = ptwXY_interpolatePoint( 273 if( ( status = ptwXY_interpolatePoint( thinned->interpolation, thinned->points[i].x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( status );
274 s = 0.5 * ( std::fabs( y ) + std::fabs 274 s = 0.5 * ( std::fabs( y ) + std::fabs( thinned->points[i].y ) );
275 dY = std::fabs( y - thinned->points[i] 275 dY = std::fabs( y - thinned->points[i].y );
276 dYR = 0; 276 dYR = 0;
277 if( s != 0 ) dYR = dY / s; 277 if( s != 0 ) dYR = dY / s;
278 if( ( dYR > dYRMax ) || ( ( dYR >= 0.9 278 if( ( dYR > dYRMax ) || ( ( dYR >= 0.9999 * dYRMax ) && ( dY > dYMax ) ) ) { /* The choice of 0.9999 is not exact science. */
279 iMax = i; 279 iMax = i;
280 if( dY > dYMax ) dYMax = dY; 280 if( dY > dYMax ) dYMax = dY;
281 if( dYR > dYRMax ) dYRMax = dYR; 281 if( dYR > dYRMax ) dYRMax = dYR;
282 } 282 }
283 } 283 }
284 if( dYRMax < accuracy ) { 284 if( dYRMax < accuracy ) {
285 for( i = i1 + 1; i < i2; i++ ) thin[i] 285 for( i = i1 + 1; i < i2; i++ ) thin[i] = 1; }
286 else { 286 else {
287 if( ( status = ptwXY_thin2( thinned, t 287 if( ( status = ptwXY_thin2( thinned, thin, accuracy, i1, iMax ) ) != nfu_Okay ) return( status );
288 status = ptwXY_thin2( thinned, thin, a 288 status = ptwXY_thin2( thinned, thin, accuracy, iMax, i2 );
289 } 289 }
290 return( status ); 290 return( status );
291 } 291 }
292 /* 292 /*
293 ********************************************** 293 ************************************************************
294 */ 294 */
295 static double ptwXY_thicken_linear_dx( int sec 295 static double ptwXY_thicken_linear_dx( int sectionSubdivideMax, double dxMax, double x1, double x2 ) {
296 296
297 int ndx; 297 int ndx;
298 double dx = x2 - x1, dndx; 298 double dx = x2 - x1, dndx;
299 299
300 if( dxMax == 0. ) { 300 if( dxMax == 0. ) {
301 dx = ( x2 - x1 ) / sectionSubdivideMax 301 dx = ( x2 - x1 ) / sectionSubdivideMax; }
302 else { 302 else {
303 dndx = dx / dxMax; 303 dndx = dx / dxMax;
304 ndx = (int) dndx; 304 ndx = (int) dndx;
305 if( ( dndx - ndx ) > 1e-6 ) ndx++; 305 if( ( dndx - ndx ) > 1e-6 ) ndx++;
306 if( ndx > sectionSubdivideMax ) ndx = 306 if( ndx > sectionSubdivideMax ) ndx = sectionSubdivideMax;
307 if( ndx > 0 ) dx /= ndx; 307 if( ndx > 0 ) dx /= ndx;
308 } 308 }
309 309
310 return( dx ); 310 return( dx );
311 } 311 }
312 /* 312 /*
313 ********************************************** 313 ************************************************************
314 */ 314 */
315 nfu_status ptwXY_trim( ptwXYPoints *ptwXY ) { 315 nfu_status ptwXY_trim( ptwXYPoints *ptwXY ) {
316 /* 316 /*
317 c Remove extra zeros at beginning and end. 317 c Remove extra zeros at beginning and end.
318 */ 318 */
319 319
320 int64_t i, i1, i2; 320 int64_t i, i1, i2;
321 nfu_status status; 321 nfu_status status;
322 322
323 if( ptwXY->status != nfu_Okay ) return( pt 323 if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
324 if( ( status = ptwXY_simpleCoalescePoints( 324 if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status );
325 for( i1 = 0; i1 < ptwXY->length; i1++ ) { 325 for( i1 = 0; i1 < ptwXY->length; i1++ ) {
326 if( ptwXY->points[i1].y != 0 ) break; 326 if( ptwXY->points[i1].y != 0 ) break;
327 } 327 }
328 if( i1 > 0 ) i1--; 328 if( i1 > 0 ) i1--;
329 for( i2 = ptwXY->length - 1; i2 >= 0; i2-- 329 for( i2 = ptwXY->length - 1; i2 >= 0; i2-- ) {
330 if( ptwXY->points[i2].y != 0 ) break; 330 if( ptwXY->points[i2].y != 0 ) break;
331 } 331 }
332 i2++; 332 i2++;
333 if( i2 < ptwXY->length ) i2++; 333 if( i2 < ptwXY->length ) i2++;
334 if( i2 > i1 ) { 334 if( i2 > i1 ) {
335 if( i1 > 0 ) { 335 if( i1 > 0 ) {
336 for( i = i1; i < i2; i++ ) ptwXY-> 336 for( i = i1; i < i2; i++ ) ptwXY->points[i - i1] = ptwXY->points[i];
337 } 337 }
338 ptwXY->length = i2 - i1; } 338 ptwXY->length = i2 - i1; }
339 else if( i2 < i1 ) { /* Rem 339 else if( i2 < i1 ) { /* Remove all zeros between endpoints. */
340 ptwXY->points[1] = ptwXY->points[ptwXY 340 ptwXY->points[1] = ptwXY->points[ptwXY->length - 1];
341 ptwXY->length = 2; 341 ptwXY->length = 2;
342 } 342 }
343 343
344 return( nfu_Okay ); 344 return( nfu_Okay );
345 } 345 }
346 /* 346 /*
347 ********************************************** 347 ************************************************************
348 */ 348 */
349 ptwXYPoints *ptwXY_union( ptwXYPoints *ptwXY1, 349 ptwXYPoints *ptwXY_union( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int unionOptions ) {
350 350
351 int64_t overflowSize, i, i1 = 0, i2 = 0, n 351 int64_t overflowSize, i, i1 = 0, i2 = 0, n1 = ptwXY1->length, n2 = ptwXY2->length, length;
352 int fillWithFirst = unionOptions & ptwXY_u 352 int fillWithFirst = unionOptions & ptwXY_union_fill, trim = unionOptions & ptwXY_union_trim;
353 ptwXYPoints *n; 353 ptwXYPoints *n;
354 double x1 = 0., x2 = 0., y1 = 0., y2 = 0., 354 double x1 = 0., x2 = 0., y1 = 0., y2 = 0., y, biSectionMax, accuracy;
355 355
356 if( ( *status = ptwXY1->status ) != nfu_Ok 356 if( ( *status = ptwXY1->status ) != nfu_Okay ) return( NULL );
357 if( ( *status = ptwXY2->status ) != nfu_Ok 357 if( ( *status = ptwXY2->status ) != nfu_Okay ) return( NULL );
358 *status = nfu_otherInterpolation; 358 *status = nfu_otherInterpolation;
359 if( ptwXY1->interpolation == ptwXY_interpo 359 if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
360 /* 360 /*
361 * Many other routines use the fact that ptwX 361 * Many other routines use the fact that ptwXY_union calls ptwXY_coalescePoints for ptwXY1 and ptwXY2 so do not change it.
362 */ 362 */
363 if( ( *status = ptwXY_simpleCoalescePoints 363 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
364 if( ( *status = ptwXY_simpleCoalescePoints 364 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY2 ) ) != nfu_Okay ) return( NULL );
365 365
366 if( ( n1 == 1 ) || ( n2 == 1 ) ) { 366 if( ( n1 == 1 ) || ( n2 == 1 ) ) {
367 *status = nfu_tooFewPoints; 367 *status = nfu_tooFewPoints;
368 return( NULL ); 368 return( NULL );
369 } 369 }
370 if( trim ) { 370 if( trim ) {
371 if( n1 > 0 ) { 371 if( n1 > 0 ) {
372 if( n2 > 0 ) { 372 if( n2 > 0 ) {
373 if( ptwXY1->points[0].x < ptwX 373 if( ptwXY1->points[0].x < ptwXY2->points[0].x ) {
374 while( i1 < n1 ) { // Loop 374 while( i1 < n1 ) { // Loop checking, 11.05.2015, T. Koi
375 if( ptwXY1->points[i1] 375 if( ptwXY1->points[i1].x >= ptwXY2->points[0].x ) break;
376 if( fillWithFirst ) { 376 if( fillWithFirst ) {
377 if( i1 < ( ptwXY1- 377 if( i1 < ( ptwXY1->length - 1 ) ) {
378 x1 = ptwXY1->p 378 x1 = ptwXY1->points[i1].x;
379 y1 = ptwXY1->p 379 y1 = ptwXY1->points[i1].y;
380 x2 = ptwXY1->p 380 x2 = ptwXY1->points[i1+1].x;
381 y2 = ptwXY1->p 381 y2 = ptwXY1->points[i1+1].y;
382 } 382 }
383 } 383 }
384 i1++; 384 i1++;
385 } } 385 } }
386 else { 386 else {
387 while( i2 < n2 ) { // Loop 387 while( i2 < n2 ) { // Loop checking, 11.06.2015, T. Koi
388 if( ptwXY2->points[i2] 388 if( ptwXY2->points[i2].x >= ptwXY1->points[0].x ) break;
389 i2++; 389 i2++;
390 } 390 }
391 } 391 }
392 if( ptwXY1->points[n1-1].x > p 392 if( ptwXY1->points[n1-1].x > ptwXY2->points[n2-1].x ) {
393 while( i1 < n1 ) { // Loop 393 while( i1 < n1 ) { // Loop checking, 11.06.2015, T. Koi
394 if( ptwXY1->points[n1- 394 if( ptwXY1->points[n1-1].x <= ptwXY2->points[n2-1].x ) break;
395 n1--; 395 n1--;
396 } } 396 } }
397 else { 397 else {
398 while( i2 < n2 ) { // Loop 398 while( i2 < n2 ) { // Loop checking, 11.06.2015, T. Koi
399 if( ptwXY2->points[n2- 399 if( ptwXY2->points[n2-1].x <= ptwXY1->points[n1-1].x ) break;
400 n2--; 400 n2--;
401 } 401 }
402 } } 402 } }
403 else { 403 else {
404 n1 = 0; 404 n1 = 0;
405 } } 405 } }
406 else { 406 else {
407 n2 = 0; 407 n2 = 0;
408 } 408 }
409 } 409 }
410 overflowSize = ptwXY1->overflowAllocatedSi 410 overflowSize = ptwXY1->overflowAllocatedSize;
411 if( overflowSize < ptwXY2->overflowAllocat 411 if( overflowSize < ptwXY2->overflowAllocatedSize ) overflowSize = ptwXY2->overflowAllocatedSize;
412 length = ( n1 - i1 ) + ( n2 - i2 ); 412 length = ( n1 - i1 ) + ( n2 - i2 );
413 if( length == 0 ) length = ptwXY_minimumSi 413 if( length == 0 ) length = ptwXY_minimumSize;
414 biSectionMax = ptwXY1->biSectionMax; 414 biSectionMax = ptwXY1->biSectionMax;
415 if( biSectionMax < ptwXY2->biSectionMax ) 415 if( biSectionMax < ptwXY2->biSectionMax ) biSectionMax = ptwXY2->biSectionMax;
416 accuracy = ptwXY1->accuracy; 416 accuracy = ptwXY1->accuracy;
417 if( accuracy < ptwXY2->accuracy ) accuracy 417 if( accuracy < ptwXY2->accuracy ) accuracy = ptwXY2->accuracy;
418 n = ptwXY_new( ptwXY1->interpolation, NULL 418 n = ptwXY_new( ptwXY1->interpolation, NULL, biSectionMax, accuracy, length, overflowSize, status, ptwXY1->userFlag );
419 if( n == NULL ) return( NULL ); 419 if( n == NULL ) return( NULL );
420 420
421 for( i = 0; ( i1 < n1 ) && ( i2 < n2 ); i+ 421 for( i = 0; ( i1 < n1 ) && ( i2 < n2 ); i++ ) {
422 y = 0.; 422 y = 0.;
423 if( ptwXY1->points[i1].x <= ptwXY2->po 423 if( ptwXY1->points[i1].x <= ptwXY2->points[i2].x ) {
424 n->points[i].x = ptwXY1->points[i1 424 n->points[i].x = ptwXY1->points[i1].x;
425 if( fillWithFirst ) { 425 if( fillWithFirst ) {
426 y = ptwXY1->points[i1].y; 426 y = ptwXY1->points[i1].y;
427 if( i1 < ( ptwXY1->length - 1 427 if( i1 < ( ptwXY1->length - 1 ) ) {
428 x1 = ptwXY1->points[i1].x; 428 x1 = ptwXY1->points[i1].x;
429 y1 = ptwXY1->points[i1].y; 429 y1 = ptwXY1->points[i1].y;
430 x2 = ptwXY1->points[i1+1]. 430 x2 = ptwXY1->points[i1+1].x;
431 y2 = ptwXY1->points[i1+1]. 431 y2 = ptwXY1->points[i1+1].y; }
432 else { 432 else {
433 y1 = 0.; 433 y1 = 0.;
434 y2 = 0.; 434 y2 = 0.;
435 } 435 }
436 } 436 }
437 if( ptwXY1->points[i1].x == ptwXY2 437 if( ptwXY1->points[i1].x == ptwXY2->points[i2].x ) i2++;
438 i1++; } 438 i1++; }
439 else { 439 else {
440 n->points[i].x = ptwXY2->points[i2 440 n->points[i].x = ptwXY2->points[i2].x;
441 if( fillWithFirst && ( ( y1 != 0. 441 if( fillWithFirst && ( ( y1 != 0. ) || ( y2 != 0. ) ) ) {
442 if( ( *status = ptwXY_interpol 442 if( ( *status = ptwXY_interpolatePoint( ptwXY1->interpolation, ptwXY2->points[i2].x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) {
443 ptwXY_free( n ); 443 ptwXY_free( n );
444 return( NULL ); 444 return( NULL );
445 } 445 }
446 } 446 }
447 i2++; 447 i2++;
448 } 448 }
449 n->points[i].y = y; 449 n->points[i].y = y;
450 } 450 }
451 451
452 y = 0.; 452 y = 0.;
453 for( ; i1 < n1; i1++, i++ ) { 453 for( ; i1 < n1; i1++, i++ ) {
454 n->points[i].x = ptwXY1->points[i1].x; 454 n->points[i].x = ptwXY1->points[i1].x;
455 if( fillWithFirst ) y = ptwXY1->points 455 if( fillWithFirst ) y = ptwXY1->points[i1].y;
456 n->points[i].y = y; 456 n->points[i].y = y;
457 } 457 }
458 for( ; i2 < n2; i2++, i++ ) { 458 for( ; i2 < n2; i2++, i++ ) {
459 n->points[i].x = ptwXY2->points[i2].x; 459 n->points[i].x = ptwXY2->points[i2].x;
460 if( fillWithFirst && trim && ( n->poin 460 if( fillWithFirst && trim && ( n->points[i].x <= x2 ) ) {
461 if( ( *status = ptwXY_interpolateP 461 if( ( *status = ptwXY_interpolatePoint( ptwXY1->interpolation, n->points[i].x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) {
462 ptwXY_free( n ); 462 ptwXY_free( n );
463 return( NULL ); 463 return( NULL );
464 } 464 }
465 } 465 }
466 n->points[i].y = y; 466 n->points[i].y = y;
467 } 467 }
468 n->length = i; 468 n->length = i;
469 469
470 if( unionOptions & ptwXY_union_mergeCloseP 470 if( unionOptions & ptwXY_union_mergeClosePoints ) {
471 if( ( *status = ptwXY_mergeClosePoints 471 if( ( *status = ptwXY_mergeClosePoints( n, 4 * DBL_EPSILON ) ) != nfu_Okay ) {
472 ptwXY_free( n ); 472 ptwXY_free( n );
473 return( NULL ); 473 return( NULL );
474 } 474 }
475 } 475 }
476 return( n ); 476 return( n );
477 } 477 }
478 /* 478 /*
479 ********************************************** 479 ************************************************************
480 */ 480 */
481 nfu_status ptwXY_scaleOffsetXAndY( ptwXYPoints 481 nfu_status ptwXY_scaleOffsetXAndY( ptwXYPoints *ptwXY, double xScale, double xOffset, double yScale, double yOffset ) {
482 482
483 int64_t i1, length = ptwXY->length; 483 int64_t i1, length = ptwXY->length;
484 ptwXYPoint *p1; 484 ptwXYPoint *p1;
485 nfu_status status; 485 nfu_status status;
486 486
487 if( ptwXY->status != nfu_Okay ) return( pt 487 if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
488 if( xScale == 0 ) return( nfu_XNotAscendin 488 if( xScale == 0 ) return( nfu_XNotAscending );
489 489
490 if( ( status = ptwXY_simpleCoalescePoints( 490 if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status );
491 491
492 for( i1 = 0, p1 = ptwXY->points; i1 < leng 492 for( i1 = 0, p1 = ptwXY->points; i1 < length; i1++, p1++ ) {
493 p1->x = xScale * p1->x + xOffset; 493 p1->x = xScale * p1->x + xOffset;
494 p1->y = yScale * p1->y + yOffset; 494 p1->y = yScale * p1->y + yOffset;
495 } 495 }
496 496
497 if( xScale < 0 ) { 497 if( xScale < 0 ) {
498 int64_t length_2 = length / 2; 498 int64_t length_2 = length / 2;
499 ptwXYPoint tmp, *p2; 499 ptwXYPoint tmp, *p2;
500 500
501 for( i1 = 0, p1 = ptwXY->points, p2 = 501 for( i1 = 0, p1 = ptwXY->points, p2 = &(ptwXY->points[length-1]); i1 < length_2; i1++ ) {
502 tmp = *p1; 502 tmp = *p1;
503 *p1 = *p2; 503 *p1 = *p2;
504 *p2 = tmp; 504 *p2 = tmp;
505 } 505 }
506 } 506 }
507 507
508 return( ptwXY->status ); 508 return( ptwXY->status );
509 } 509 }
510 510
511 #if defined __cplusplus 511 #if defined __cplusplus
512 } 512 }
513 #endif 513 #endif
514 514