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 7 8 #include "ptwXY.h" 8 #include "ptwXY.h" 9 9 10 #if defined __cplusplus 10 #if defined __cplusplus 11 #include "G4Exp.hh" 11 #include "G4Exp.hh" 12 #include "G4Pow.hh" 12 #include "G4Pow.hh" 13 namespace GIDI { 13 namespace GIDI { 14 using namespace GIDI; 14 using namespace GIDI; 15 #endif 15 #endif 16 16 17 static nfu_status ptwXY_pow_callback( ptwXYPoi 17 static nfu_status ptwXY_pow_callback( ptwXYPoint *point, void *argList ); 18 static nfu_status ptwXY_exp_s( ptwXYPoints *pt 18 static nfu_status ptwXY_exp_s( ptwXYPoints *ptwXY, double x1, double y1, double z1, double x2, double y2, double z2, int level ); 19 static nfu_status ptwXY_convolution2( ptwXYPoi 19 static nfu_status ptwXY_convolution2( ptwXYPoints *f1, ptwXYPoints *f2, double y, double yMin, double *c ); 20 static nfu_status ptwXY_convolution3( ptwXYPoi 20 static nfu_status ptwXY_convolution3( ptwXYPoints *convolute, ptwXYPoints *f1, ptwXYPoints *f2, double y1, double c1, double y2, double c2, double yMin ); 21 /* 21 /* 22 ********************************************** 22 ************************************************************ 23 */ 23 */ 24 nfu_status ptwXY_pow( ptwXYPoints *ptwXY, doub 24 nfu_status ptwXY_pow( ptwXYPoints *ptwXY, double v ) { 25 25 26 return( ptwXY_applyFunction( ptwXY, ptwXY_ 26 return( ptwXY_applyFunction( ptwXY, ptwXY_pow_callback, (void *) &v, 0 ) ); 27 } 27 } 28 /* 28 /* 29 ********************************************** 29 ************************************************************ 30 */ 30 */ 31 static nfu_status ptwXY_pow_callback( ptwXYPoi 31 static nfu_status ptwXY_pow_callback( ptwXYPoint *point, void *argList ) { 32 32 33 nfu_status status = nfu_Okay; 33 nfu_status status = nfu_Okay; 34 double *v = (double *) argList; 34 double *v = (double *) argList; 35 35 36 point->y = G4Pow::GetInstance()->powA( poi 36 point->y = G4Pow::GetInstance()->powA( point->y, *v ); 37 /* ???? Need to test for valid y-value. */ 37 /* ???? Need to test for valid y-value. */ 38 return( status ); 38 return( status ); 39 } 39 } 40 /* 40 /* 41 ********************************************** 41 ************************************************************ 42 */ 42 */ 43 nfu_status ptwXY_exp( ptwXYPoints *ptwXY, doub 43 nfu_status ptwXY_exp( ptwXYPoints *ptwXY, double a ) { 44 44 45 int64_t i, length; 45 int64_t i, length; 46 nfu_status status; 46 nfu_status status; 47 double x1, y1, z1, x2, y2, z2; 47 double x1, y1, z1, x2, y2, z2; 48 48 49 length = ptwXY->length; 49 length = ptwXY->length; 50 if( length < 1 ) return( ptwXY->status ); 50 if( length < 1 ) return( ptwXY->status ); 51 if( ptwXY->interpolation == ptwXY_interpol 51 if( ptwXY->interpolation == ptwXY_interpolationFlat ) return( nfu_invalidInterpolation ); 52 if( ptwXY->interpolation == ptwXY_interpol 52 if( ptwXY->interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation ); 53 53 54 if( ( status = ptwXY_simpleCoalescePoints( 54 if( ( status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( status ); 55 x2 = ptwXY->points[length-1].x; 55 x2 = ptwXY->points[length-1].x; 56 y2 = a * ptwXY->points[length-1].y; 56 y2 = a * ptwXY->points[length-1].y; 57 z2 = ptwXY->points[length-1].y = G4Exp( y2 57 z2 = ptwXY->points[length-1].y = G4Exp( y2 ); 58 for( i = length - 2; i >= 0; i-- ) { 58 for( i = length - 2; i >= 0; i-- ) { 59 x1 = ptwXY->points[i].x; 59 x1 = ptwXY->points[i].x; 60 y1 = a * ptwXY->points[i].y; 60 y1 = a * ptwXY->points[i].y; 61 z1 = ptwXY->points[i].y = G4Exp( y1 ); 61 z1 = ptwXY->points[i].y = G4Exp( y1 ); 62 if( ( status = ptwXY_exp_s( ptwXY, x1, 62 if( ( status = ptwXY_exp_s( ptwXY, x1, y1, z1, x2, y2, z2, 0 ) ) != nfu_Okay ) return( status ); 63 x2 = x1; 63 x2 = x1; 64 y2 = y1; 64 y2 = y1; 65 } 65 } 66 return( status ); 66 return( status ); 67 } 67 } 68 /* 68 /* 69 ********************************************** 69 ************************************************************ 70 */ 70 */ 71 static nfu_status ptwXY_exp_s( ptwXYPoints *pt 71 static nfu_status ptwXY_exp_s( ptwXYPoints *ptwXY, double x1, double y1, double z1, double x2, double y2, double z2, int level ) { 72 72 73 nfu_status status; 73 nfu_status status; 74 double x, y, dx, dy, z, zp, s; 74 double x, y, dx, dy, z, zp, s; 75 75 76 if( ( x1 == x2 ) || ( y1 == y2 ) ) return( 76 if( ( x1 == x2 ) || ( y1 == y2 ) ) return( nfu_Okay ); 77 if( level >= ptwXY->biSectionMax ) return( 77 if( level >= ptwXY->biSectionMax ) return( nfu_Okay ); 78 level++; 78 level++; 79 dx = x2 - x1; 79 dx = x2 - x1; 80 dy = y2 - y1; 80 dy = y2 - y1; 81 s = dy / dx; 81 s = dy / dx; 82 x = 1. / s + x2 - z2 * dx / ( z2 - z1 ); 82 x = 1. / s + x2 - z2 * dx / ( z2 - z1 ); 83 y = ( y1 * ( x2 - x ) + y2 * ( x - x1 ) ) 83 y = ( y1 * ( x2 - x ) + y2 * ( x - x1 ) ) / dx; 84 z = z1 * G4Exp( 1 - dy / ( G4Exp( dy ) - 1 84 z = z1 * G4Exp( 1 - dy / ( G4Exp( dy ) - 1 ) ); 85 zp = ( z2 - z1 ) / ( y2 - y1 ); 85 zp = ( z2 - z1 ) / ( y2 - y1 ); 86 86 87 if( std::fabs( z - zp ) < std::fabs( z * p 87 if( std::fabs( z - zp ) < std::fabs( z * ptwXY->accuracy ) ) return( nfu_Okay ); 88 if( ( status = ptwXY_setValueAtX( ptwXY, x 88 if( ( status = ptwXY_setValueAtX( ptwXY, x, z ) ) != nfu_Okay ) return( status ); 89 if( ( status = ptwXY_exp_s( ptwXY, x, y, z 89 if( ( status = ptwXY_exp_s( ptwXY, x, y, z, x2, y2, z2, level ) ) != nfu_Okay ) return( status ); 90 if( ( status = ptwXY_exp_s( ptwXY, x1, y1, 90 if( ( status = ptwXY_exp_s( ptwXY, x1, y1, z1, x, y, z, level ) ) != nfu_Okay ) return( status ); 91 return( status ); 91 return( status ); 92 } 92 } 93 /* 93 /* 94 ********************************************** 94 ************************************************************ 95 */ 95 */ 96 ptwXYPoints *ptwXY_convolution( ptwXYPoints *p 96 ptwXYPoints *ptwXY_convolution( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int mode ) { 97 /* 97 /* 98 * Currently, only supports linear-linear int 98 * Currently, only supports linear-linear interpolation. 99 * 99 * 100 * This function calculates c(y) = int 100 * This function calculates c(y) = integral dx f1(x) * f2(y-x) 101 * 101 * 102 */ 102 */ 103 int64_t i1, i2, n1, n2, n; 103 int64_t i1, i2, n1, n2, n; 104 ptwXYPoints *f1 = ptwXY1, *f2 = ptwXY2, *c 104 ptwXYPoints *f1 = ptwXY1, *f2 = ptwXY2, *convolute; 105 double accuracy = ptwXY1->accuracy, yMin, 105 double accuracy = ptwXY1->accuracy, yMin, yMax, c, y, dy; 106 106 107 if( ( *status = ptwXY_simpleCoalescePoints 107 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL ); 108 if( ( *status = ptwXY_simpleCoalescePoints 108 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY2 ) ) != nfu_Okay ) return( NULL ); 109 109 110 *status = nfu_unsupportedInterpolation; 110 *status = nfu_unsupportedInterpolation; 111 if( ( ptwXY1->interpolation != ptwXY_inter 111 if( ( ptwXY1->interpolation != ptwXY_interpolationLinLin ) || ( ptwXY2->interpolation != ptwXY_interpolationLinLin ) ) return( NULL ); 112 *status = nfu_Okay; 112 *status = nfu_Okay; 113 113 114 n1 = f1->length; 114 n1 = f1->length; 115 n2 = f2->length; 115 n2 = f2->length; 116 116 117 if( ( n1 == 0 ) || ( n2 == 0 ) ) { 117 if( ( n1 == 0 ) || ( n2 == 0 ) ) { 118 convolute = ptwXY_new( ptwXY_interpola 118 convolute = ptwXY_new( ptwXY_interpolationLinLin, NULL, 1., accuracy, 0, 0, status, 0 ); 119 return( convolute ); 119 return( convolute ); 120 } 120 } 121 121 122 if( ( n1 == 1 ) || ( n2 == 1 ) ) { 122 if( ( n1 == 1 ) || ( n2 == 1 ) ) { 123 *status = nfu_tooFewPoints; 123 *status = nfu_tooFewPoints; 124 return( NULL ); 124 return( NULL ); 125 } 125 } 126 126 127 if( accuracy < ptwXY2->accuracy ) accuracy 127 if( accuracy < ptwXY2->accuracy ) accuracy = ptwXY2->accuracy; 128 n = n1 * n2; 128 n = n1 * n2; 129 if( mode == 0 ) { 129 if( mode == 0 ) { 130 mode = 1; 130 mode = 1; 131 if( n > 1000 ) mode = -1; 131 if( n > 1000 ) mode = -1; 132 } 132 } 133 if( n > 100000 ) mode = -1; 133 if( n > 100000 ) mode = -1; 134 if( ( convolute = ptwXY_new( ptwXY_interpo 134 if( ( convolute = ptwXY_new( ptwXY_interpolationLinLin, NULL, 1., accuracy, 400, 40, status, 0 ) ) == NULL ) return( NULL ); 135 135 136 yMin = f1->points[0].x + f2->points[0].x; 136 yMin = f1->points[0].x + f2->points[0].x; 137 yMax = f1->points[n1 - 1].x + f2->points[n 137 yMax = f1->points[n1 - 1].x + f2->points[n2 - 1].x; 138 138 139 if( ( *status = ptwXY_setValueAtX( convolu 139 if( ( *status = ptwXY_setValueAtX( convolute, yMin, 0. ) ) != nfu_Okay ) goto Err; 140 140 141 if( mode < 0 ) { 141 if( mode < 0 ) { 142 dy = ( yMax - yMin ) / 400; 142 dy = ( yMax - yMin ) / 400; 143 for( y = yMin + dy; y < yMax; y += dy 143 for( y = yMin + dy; y < yMax; y += dy ) { 144 if( ( *status = ptwXY_convolution2 144 if( ( *status = ptwXY_convolution2( f1, f2, y, yMin, &c ) ) != nfu_Okay ) goto Err; 145 if( ( *status = ptwXY_setValueAtX( 145 if( ( *status = ptwXY_setValueAtX( convolute, y, c ) ) != nfu_Okay ) goto Err; 146 } } 146 } } 147 else { 147 else { 148 for( i1 = 0; i1 < n1; i1++ ) { 148 for( i1 = 0; i1 < n1; i1++ ) { 149 for( i2 = 0; i2 < n2; i2++ ) { 149 for( i2 = 0; i2 < n2; i2++ ) { 150 y = yMin + ( f1->points[i1].x 150 y = yMin + ( f1->points[i1].x - f1->points[0].x ) + ( f2->points[i2].x - f2->points[0].x ); 151 if( y <= yMin ) continue; 151 if( y <= yMin ) continue; 152 if( y >= yMax ) continue; 152 if( y >= yMax ) continue; 153 if( ( *status = ptwXY_convolut 153 if( ( *status = ptwXY_convolution2( f1, f2, y, yMin, &c ) ) != nfu_Okay ) goto Err; 154 if( ( *status = ptwXY_setValue 154 if( ( *status = ptwXY_setValueAtX( convolute, y, c ) ) != nfu_Okay ) goto Err; 155 } 155 } 156 } 156 } 157 } 157 } 158 if( ( *status = ptwXY_setValueAtX( convolu 158 if( ( *status = ptwXY_setValueAtX( convolute, yMax, 0. ) ) != nfu_Okay ) goto Err; 159 if( ( *status = ptwXY_simpleCoalescePoints 159 if( ( *status = ptwXY_simpleCoalescePoints( convolute ) ) != nfu_Okay ) goto Err; 160 for( i1 = convolute->length - 1; i1 > 0; i 160 for( i1 = convolute->length - 1; i1 > 0; i1-- ) { 161 if( ( *status = ptwXY_convolution3( co 161 if( ( *status = ptwXY_convolution3( convolute, f1, f2, convolute->points[i1 - 1].x, convolute->points[i1 - 1].y, 162 convolute->points[i1].x, convolute 162 convolute->points[i1].x, convolute->points[i1].y, yMin ) ) != nfu_Okay ) goto Err; 163 } 163 } 164 164 165 return( convolute ); 165 return( convolute ); 166 166 167 Err: 167 Err: 168 ptwXY_free( convolute ); 168 ptwXY_free( convolute ); 169 return( NULL ); 169 return( NULL ); 170 } 170 } 171 /* 171 /* 172 ********************************************** 172 ************************************************************ 173 */ 173 */ 174 static nfu_status ptwXY_convolution2( ptwXYPoi 174 static nfu_status ptwXY_convolution2( ptwXYPoints *f1, ptwXYPoints *f2, double y, double yMin, double *c ) { 175 175 176 int64_t i1 = 0, i2 = 0, n1 = f1->length, n 176 int64_t i1 = 0, i2 = 0, n1 = f1->length, n2 = f2->length, mode; 177 double dx1, dx2, x1MinP, x1Min, x2Max; 177 double dx1, dx2, x1MinP, x1Min, x2Max; 178 double f1x1 = 0, f1y1 = 0, f1x2 = 0, f1 178 double f1x1 = 0, f1y1 = 0, f1x2 = 0, f1y2 = 0, f2x1 = 0, f2y1 = 0, f2x2 = 0, f2y2 = 0; 179 double f1x1p, f1y1p, f1x2p, f1y2p, f2x1p, 179 double f1x1p, f1y1p, f1x2p, f1y2p, f2x1p, f2y1p, f2x2p, f2y2p; 180 ptwXY_lessEqualGreaterX legx; 180 ptwXY_lessEqualGreaterX legx; 181 ptwXYOverflowPoint lessThanEqualXPoint, gr 181 ptwXYOverflowPoint lessThanEqualXPoint, greaterThanXPoint; 182 nfu_status status; 182 nfu_status status; 183 183 184 x2Max = f2->points[0].x + ( y - yMin ); 184 x2Max = f2->points[0].x + ( y - yMin ); 185 if( x2Max > f2->points[n2 - 1].x ) x2Max = 185 if( x2Max > f2->points[n2 - 1].x ) x2Max = f2->points[n2 - 1].x; 186 x1Min = f1->points[0].x; 186 x1Min = f1->points[0].x; 187 x1MinP = y - f2->points[n2 - 1].x; 187 x1MinP = y - f2->points[n2 - 1].x; 188 if( x1Min < x1MinP ) x1Min = x1MinP; 188 if( x1Min < x1MinP ) x1Min = x1MinP; 189 *c = 0.; 189 *c = 0.; 190 190 191 switch( legx = ptwXY_getPointsAroundX( f1, 191 switch( legx = ptwXY_getPointsAroundX( f1, x1Min, &lessThanEqualXPoint, &greaterThanXPoint ) ) { 192 case ptwXY_lessEqualGreaterX_empty : 192 case ptwXY_lessEqualGreaterX_empty : /* These three should not happen. */ 193 case ptwXY_lessEqualGreaterX_lessThan : 193 case ptwXY_lessEqualGreaterX_lessThan : 194 case ptwXY_lessEqualGreaterX_greater : 194 case ptwXY_lessEqualGreaterX_greater : 195 return( nfu_Okay ); 195 return( nfu_Okay ); 196 case ptwXY_lessEqualGreaterX_equal : 196 case ptwXY_lessEqualGreaterX_equal : 197 case ptwXY_lessEqualGreaterX_between : 197 case ptwXY_lessEqualGreaterX_between : 198 i1 = lessThanEqualXPoint.index; 198 i1 = lessThanEqualXPoint.index; 199 f1x1 = f1->points[i1].x; 199 f1x1 = f1->points[i1].x; 200 f1y1p = f1y1 = f1->points[i1].y; 200 f1y1p = f1y1 = f1->points[i1].y; 201 i1++; 201 i1++; 202 if( i1 == n1 ) return( nfu_Okay ); 202 if( i1 == n1 ) return( nfu_Okay ); 203 f1x2 = f1->points[i1].x; 203 f1x2 = f1->points[i1].x; 204 f1y2 = f1->points[i1].y; 204 f1y2 = f1->points[i1].y; 205 if( legx == ptwXY_lessEqualGreaterX_be 205 if( legx == ptwXY_lessEqualGreaterX_between ) { 206 if( ( status = ptwXY_interpolatePo 206 if( ( status = ptwXY_interpolatePoint( f1->interpolation, x1Min, &f1y1p, f1x1, f1y1, f1x2, f1y2 ) ) != nfu_Okay ) return( status ); 207 } 207 } 208 break; 208 break; 209 } 209 } 210 210 211 switch( legx = ptwXY_getPointsAroundX( f2, 211 switch( legx = ptwXY_getPointsAroundX( f2, x2Max, &lessThanEqualXPoint, &greaterThanXPoint ) ) { 212 case ptwXY_lessEqualGreaterX_empty : 212 case ptwXY_lessEqualGreaterX_empty : /* These three should not happen. */ 213 case ptwXY_lessEqualGreaterX_lessThan : 213 case ptwXY_lessEqualGreaterX_lessThan : 214 case ptwXY_lessEqualGreaterX_greater : 214 case ptwXY_lessEqualGreaterX_greater : 215 return( nfu_Okay ); 215 return( nfu_Okay ); 216 case ptwXY_lessEqualGreaterX_equal : 216 case ptwXY_lessEqualGreaterX_equal : 217 case ptwXY_lessEqualGreaterX_between : 217 case ptwXY_lessEqualGreaterX_between : 218 i2 = lessThanEqualXPoint.index; 218 i2 = lessThanEqualXPoint.index; 219 if( i2 < f2->length - 1 ) i2++; 219 if( i2 < f2->length - 1 ) i2++; 220 f2x2 = f2->points[i2].x; 220 f2x2 = f2->points[i2].x; 221 f2y2p = f2y2 = f2->points[i2].y; 221 f2y2p = f2y2 = f2->points[i2].y; 222 i2--; 222 i2--; 223 f2x1 = f2->points[i2].x; 223 f2x1 = f2->points[i2].x; 224 f2y1 = f2->points[i2].y; 224 f2y1 = f2->points[i2].y; 225 if( legx == ptwXY_lessEqualGreaterX_be 225 if( legx == ptwXY_lessEqualGreaterX_between ) { 226 if( ( status = ptwXY_interpolatePo 226 if( ( status = ptwXY_interpolatePoint( f2->interpolation, x2Max, &f2y2p, f2x1, f2y1, f2x2, f2y2 ) ) != nfu_Okay ) return( status ); 227 } 227 } 228 break; 228 break; 229 } 229 } 230 230 231 f1x1p = x1Min; 231 f1x1p = x1Min; 232 f2x2p = x2Max; 232 f2x2p = x2Max; 233 f1y2p = f1y2; 233 f1y2p = f1y2; 234 f2y1p = f2y1; 234 f2y1p = f2y1; 235 while( ( i1 < n1 ) && ( i2 >= 0 ) ) { // L 235 while( ( i1 < n1 ) && ( i2 >= 0 ) ) { // Loop checking, 11.06.2015, T. Koi 236 dx1 = f1x2 - f1x1p; 236 dx1 = f1x2 - f1x1p; 237 dx2 = f2x2p - f2x1; 237 dx2 = f2x2p - f2x1; 238 mode = 2; 238 mode = 2; 239 if( i1 < n1 ) { 239 if( i1 < n1 ) { 240 if( dx1 < dx2 ) mode = 1; 240 if( dx1 < dx2 ) mode = 1; 241 } 241 } 242 if( mode == 1 ) { 242 if( mode == 1 ) { /* Point in f1 is limiting dx step size. */ 243 f2x1p = f2x2p - dx1; 243 f2x1p = f2x2p - dx1; 244 if( f2x1p < f2->points[i2].x ) { 244 if( f2x1p < f2->points[i2].x ) { /* Round off issue may cause this. */ 245 f2x1p = f2x2; 245 f2x1p = f2x2; 246 f2y1p = f2y2; } 246 f2y1p = f2y2; } 247 else { 247 else { 248 if( ( status = ptwXY_interpola 248 if( ( status = ptwXY_interpolatePoint( f2->interpolation, f2x1p, &f2y1p, f2x1, f2y1, f2x2, f2y2 ) ) != nfu_Okay ) return( status ); 249 } 249 } 250 *c += ( ( f1y1p + f1y2p ) * ( f2y1 250 *c += ( ( f1y1p + f1y2p ) * ( f2y1p + f2y2p ) + f1y1p * f2y2p + f1y2p * f2y1p ) * dx1; /* Note the reversing of f2y1p and f2y2p. */ 251 i1++; 251 i1++; 252 if( i1 == n1 ) break; 252 if( i1 == n1 ) break; 253 f1x1p = f1x1 = f1x2; 253 f1x1p = f1x1 = f1x2; 254 f1y1p = f1y1 = f1y2; 254 f1y1p = f1y1 = f1y2; 255 f1x2 = f1->points[i1].x; 255 f1x2 = f1->points[i1].x; 256 f1y2p = f1y2 = f1->points[i1].y; 256 f1y2p = f1y2 = f1->points[i1].y; 257 f2x2p = f2x1p; 257 f2x2p = f2x1p; 258 f2y2p = f2y1p; 258 f2y2p = f2y1p; 259 f2y1p = f2y1; } 259 f2y1p = f2y1; } 260 else { 260 else { 261 f1x2p = f1x1p + dx2; 261 f1x2p = f1x1p + dx2; 262 if( ( f1x2p > f1->points[i1].x ) | 262 if( ( f1x2p > f1->points[i1].x ) || ( dx1 == dx2 ) ) { /* Round off issue may cause first test to trip. */ 263 f1x2p = f1x2; 263 f1x2p = f1x2; 264 f1y2p = f1y2; } 264 f1y2p = f1y2; } 265 else { 265 else { 266 if( ( status = ptwXY_interpola 266 if( ( status = ptwXY_interpolatePoint( f1->interpolation, f1x2p, &f1y2p, f1x1, f1y1, f1x2, f1y2 ) ) != nfu_Okay ) return( status ); 267 } 267 } 268 *c += ( ( f1y1p + f1y2p ) * ( f2y1 268 *c += ( ( f1y1p + f1y2p ) * ( f2y1p + f2y2p ) + f1y1p * f2y2p + f1y2p * f2y1p ) * dx2; /* Note the reversing of f2y1p and f2y2p. */ 269 if( i2 == 0 ) break; 269 if( i2 == 0 ) break; 270 i2--; 270 i2--; 271 f2x2p = f2x2 = f2x1; 271 f2x2p = f2x2 = f2x1; 272 f2y2p = f2y2 = f2y1; 272 f2y2p = f2y2 = f2y1; 273 f2x1 = f2->points[i2].x; 273 f2x1 = f2->points[i2].x; 274 f2y1p = f2y1 = f2->points[i2].y; 274 f2y1p = f2y1 = f2->points[i2].y; 275 f1x1p = f1x2p; 275 f1x1p = f1x2p; 276 if( dx1 == dx2 ) { 276 if( dx1 == dx2 ) { 277 f1x1p = f1x1 = f1x2; 277 f1x1p = f1x1 = f1x2; 278 f1y1p = f1y1 = f1y2; 278 f1y1p = f1y1 = f1y2; 279 i1++; 279 i1++; 280 f1x2 = f1->points[i1].x; 280 f1x2 = f1->points[i1].x; 281 f1y2p = f1y2 = f1->points[i1]. 281 f1y2p = f1y2 = f1->points[i1].y; } 282 else { 282 else { 283 f1y1p = f1y2p; 283 f1y1p = f1y2p; 284 f1y2p = f1->points[i1].y; 284 f1y2p = f1->points[i1].y; 285 } 285 } 286 } 286 } 287 } 287 } 288 *c /= 6.; 288 *c /= 6.; 289 return( nfu_Okay ); 289 return( nfu_Okay ); 290 } 290 } 291 /* 291 /* 292 ********************************************** 292 ************************************************************ 293 */ 293 */ 294 static nfu_status ptwXY_convolution3( ptwXYPoi 294 static nfu_status ptwXY_convolution3( ptwXYPoints *convolute, ptwXYPoints *f1, ptwXYPoints *f2, double y1, double c1, double y2, double c2, double yMin ) { 295 295 296 nfu_status status; 296 nfu_status status; 297 double yMid = 0.5 * ( y1 + y2 ), cMid = 0. 297 double yMid = 0.5 * ( y1 + y2 ), cMid = 0.5 * ( c1 + c2 ), c; 298 298 299 if( ( y2 - yMid ) <= 1e-5 * ( ptwXY_getXMa 299 if( ( y2 - yMid ) <= 1e-5 * ( ptwXY_getXMax( convolute ) - ptwXY_getXMin( convolute ) ) ) return( nfu_Okay ); 300 if( ( status = ptwXY_convolution2( f1, f2, 300 if( ( status = ptwXY_convolution2( f1, f2, yMid, yMin, &c ) ) != nfu_Okay ) return( status ); 301 if( std::fabs( c - cMid ) <= convolute->ac 301 if( std::fabs( c - cMid ) <= convolute->accuracy * 0.5 * ( std::fabs( c ) + std::fabs( cMid ) ) ) return( nfu_Okay ); 302 if( ( status = ptwXY_setValueAtX( convolut 302 if( ( status = ptwXY_setValueAtX( convolute, yMid, c ) ) != nfu_Okay ) return( status ); 303 if( ( status = ptwXY_convolution3( convolu 303 if( ( status = ptwXY_convolution3( convolute, f1, f2, y1, c1, yMid, c, yMin ) ) != nfu_Okay ) return( status ); 304 return( ptwXY_convolution3( convolute, f1, 304 return( ptwXY_convolution3( convolute, f1, f2, yMid, c, y2, c2, yMin ) ); 305 } 305 } 306 306 307 #if defined __cplusplus 307 #if defined __cplusplus 308 } 308 } 309 #endif 309 #endif 310 310