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