| Geant4 Cross Reference |
1 /* 1 /*
2 # <<BEGIN-copyright>> 2 # <<BEGIN-copyright>>
3 # <<END-copyright>> 3 # <<END-copyright>>
4 */ 4 */
5 5
6 #include <stdlib.h> 6 #include <stdlib.h>
7 #include <cmath> 7 #include <cmath>
8 #include <float.h> 8 #include <float.h>
9 9
10 #include "ptwXY.h" 10 #include "ptwXY.h"
11 11
12 #if defined __cplusplus 12 #if defined __cplusplus
13 #include <cmath> 13 #include <cmath>
14 #include "G4Exp.hh" 14 #include "G4Exp.hh"
15 #include "G4Log.hh" 15 #include "G4Log.hh"
16 namespace GIDI { 16 namespace GIDI {
17 using namespace GIDI; 17 using namespace GIDI;
18 #endif 18 #endif
19 19
20 static nfu_status ptwXY_createGaussianCentered 20 static nfu_status ptwXY_createGaussianCenteredSigma1_2( ptwXYPoints *ptwXY, double x1, double y1, double x2, double y2, int addX1Point );
21 /* 21 /*
22 ********************************************** 22 ************************************************************
23 */ 23 */
24 ptwXPoints *ptwXY_getXArray( ptwXYPoints *ptwX 24 ptwXPoints *ptwXY_getXArray( ptwXYPoints *ptwXY, nfu_status *status ) {
25 25
26 int64_t i, n; 26 int64_t i, n;
27 ptwXPoints *xArray; 27 ptwXPoints *xArray;
28 28
29 if( ( *status = ptwXY->status ) != nfu_Oka 29 if( ( *status = ptwXY->status ) != nfu_Okay ) return( NULL );
30 n = ptwXY->length; 30 n = ptwXY->length;
31 31
32 if( ( *status = ptwXY_simpleCoalescePoints 32 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( NULL );
33 if( ( xArray = ptwX_new( n, status ) ) == 33 if( ( xArray = ptwX_new( n, status ) ) == NULL ) return( NULL );
34 for( i = 0; i < n; i++ ) xArray->points[i] 34 for( i = 0; i < n; i++ ) xArray->points[i] = ptwXY->points[i].x;
35 xArray->length = n; 35 xArray->length = n;
36 36
37 return( xArray ); 37 return( xArray );
38 } 38 }
39 /* 39 /*
40 ********************************************** 40 ************************************************************
41 */ 41 */
42 nfu_status ptwXY_dullEdges( ptwXYPoints *ptwXY 42 nfu_status ptwXY_dullEdges( ptwXYPoints *ptwXY, double lowerEps, double upperEps, int positiveXOnly ) {
43 43
44 #define minEps 5e-16 44 #define minEps 5e-16
45 45
46 nfu_status status; 46 nfu_status status;
47 double xm, xp, dx, y, x1, y1, x2, y2, sign 47 double xm, xp, dx, y, x1, y1, x2, y2, sign;
48 ptwXYPoint *p; 48 ptwXYPoint *p;
49 49
50 /* This routine can only be used for linear in 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. */ 51 This needs to be fixed and documented. */
52 if( ( status = ptwXY->status ) != nfu_Okay 52 if( ( status = ptwXY->status ) != nfu_Okay ) return( status );
53 if( ptwXY->interpolation == ptwXY_interpol 53 if( ptwXY->interpolation == ptwXY_interpolationFlat ) return( nfu_invalidInterpolation );
54 if( ptwXY->interpolation == ptwXY_interpol 54 if( ptwXY->interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
55 55
56 if( ptwXY->length < 2 ) return( nfu_Okay ) 56 if( ptwXY->length < 2 ) return( nfu_Okay );
57 57
58 if( lowerEps != 0. ) { 58 if( lowerEps != 0. ) {
59 if( std::fabs( lowerEps ) < minEps ) { 59 if( std::fabs( lowerEps ) < minEps ) {
60 sign = 1; 60 sign = 1;
61 if( lowerEps < 0. ) sign = -1; 61 if( lowerEps < 0. ) sign = -1;
62 lowerEps = sign * minEps; 62 lowerEps = sign * minEps;
63 } 63 }
64 64
65 p = ptwXY_getPointAtIndex_Unsafely( pt 65 p = ptwXY_getPointAtIndex_Unsafely( ptwXY, 0 );
66 x1 = p->x; 66 x1 = p->x;
67 y1 = p->y; 67 y1 = p->y;
68 p = ptwXY_getPointAtIndex_Unsafely( pt 68 p = ptwXY_getPointAtIndex_Unsafely( ptwXY, 1 );
69 x2 = p->x; 69 x2 = p->x;
70 y2 = p->y; 70 y2 = p->y;
71 71
72 if( y1 != 0. ) { 72 if( y1 != 0. ) {
73 dx = std::fabs( x1 * lowerEps ); 73 dx = std::fabs( x1 * lowerEps );
74 if( x1 == 0 ) dx = std::fabs( lowe 74 if( x1 == 0 ) dx = std::fabs( lowerEps );
75 xm = x1 - dx; 75 xm = x1 - dx;
76 xp = x1 + dx; 76 xp = x1 + dx;
77 if( ( xp + dx ) < x2 ) { 77 if( ( xp + dx ) < x2 ) {
78 if( ( status = ptwXY_getValueA 78 if( ( status = ptwXY_getValueAtX( ptwXY, xp, &y ) ) != nfu_Okay ) return( status );
79 if( ( status = ptwXY_setValueA 79 if( ( status = ptwXY_setValueAtX( ptwXY, xp, y ) ) != nfu_Okay ) return( status ); }
80 else { 80 else {
81 xp = x2; 81 xp = x2;
82 y = y2; 82 y = y2;
83 } 83 }
84 if( lowerEps > 0 ) { 84 if( lowerEps > 0 ) {
85 if( ( status = ptwXY_setValueA 85 if( ( status = ptwXY_setValueAtX( ptwXY, x1, 0. ) ) != nfu_Okay ) return( status ); }
86 else { 86 else {
87 if( ( xm < 0. ) && ( x1 >= 0. 87 if( ( xm < 0. ) && ( x1 >= 0. ) && positiveXOnly ) {
88 if( ( status = ptwXY_setVa 88 if( ( status = ptwXY_setValueAtX( ptwXY, x1, 0. ) ) != nfu_Okay ) return( status ); }
89 else { 89 else {
90 if( ( status = ptwXY_setVa 90 if( ( status = ptwXY_setValueAtX( ptwXY, xm, 0. ) ) != nfu_Okay ) return( status );
91 if( ( status = ptwXY_inter 91 if( ( status = ptwXY_interpolatePoint( ptwXY->interpolation, x1, &y, xm, 0., xp, y ) ) != nfu_Okay ) return( status );
92 if( ( status = ptwXY_setVa 92 if( ( status = ptwXY_setValueAtX( ptwXY, x1, y ) ) != nfu_Okay ) return( status );
93 } 93 }
94 } 94 }
95 } 95 }
96 } 96 }
97 97
98 if( upperEps != 0. ) { 98 if( upperEps != 0. ) {
99 if( std::fabs( upperEps ) < minEps ) { 99 if( std::fabs( upperEps ) < minEps ) {
100 sign = 1; 100 sign = 1;
101 if( upperEps < 0. ) sign = -1; 101 if( upperEps < 0. ) sign = -1;
102 upperEps = sign * minEps; 102 upperEps = sign * minEps;
103 } 103 }
104 104
105 p = ptwXY_getPointAtIndex_Unsafely( pt 105 p = ptwXY_getPointAtIndex_Unsafely( ptwXY, ptwXY->length - 2 );
106 x1 = p->x; 106 x1 = p->x;
107 y1 = p->y; 107 y1 = p->y;
108 p = ptwXY_getPointAtIndex_Unsafely( pt 108 p = ptwXY_getPointAtIndex_Unsafely( ptwXY, ptwXY->length - 1 );
109 x2 = p->x; 109 x2 = p->x;
110 y2 = p->y; 110 y2 = p->y;
111 111
112 if( y2 != 0. ) { 112 if( y2 != 0. ) {
113 dx = std::fabs( x2 * upperEps ); 113 dx = std::fabs( x2 * upperEps );
114 if( x2 == 0 ) dx = std::fabs( uppe 114 if( x2 == 0 ) dx = std::fabs( upperEps );
115 xm = x2 - dx; 115 xm = x2 - dx;
116 xp = x2 + dx; 116 xp = x2 + dx;
117 if( ( xm - dx ) > x1 ) { 117 if( ( xm - dx ) > x1 ) {
118 if( ( status = ptwXY_getValueA 118 if( ( status = ptwXY_getValueAtX( ptwXY, xm, &y ) ) != nfu_Okay ) return( status );
119 if( ( status = ptwXY_setValueA 119 if( ( status = ptwXY_setValueAtX( ptwXY, xm, y ) ) != nfu_Okay ) return( status ); }
120 else { 120 else {
121 xm = x1; 121 xm = x1;
122 y = y1; 122 y = y1;
123 } 123 }
124 if( upperEps < 0 ) { 124 if( upperEps < 0 ) {
125 if( ( status = ptwXY_setValueA 125 if( ( status = ptwXY_setValueAtX( ptwXY, x2, 0. ) ) != nfu_Okay ) return( status ); }
126 else { 126 else {
127 if( ( status = ptwXY_setValueA 127 if( ( status = ptwXY_setValueAtX( ptwXY, xp, 0. ) ) != nfu_Okay ) return( status );
128 if( ( status = ptwXY_interpola 128 if( ( status = ptwXY_interpolatePoint( ptwXY->interpolation, x2, &y, xm, y, xp, 0. ) ) != nfu_Okay ) return( status );
129 if( ( status = ptwXY_setValueA 129 if( ( status = ptwXY_setValueAtX( ptwXY, x2, y ) ) != nfu_Okay ) return( status );
130 } 130 }
131 } 131 }
132 } 132 }
133 133
134 return( ptwXY->status ); 134 return( ptwXY->status );
135 135
136 #undef minEps 136 #undef minEps
137 } 137 }
138 /* 138 /*
139 ********************************************** 139 ************************************************************
140 */ 140 */
141 nfu_status ptwXY_mergeClosePoints( ptwXYPoints 141 nfu_status ptwXY_mergeClosePoints( ptwXYPoints *ptwXY, double epsilon ) {
142 142
143 int64_t i, i1, j, k, n = ptwXY->length; 143 int64_t i, i1, j, k, n = ptwXY->length;
144 double x, y; 144 double x, y;
145 ptwXYPoint *p1, *p2; 145 ptwXYPoint *p1, *p2;
146 146
147 if( n < 2 ) return( ptwXY->status ); 147 if( n < 2 ) return( ptwXY->status );
148 if( epsilon < 4 * DBL_EPSILON ) epsilon = 148 if( epsilon < 4 * DBL_EPSILON ) epsilon = 4 * DBL_EPSILON;
149 if( ptwXY_simpleCoalescePoints( ptwXY ) != 149 if( ptwXY_simpleCoalescePoints( ptwXY ) != nfu_Okay ) return( ptwXY->status );
150 150
151 p2 = ptwXY->points; 151 p2 = ptwXY->points;
152 x = p2->x; 152 x = p2->x;
153 for( i1 = 1, p2++; i1 < ( n - 1 ); i1++, p 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 * ( 154 if( ( p2->x - x ) > 0.5 * epsilon * ( std::fabs( x ) + std::fabs( p2->x ) ) ) break;
155 } 155 }
156 if( i1 != 1 ) { 156 if( i1 != 1 ) {
157 for( i = i1, p1 = &(ptwXY->points[1]); 157 for( i = i1, p1 = &(ptwXY->points[1]); i < n; i++, p1++, p2++ ) *p1 = *p2;
158 n = ptwXY->length = ptwXY->length - i1 158 n = ptwXY->length = ptwXY->length - i1 + 1;
159 } 159 }
160 160
161 p1 = &(ptwXY->points[n-1]); 161 p1 = &(ptwXY->points[n-1]);
162 x = p1->x; 162 x = p1->x;
163 for( i1 = n - 2, p1--; i1 > 0; i1--, p1-- 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: 164 if( x - p1->x > 0.5 * epsilon * ( std::fabs( x ) + std::fabs( p1->x ) ) ) break;
165 } 165 }
166 if( i1 != ( n - 2 ) ) { 166 if( i1 != ( n - 2 ) ) {
167 ptwXY->points[i1 + 1] = ptwXY->points[ 167 ptwXY->points[i1 + 1] = ptwXY->points[n - 1];
168 n = ptwXY->length = i1 + 2; 168 n = ptwXY->length = i1 + 2;
169 } 169 }
170 170
171 for( i = 1; i < n - 1; i++ ) { 171 for( i = 1; i < n - 1; i++ ) {
172 p1 = &(ptwXY->points[i]); 172 p1 = &(ptwXY->points[i]);
173 x = p1->x; 173 x = p1->x;
174 y = p1->y; 174 y = p1->y;
175 for( j = i + 1, p2 = &(ptwXY->points[i 175 for( j = i + 1, p2 = &(ptwXY->points[i+1]); j < n - 1; j++, p2++ ) {
176 if( ( p2->x - p1->x ) > 0.5 * epsi 176 if( ( p2->x - p1->x ) > 0.5 * epsilon * ( std::fabs( p2->x ) + std::fabs( p1->x ) ) ) break;
177 x += p2->x; 177 x += p2->x;
178 y += p2->y; 178 y += p2->y;
179 } 179 }
180 if( ( k = ( j - i ) ) > 1 ) { 180 if( ( k = ( j - i ) ) > 1 ) {
181 p1->x = x / k; 181 p1->x = x / k;
182 p1->y = y / k; 182 p1->y = y / k;
183 for( p1 = &(ptwXY->points[i+1]); j 183 for( p1 = &(ptwXY->points[i+1]); j < n; j++, p1++, p2++ ) *p1 = *p2;
184 n -= ( k - 1 ); 184 n -= ( k - 1 );
185 } 185 }
186 } 186 }
187 ptwXY->length = n; 187 ptwXY->length = n;
188 188
189 return( ptwXY->status ); 189 return( ptwXY->status );
190 } 190 }
191 /* 191 /*
192 ********************************************** 192 ************************************************************
193 */ 193 */
194 ptwXYPoints *ptwXY_intersectionWith_ptwX( ptwX 194 ptwXYPoints *ptwXY_intersectionWith_ptwX( ptwXYPoints *ptwXY, ptwXPoints *ptwX, nfu_status *status ) {
195 195
196 int64_t i, i1, i2, lengthX = ptwX_length( 196 int64_t i, i1, i2, lengthX = ptwX_length( ptwX );
197 double x, y, xMin, xMax; 197 double x, y, xMin, xMax;
198 ptwXYPoints *n = NULL; 198 ptwXYPoints *n = NULL;
199 199
200 if( ( *status = ptwXY->status ) != nfu_Oka 200 if( ( *status = ptwXY->status ) != nfu_Okay ) return( NULL );
201 if( ( *status = ptwX->status ) != nfu_Okay 201 if( ( *status = ptwX->status ) != nfu_Okay ) return( NULL );
202 if( ( *status = ptwXY_simpleCoalescePoints 202 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) goto Err;
203 *status = nfu_otherInterpolation; 203 *status = nfu_otherInterpolation;
204 if( ptwXY->interpolation == ptwXY_interpol 204 if( ptwXY->interpolation == ptwXY_interpolationOther ) return( NULL );
205 if( ( n = ptwXY_clone( ptwXY, status ) ) = 205 if( ( n = ptwXY_clone( ptwXY, status ) ) == NULL ) return( NULL );
206 if( ptwXY->length == 0 ) return( n ); 206 if( ptwXY->length == 0 ) return( n );
207 xMin = ptwXY->points[0].x; 207 xMin = ptwXY->points[0].x;
208 xMax = ptwXY->points[ptwXY->length - 1].x; 208 xMax = ptwXY->points[ptwXY->length - 1].x;
209 209
210 if( ( xMin >= ptwX->points[lengthX-1] ) || 210 if( ( xMin >= ptwX->points[lengthX-1] ) || ( xMax <= ptwX->points[0] ) ) { /* No overlap. */
211 n->length = 0; 211 n->length = 0;
212 return( n ); 212 return( n );
213 } 213 }
214 214
215 for( i = 0; i < lengthX; i++ ) { /* 215 for( i = 0; i < lengthX; i++ ) { /* Fill in ptwXY at x-points in ptwX. */
216 x = ptwX->points[i]; 216 x = ptwX->points[i];
217 if( x <= xMin ) continue; 217 if( x <= xMin ) continue;
218 if( x >= xMax ) break; 218 if( x >= xMax ) break;
219 if( ( *status = ptwXY_getValueAtX( ptw 219 if( ( *status = ptwXY_getValueAtX( ptwXY, x, &y ) ) != nfu_Okay ) goto Err;
220 if( ( *status = ptwXY_setValueAtX( n, 220 if( ( *status = ptwXY_setValueAtX( n, x, y ) ) != nfu_Okay ) goto Err;
221 } 221 }
222 if( ( *status = ptwXY_simpleCoalescePoints 222 if( ( *status = ptwXY_simpleCoalescePoints( n ) ) != nfu_Okay ) goto Err;
223 223
224 i1 = 0; 224 i1 = 0;
225 i2 = n->length - 1; 225 i2 = n->length - 1;
226 if( lengthX > 0 ) { 226 if( lengthX > 0 ) {
227 x = ptwX->points[0]; 227 x = ptwX->points[0];
228 if( x > n->points[i1].x ) { 228 if( x > n->points[i1].x ) {
229 for( ; i1 < n->length; i1++ ) { 229 for( ; i1 < n->length; i1++ ) {
230 if( n->points[i1].x == x ) bre 230 if( n->points[i1].x == x ) break;
231 } 231 }
232 } 232 }
233 233
234 x = ptwX->points[lengthX - 1]; 234 x = ptwX->points[lengthX - 1];
235 if( x < n->points[i2].x ) { 235 if( x < n->points[i2].x ) {
236 for( ; i2 > i1; i2-- ) { 236 for( ; i2 > i1; i2-- ) {
237 if( n->points[i2].x == x ) bre 237 if( n->points[i2].x == x ) break;
238 } 238 }
239 } 239 }
240 } 240 }
241 i2++; 241 i2++;
242 242
243 if( i1 != 0 ) { 243 if( i1 != 0 ) {
244 for( i = i1; i < i2; i++ ) n->points[i 244 for( i = i1; i < i2; i++ ) n->points[i - i1] = n->points[i];
245 } 245 }
246 n->length = i2 - i1; 246 n->length = i2 - i1;
247 247
248 return( n ); 248 return( n );
249 249
250 Err: 250 Err:
251 ptwXY_free( n ); 251 ptwXY_free( n );
252 return( NULL ); 252 return( NULL );
253 } 253 }
254 /* 254 /*
255 ********************************************** 255 ************************************************************
256 */ 256 */
257 nfu_status ptwXY_areDomainsMutual( ptwXYPoints 257 nfu_status ptwXY_areDomainsMutual( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2 ) {
258 258
259 nfu_status status; 259 nfu_status status;
260 int64_t n1 = ptwXY1->length, n2 = ptwXY2-> 260 int64_t n1 = ptwXY1->length, n2 = ptwXY2->length;
261 ptwXYPoint *xy1, *xy2; 261 ptwXYPoint *xy1, *xy2;
262 262
263 if( ( status = ptwXY1->status ) != nfu_Oka 263 if( ( status = ptwXY1->status ) != nfu_Okay ) return( status );
264 if( ( status = ptwXY2->status ) != nfu_Oka 264 if( ( status = ptwXY2->status ) != nfu_Okay ) return( status );
265 if( n1 == 0 ) return( nfu_empty ); 265 if( n1 == 0 ) return( nfu_empty );
266 if( n2 == 0 ) return( nfu_empty ); 266 if( n2 == 0 ) return( nfu_empty );
267 if( n1 < 2 ) { 267 if( n1 < 2 ) {
268 status = nfu_tooFewPoints; } 268 status = nfu_tooFewPoints; }
269 else if( n2 < 2 ) { 269 else if( n2 < 2 ) {
270 status = nfu_tooFewPoints; } 270 status = nfu_tooFewPoints; }
271 else { 271 else {
272 xy1 = ptwXY_getPointAtIndex_Unsafely( 272 xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, 0 );
273 xy2 = ptwXY_getPointAtIndex_Unsafely( 273 xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, 0 );
274 if( xy1->x < xy2->x ) { 274 if( xy1->x < xy2->x ) {
275 if( xy2->y != 0. ) status = nfu_do 275 if( xy2->y != 0. ) status = nfu_domainsNotMutual; }
276 else if( xy1->x > xy2->x ) { 276 else if( xy1->x > xy2->x ) {
277 if( xy1->y != 0. ) status = nfu_do 277 if( xy1->y != 0. ) status = nfu_domainsNotMutual;
278 } 278 }
279 279
280 if( status == nfu_Okay ) { 280 if( status == nfu_Okay ) {
281 xy1 = ptwXY_getPointAtIndex_Unsafe 281 xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, n1 - 1 );
282 xy2 = ptwXY_getPointAtIndex_Unsafe 282 xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, n2 - 1 );
283 if( xy1->x < xy2->x ) { 283 if( xy1->x < xy2->x ) {
284 if( xy1->y != 0. ) status = nf 284 if( xy1->y != 0. ) status = nfu_domainsNotMutual; }
285 else if( xy1->x > xy2->x ) { 285 else if( xy1->x > xy2->x ) {
286 if( xy2->y != 0. ) status = nf 286 if( xy2->y != 0. ) status = nfu_domainsNotMutual;
287 } 287 }
288 } 288 }
289 } 289 }
290 return( status ); 290 return( status );
291 } 291 }
292 /* 292 /*
293 ********************************************** 293 ************************************************************
294 */ 294 */
295 nfu_status ptwXY_tweakDomainsToMutualify( ptwX 295 nfu_status ptwXY_tweakDomainsToMutualify( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, int epsilonFactor, double epsilon ) {
296 296
297 nfu_status status; 297 nfu_status status;
298 int64_t n1 = ptwXY1->length, n2 = ptwXY2-> 298 int64_t n1 = ptwXY1->length, n2 = ptwXY2->length;
299 double sum, diff; 299 double sum, diff;
300 ptwXYPoint *xy1, *xy2; 300 ptwXYPoint *xy1, *xy2;
301 301
302 epsilon = std::fabs( epsilon ) + std::fabs 302 epsilon = std::fabs( epsilon ) + std::fabs( epsilonFactor * DBL_EPSILON );
303 303
304 if( ( status = ptwXY1->status ) != nfu_Oka 304 if( ( status = ptwXY1->status ) != nfu_Okay ) return( status );
305 if( ( status = ptwXY2->status ) != nfu_Oka 305 if( ( status = ptwXY2->status ) != nfu_Okay ) return( status );
306 if( n1 == 0 ) return( nfu_empty ); 306 if( n1 == 0 ) return( nfu_empty );
307 if( n2 == 0 ) return( nfu_empty ); 307 if( n2 == 0 ) return( nfu_empty );
308 if( n1 < 2 ) { 308 if( n1 < 2 ) {
309 status = nfu_tooFewPoints; } 309 status = nfu_tooFewPoints; }
310 else if( n2 < 2 ) { 310 else if( n2 < 2 ) {
311 status = nfu_tooFewPoints; } 311 status = nfu_tooFewPoints; }
312 else { 312 else {
313 xy1 = ptwXY_getPointAtIndex_Unsafely( 313 xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, 0 );
314 xy2 = ptwXY_getPointAtIndex_Unsafely( 314 xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, 0 );
315 if( xy1->x < xy2->x ) { 315 if( xy1->x < xy2->x ) {
316 if( xy2->y != 0. ) { 316 if( xy2->y != 0. ) {
317 sum = std::fabs( xy1->x ) + st 317 sum = std::fabs( xy1->x ) + std::fabs( xy2->x );
318 diff = std::fabs( xy2->x - xy1 318 diff = std::fabs( xy2->x - xy1->x );
319 if( diff > epsilon * sum ) { 319 if( diff > epsilon * sum ) {
320 status = nfu_domainsNotMut 320 status = nfu_domainsNotMutual; }
321 else { 321 else {
322 xy1->x = xy2->x; 322 xy1->x = xy2->x;
323 } 323 }
324 } } 324 } }
325 else if( xy1->x > xy2->x ) { 325 else if( xy1->x > xy2->x ) {
326 if( xy1->y != 0. ) { 326 if( xy1->y != 0. ) {
327 sum = std::fabs( xy1->x ) + st 327 sum = std::fabs( xy1->x ) + std::fabs( xy2->x );
328 diff = std::fabs( xy2->x - xy1 328 diff = std::fabs( xy2->x - xy1->x );
329 if( diff > epsilon * sum ) { 329 if( diff > epsilon * sum ) {
330 status = nfu_domainsNotMut 330 status = nfu_domainsNotMutual; }
331 else { 331 else {
332 xy2->x = xy1->x; 332 xy2->x = xy1->x;
333 } 333 }
334 } 334 }
335 } 335 }
336 336
337 if( status == nfu_Okay ) { 337 if( status == nfu_Okay ) {
338 xy1 = ptwXY_getPointAtIndex_Unsafe 338 xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, n1 - 1 );
339 xy2 = ptwXY_getPointAtIndex_Unsafe 339 xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, n2 - 1 );
340 if( xy1->x < xy2->x ) { 340 if( xy1->x < xy2->x ) {
341 if( xy1->y != 0. ) { 341 if( xy1->y != 0. ) {
342 sum = std::fabs( xy1->x ) 342 sum = std::fabs( xy1->x ) + std::fabs( xy2->x );
343 diff = std::fabs( xy2->x - 343 diff = std::fabs( xy2->x - xy1->x );
344 if( diff > epsilon * sum ) 344 if( diff > epsilon * sum ) {
345 status = nfu_domainsNo 345 status = nfu_domainsNotMutual; }
346 else { 346 else {
347 xy2->x = xy1->x; 347 xy2->x = xy1->x;
348 } 348 }
349 } } 349 } }
350 else if( xy1->x > xy2->x ) { 350 else if( xy1->x > xy2->x ) {
351 if( xy2->y != 0. ) { 351 if( xy2->y != 0. ) {
352 sum = std::fabs( xy1->x ) 352 sum = std::fabs( xy1->x ) + std::fabs( xy2->x );
353 diff = std::fabs( xy2->x - 353 diff = std::fabs( xy2->x - xy1->x );
354 if( diff > epsilon * sum ) 354 if( diff > epsilon * sum ) {
355 status = nfu_domainsNo 355 status = nfu_domainsNotMutual; }
356 else { 356 else {
357 xy1->x = xy2->x; 357 xy1->x = xy2->x;
358 } 358 }
359 } 359 }
360 } 360 }
361 } 361 }
362 } 362 }
363 return( status ); 363 return( status );
364 } 364 }
365 /* 365 /*
366 ********************************************** 366 ************************************************************
367 */ 367 */
368 nfu_status ptwXY_mutualifyDomains( ptwXYPoints 368 nfu_status ptwXY_mutualifyDomains( ptwXYPoints *ptwXY1, double lowerEps1, double upperEps1, int positiveXOnly1,
369 ptwX 369 ptwXYPoints *ptwXY2, double lowerEps2, double upperEps2, int positiveXOnly2 ) {
370 370
371 nfu_status status; 371 nfu_status status;
372 int64_t n1 = ptwXY1->length, n2 = ptwXY2-> 372 int64_t n1 = ptwXY1->length, n2 = ptwXY2->length;
373 ptwXYPoint *xy1, *xy2; 373 ptwXYPoint *xy1, *xy2;
374 374
375 switch( status = ptwXY_areDomainsMutual( p 375 switch( status = ptwXY_areDomainsMutual( ptwXY1, ptwXY2 ) ) {
376 case nfu_Okay : 376 case nfu_Okay :
377 case nfu_empty : 377 case nfu_empty :
378 return( nfu_Okay ); 378 return( nfu_Okay );
379 case nfu_domainsNotMutual : 379 case nfu_domainsNotMutual :
380 break; 380 break;
381 default : 381 default :
382 return( status ); 382 return( status );
383 } 383 }
384 384
385 if( ptwXY1->interpolation == ptwXY_interpo 385 if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
386 if( ptwXY2->interpolation == ptwXY_interpo 386 if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
387 if( ptwXY1->interpolation == ptwXY_interpo 387 if( ptwXY1->interpolation == ptwXY_interpolationFlat ) return( nfu_invalidInterpolation );
388 if( ptwXY2->interpolation == ptwXY_interpo 388 if( ptwXY2->interpolation == ptwXY_interpolationFlat ) return( nfu_invalidInterpolation );
389 389
390 xy1 = ptwXY_getPointAtIndex_Unsafely( ptwX 390 xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, 0 );
391 xy2 = ptwXY_getPointAtIndex_Unsafely( ptwX 391 xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, 0 );
392 if( xy1->x < xy2->x ) { 392 if( xy1->x < xy2->x ) {
393 lowerEps1 = 0.; 393 lowerEps1 = 0.;
394 if( xy2->y == 0. ) lowerEps2 = 0.; } 394 if( xy2->y == 0. ) lowerEps2 = 0.; }
395 else if( xy1->x > xy2->x ) { 395 else if( xy1->x > xy2->x ) {
396 lowerEps2 = 0.; 396 lowerEps2 = 0.;
397 if( xy1->y == 0. ) lowerEps1 = 0.; } 397 if( xy1->y == 0. ) lowerEps1 = 0.; }
398 else { 398 else {
399 lowerEps1 = lowerEps2 = 0.; 399 lowerEps1 = lowerEps2 = 0.;
400 } 400 }
401 401
402 xy1 = ptwXY_getPointAtIndex_Unsafely( ptwX 402 xy1 = ptwXY_getPointAtIndex_Unsafely( ptwXY1, n1 - 1 );
403 xy2 = ptwXY_getPointAtIndex_Unsafely( ptwX 403 xy2 = ptwXY_getPointAtIndex_Unsafely( ptwXY2, n2 - 1 );
404 if( xy1->x < xy2->x ) { 404 if( xy1->x < xy2->x ) {
405 upperEps2 = 0.; 405 upperEps2 = 0.;
406 if( xy1->y == 0. ) upperEps1 = 0.; } 406 if( xy1->y == 0. ) upperEps1 = 0.; }
407 else if( xy1->x > xy2->x ) { 407 else if( xy1->x > xy2->x ) {
408 upperEps1 = 0.; 408 upperEps1 = 0.;
409 if( xy2->y == 0. ) upperEps2 = 0.; } 409 if( xy2->y == 0. ) upperEps2 = 0.; }
410 else { 410 else {
411 upperEps1 = upperEps2 = 0.; 411 upperEps1 = upperEps2 = 0.;
412 } 412 }
413 413
414 if( ( lowerEps1 != 0. ) || ( upperEps1 != 414 if( ( lowerEps1 != 0. ) || ( upperEps1 != 0. ) )
415 if( ( status = ptwXY_dullEdges( ptwXY1 415 if( ( status = ptwXY_dullEdges( ptwXY1, lowerEps1, upperEps1, positiveXOnly1 ) ) != nfu_Okay ) return( status );
416 if( ( lowerEps2 != 0. ) || ( upperEps2 != 416 if( ( lowerEps2 != 0. ) || ( upperEps2 != 0. ) )
417 if( ( status = ptwXY_dullEdges( ptwXY2 417 if( ( status = ptwXY_dullEdges( ptwXY2, lowerEps2, upperEps2, positiveXOnly2 ) ) != nfu_Okay ) return( status );
418 418
419 return( status ); 419 return( status );
420 } 420 }
421 /* 421 /*
422 ********************************************** 422 ************************************************************
423 */ 423 */
424 nfu_status ptwXY_copyToC_XY( ptwXYPoints *ptwX 424 nfu_status ptwXY_copyToC_XY( ptwXYPoints *ptwXY, int64_t index1, int64_t index2, int64_t allocatedSize, int64_t *numberOfPoints, double *xys ) {
425 425
426 int64_t i; 426 int64_t i;
427 double *d = xys; 427 double *d = xys;
428 nfu_status status; 428 nfu_status status;
429 ptwXYPoint *pointFrom; 429 ptwXYPoint *pointFrom;
430 430
431 if( ptwXY->status != nfu_Okay ) return( pt 431 if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
432 if( ( status = ptwXY_simpleCoalescePoints( 432 if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status );
433 433
434 if( index1 < 0 ) index1 = 0; 434 if( index1 < 0 ) index1 = 0;
435 if( index2 > ptwXY->length ) index2 = ptwX 435 if( index2 > ptwXY->length ) index2 = ptwXY->length;
436 if( index2 < index1 ) index2 = index1; 436 if( index2 < index1 ) index2 = index1;
437 *numberOfPoints = index2 - index1; 437 *numberOfPoints = index2 - index1;
438 if( allocatedSize < ( index2 - index1 ) ) 438 if( allocatedSize < ( index2 - index1 ) ) return( nfu_insufficientMemory );
439 439
440 for( i = index1, pointFrom = ptwXY->points 440 for( i = index1, pointFrom = ptwXY->points; i < index2; i++, pointFrom++ ) {
441 *(d++) = pointFrom->x; 441 *(d++) = pointFrom->x;
442 *(d++) = pointFrom->y; 442 *(d++) = pointFrom->y;
443 } 443 }
444 444
445 return( nfu_Okay ); 445 return( nfu_Okay );
446 } 446 }
447 /* 447 /*
448 ********************************************** 448 ************************************************************
449 */ 449 */
450 nfu_status ptwXY_valueTo_ptwXAndY( ptwXYPoints 450 nfu_status ptwXY_valueTo_ptwXAndY( ptwXYPoints *ptwXY, double **xs, double **ys ) {
451 451
452 int64_t i1, length = ptwXY_length( ptwXY ) 452 int64_t i1, length = ptwXY_length( ptwXY );
453 double *xps, *yps; 453 double *xps, *yps;
454 ptwXYPoint *pointFrom; 454 ptwXYPoint *pointFrom;
455 nfu_status status; 455 nfu_status status;
456 456
457 if( ptwXY->status != nfu_Okay ) return( pt 457 if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
458 if( ( status = ptwXY_simpleCoalescePoints( 458 if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status );
459 459
460 if( ( *xs = (double *) malloc( length * si 460 if( ( *xs = (double *) malloc( length * sizeof( double ) ) ) == NULL ) return( nfu_mallocError );
461 if( ( *ys = (double *) malloc( length * si 461 if( ( *ys = (double *) malloc( length * sizeof( double ) ) ) == NULL ) {
462 free( *xs ); 462 free( *xs );
463 *xs = NULL; 463 *xs = NULL;
464 return( nfu_mallocError ); 464 return( nfu_mallocError );
465 } 465 }
466 466
467 for( i1 = 0, pointFrom = ptwXY->points, xp 467 for( i1 = 0, pointFrom = ptwXY->points, xps = *xs, yps = *ys; i1 < length; ++i1, ++pointFrom, ++xps, ++yps ) {
468 *xps = pointFrom->x; 468 *xps = pointFrom->x;
469 *yps = pointFrom->y; 469 *yps = pointFrom->y;
470 } 470 }
471 471
472 return( nfu_Okay ); 472 return( nfu_Okay );
473 } 473 }
474 /* 474 /*
475 ********************************************** 475 ************************************************************
476 */ 476 */
477 ptwXYPoints *ptwXY_valueTo_ptwXY( double x1, d 477 ptwXYPoints *ptwXY_valueTo_ptwXY( double x1, double x2, double y, nfu_status *status ) {
478 478
479 ptwXYPoints *n; 479 ptwXYPoints *n;
480 480
481 *status = nfu_XNotAscending; 481 *status = nfu_XNotAscending;
482 if( x1 >= x2 ) return( NULL ); 482 if( x1 >= x2 ) return( NULL );
483 *status = nfu_Okay; 483 *status = nfu_Okay;
484 if( ( n = ptwXY_new( ptwXY_interpolationLi 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 ); 485 ptwXY_setValueAtX( n, x1, y );
486 ptwXY_setValueAtX( n, x2, y ); 486 ptwXY_setValueAtX( n, x2, y );
487 return( n ); 487 return( n );
488 } 488 }
489 /* 489 /*
490 ********************************************** 490 ************************************************************
491 */ 491 */
492 ptwXYPoints *ptwXY_createGaussianCenteredSigma 492 ptwXYPoints *ptwXY_createGaussianCenteredSigma1( double accuracy, nfu_status *status ) {
493 493
494 int64_t i, n; 494 int64_t i, n;
495 ptwXYPoint *pm, *pp; 495 ptwXYPoint *pm, *pp;
496 double x1, y1, x2, y2, accuracy2, yMin = 1 496 double x1, y1, x2, y2, accuracy2, yMin = 1e-10;
497 ptwXYPoints *gaussian; 497 ptwXYPoints *gaussian;
498 498
499 if( accuracy < 1e-5 ) accuracy = 1e-5; 499 if( accuracy < 1e-5 ) accuracy = 1e-5;
500 if( accuracy > 1e-1 ) accuracy = 1e-1; 500 if( accuracy > 1e-1 ) accuracy = 1e-1;
501 if( ( gaussian = ptwXY_new( ptwXY_interpol 501 if( ( gaussian = ptwXY_new( ptwXY_interpolationLinLin, NULL, 1., accuracy, 200, 100, status, 0 ) ) == NULL ) return( NULL );
502 accuracy2 = accuracy = gaussian->accuracy; 502 accuracy2 = accuracy = gaussian->accuracy;
503 if( accuracy2 > 5e-3 ) accuracy2 = 5e-3; 503 if( accuracy2 > 5e-3 ) accuracy2 = 5e-3;
504 504
505 x1 = -std::sqrt( -2. * G4Log( yMin ) ); 505 x1 = -std::sqrt( -2. * G4Log( yMin ) );
506 y1 = yMin; 506 y1 = yMin;
507 x2 = -5.2; 507 x2 = -5.2;
508 y2 = G4Exp( -0.5 * x2 * x2 ); 508 y2 = G4Exp( -0.5 * x2 * x2 );
509 if( ( *status = ptwXY_setValueAtX( gaussia 509 if( ( *status = ptwXY_setValueAtX( gaussian, x1, y1 ) ) != nfu_Okay ) goto Err;
510 gaussian->accuracy = 20 * accuracy2; 510 gaussian->accuracy = 20 * accuracy2;
511 if( ( *status = ptwXY_createGaussianCenter 511 if( ( *status = ptwXY_createGaussianCenteredSigma1_2( gaussian, x1, y1, x2, y2, 1 ) ) != nfu_Okay ) goto Err;
512 x1 = x2; 512 x1 = x2;
513 y1 = y2; 513 y1 = y2;
514 x2 = -4.; 514 x2 = -4.;
515 y2 = G4Exp( -0.5 * x2 * x2 ); 515 y2 = G4Exp( -0.5 * x2 * x2 );
516 gaussian->accuracy = 5 * accuracy2; 516 gaussian->accuracy = 5 * accuracy2;
517 if( ( *status = ptwXY_createGaussianCenter 517 if( ( *status = ptwXY_createGaussianCenteredSigma1_2( gaussian, x1, y1, x2, y2, 1 ) ) != nfu_Okay ) goto Err;
518 x1 = x2; 518 x1 = x2;
519 y1 = y2; 519 y1 = y2;
520 x2 = -1; 520 x2 = -1;
521 y2 = G4Exp( -0.5 * x2 * x2 ); 521 y2 = G4Exp( -0.5 * x2 * x2 );
522 gaussian->accuracy = accuracy; 522 gaussian->accuracy = accuracy;
523 if( ( *status = ptwXY_createGaussianCenter 523 if( ( *status = ptwXY_createGaussianCenteredSigma1_2( gaussian, x1, y1, x2, y2, 1 ) ) != nfu_Okay ) goto Err;
524 x1 = x2; 524 x1 = x2;
525 y1 = y2; 525 y1 = y2;
526 x2 = 0; 526 x2 = 0;
527 y2 = G4Exp( -0.5 * x2 * x2 ); 527 y2 = G4Exp( -0.5 * x2 * x2 );
528 if( ( *status = ptwXY_createGaussianCenter 528 if( ( *status = ptwXY_createGaussianCenteredSigma1_2( gaussian, x1, y1, x2, y2, 1 ) ) != nfu_Okay ) goto Err;
529 529
530 n = gaussian->length; 530 n = gaussian->length;
531 if( ( *status = ptwXY_coalescePoints( gaus 531 if( ( *status = ptwXY_coalescePoints( gaussian, 2 * n + 1, NULL, 0 ) ) != nfu_Okay ) goto Err;
532 if( ( *status = ptwXY_setValueAtX( gaussia 532 if( ( *status = ptwXY_setValueAtX( gaussian, 0., 1. ) ) != nfu_Okay ) goto Err;
533 pp = &(gaussian->points[gaussian->length]) 533 pp = &(gaussian->points[gaussian->length]);
534 for( i = 0, pm = pp - 2; i < n; i++, pp++, 534 for( i = 0, pm = pp - 2; i < n; i++, pp++, pm-- ) {
535 *pp = *pm; 535 *pp = *pm;
536 pp->x *= -1; 536 pp->x *= -1;
537 } 537 }
538 gaussian->length = 2 * n + 1; 538 gaussian->length = 2 * n + 1;
539 539
540 return( gaussian ); 540 return( gaussian );
541 541
542 Err: 542 Err:
543 ptwXY_free( gaussian ); 543 ptwXY_free( gaussian );
544 return( NULL ); 544 return( NULL );
545 } 545 }
546 /* 546 /*
547 ********************************************** 547 ************************************************************
548 */ 548 */
549 static nfu_status ptwXY_createGaussianCentered 549 static nfu_status ptwXY_createGaussianCenteredSigma1_2( ptwXYPoints *ptwXY, double x1, double y1, double x2, double y2, int addX1Point ) {
550 550
551 nfu_status status = nfu_Okay; 551 nfu_status status = nfu_Okay;
552 int morePoints = 0; 552 int morePoints = 0;
553 double x = 0.5 * ( x1 + x2 ); 553 double x = 0.5 * ( x1 + x2 );
554 double y = G4Exp( -x * x / 2 ), yMin = ( y 554 double y = G4Exp( -x * x / 2 ), yMin = ( y1 * ( x2 - x ) + y2 * ( x - x1 ) ) / ( x2 - x1 );
555 555
556 if( std::fabs( y - yMin ) > y * ptwXY->acc 556 if( std::fabs( y - yMin ) > y * ptwXY->accuracy ) morePoints = 1;
557 if( morePoints && ( status = ptwXY_createG 557 if( morePoints && ( status = ptwXY_createGaussianCenteredSigma1_2( ptwXY, x, y, x2, y2, 0 ) ) != nfu_Okay ) return( status );
558 if( ( status = ptwXY_setValueAtX( ptwXY, x 558 if( ( status = ptwXY_setValueAtX( ptwXY, x, y ) ) != nfu_Okay ) return( status );
559 if( morePoints && ( status = ptwXY_createG 559 if( morePoints && ( status = ptwXY_createGaussianCenteredSigma1_2( ptwXY, x1, y1, x, y, 0 ) ) != nfu_Okay ) return( status );
560 if( addX1Point ) status = ptwXY_setValueAt 560 if( addX1Point ) status = ptwXY_setValueAtX( ptwXY, x1, y1 );
561 return( status ); 561 return( status );
562 } 562 }
563 /* 563 /*
564 ********************************************** 564 ************************************************************
565 */ 565 */
566 ptwXYPoints *ptwXY_createGaussian( double accu 566 ptwXYPoints *ptwXY_createGaussian( double accuracy, double xCenter, double sigma, double amplitude, double xMin, double xMax,
567 double /*dullEps*/, nfu_status *status 567 double /*dullEps*/, nfu_status *status ) {
568 568
569 int64_t i; 569 int64_t i;
570 ptwXYPoints *gaussian, *sliced; 570 ptwXYPoints *gaussian, *sliced;
571 ptwXYPoint *point; 571 ptwXYPoint *point;
572 572
573 if( ( gaussian = ptwXY_createGaussianCente 573 if( ( gaussian = ptwXY_createGaussianCenteredSigma1( accuracy, status ) ) == NULL ) return( NULL );
574 for( i = 0, point = gaussian->points; i < 574 for( i = 0, point = gaussian->points; i < gaussian->length; i++, point++ ) {
575 point->x = point->x * sigma + xCenter; 575 point->x = point->x * sigma + xCenter;
576 point->y *= amplitude; 576 point->y *= amplitude;
577 } 577 }
578 if( ( gaussian->points[0].x < xMin ) || ( 578 if( ( gaussian->points[0].x < xMin ) || ( gaussian->points[gaussian->length - 1].x > xMax ) ) {
579 if( ( sliced = ptwXY_xSlice( gaussian, 579 if( ( sliced = ptwXY_xSlice( gaussian, xMin, xMax, 10, 1, status ) ) == NULL ) goto Err;
580 ptwXY_free( gaussian ); 580 ptwXY_free( gaussian );
581 gaussian = sliced; 581 gaussian = sliced;
582 } 582 }
583 583
584 return( gaussian ); 584 return( gaussian );
585 585
586 Err: 586 Err:
587 ptwXY_free( gaussian ); 587 ptwXY_free( gaussian );
588 return( NULL ); 588 return( NULL );
589 } 589 }
590 590
591 #if defined __cplusplus 591 #if defined __cplusplus
592 } 592 }
593 #endif 593 #endif
594 594