10 #if defined __cplusplus
34 double *
v = (
double *) argList;
50 if( length < 1 )
return( ptwXY->
status );
56 y2 = a * ptwXY->
points[length-1].
y;
58 for( i = length - 2; i >= 0; i-- ) {
62 if( ( status =
ptwXY_exp_s( ptwXY, x1, y1, z1, x2, y2, z2, 0 ) ) !=
nfu_Okay )
return( status );
76 if( ( x1 == x2 ) || ( y1 == y2 ) )
return(
nfu_Okay );
82 x = 1. / s + x2 - z2 * dx / ( z2 -
z1 );
83 y = ( y1 * ( x2 -
x ) + y2 * ( x - x1 ) ) / dx;
84 z = z1 *
G4Exp( 1 - dy / (
G4Exp( dy ) - 1 ) );
85 zp = ( z2 -
z1 ) / ( y2 - y1 );
87 if( std::fabs( z - zp ) < std::fabs( z * ptwXY->
accuracy ) )
return(
nfu_Okay );
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, z1, x, y, z, level ) ) !=
nfu_Okay )
return( status );
103 int64_t i1, i2,
n1, n2,
n;
105 double accuracy = ptwXY1->
accuracy, yMin, yMax,
c,
y,
dy;
117 if( ( n1 == 0 ) || ( n2 == 0 ) ) {
122 if( ( n1 == 1 ) || ( n2 == 1 ) ) {
127 if( accuracy < ptwXY2->accuracy ) accuracy = ptwXY2->
accuracy;
131 if( n > 1000 ) mode = -1;
133 if( n > 100000 ) mode = -1;
136 yMin = f1->
points[0].
x + f2->points[0].x;
137 yMax = f1->
points[n1 - 1].
x + f2->points[n2 - 1].x;
142 dy = ( yMax - yMin ) / 400;
143 for( y = yMin + dy; y < yMax; y +=
dy ) {
148 for( i1 = 0; i1 <
n1; i1++ ) {
149 for( i2 = 0; i2 < n2; i2++ ) {
150 y = yMin + ( f1->
points[i1].
x - f1->
points[0].
x ) + ( f2->points[i2].x - f2->points[0].x );
151 if( y <= yMin )
continue;
152 if( y >= yMax )
continue;
160 for( i1 = convolute->length - 1; i1 > 0; i1-- ) {
162 convolute->points[i1].x, convolute->points[i1].y, yMin ) ) !=
nfu_Okay )
goto Err;
177 double dx1, dx2, x1MinP, x1Min, x2Max;
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, f2y1p, f2x2p, f2y2p;
184 x2Max = f2->
points[0].
x + ( y - yMin );
187 x1MinP = y - f2->
points[n2 - 1].
x;
188 if( x1Min < x1MinP ) x1Min = x1MinP;
198 i1 = lessThanEqualXPoint.
index;
200 f1y1p = f1y1 = f1->
points[i1].
y;
218 i2 = lessThanEqualXPoint.
index;
219 if( i2 < f2->
length - 1 ) i2++;
221 f2y2p = f2y2 = f2->
points[i2].
y;
235 while( ( i1 <
n1 ) && ( i2 >= 0 ) ) {
240 if( dx1 < dx2 ) mode = 1;
244 if( f2x1p < f2->points[i2].
x ) {
250 *c += ( ( f1y1p + f1y2p ) * ( f2y1p + f2y2p ) + f1y1p * f2y2p + f1y2p * f2y1p ) * dx1;
252 if( i1 ==
n1 )
break;
256 f1y2p = f1y2 = f1->
points[i1].
y;
262 if( ( f1x2p > f1->
points[i1].
x ) || ( dx1 == dx2 ) ) {
268 *c += ( ( f1y1p + f1y2p ) * ( f2y1p + f2y2p ) + f1y1p * f2y2p + f1y2p * f2y1p ) * dx2;
274 f2y1p = f2y1 = f2->
points[i2].
y;
281 f1y2p = f1y2 = f1->
points[i1].
y; }
297 double yMid = 0.5 * ( y1 +
y2 ), cMid = 0.5 * ( c1 + c2 ),
c;
301 if( std::fabs(
c - cMid ) <= convolute->
accuracy * 0.5 * ( std::fabs(
c ) + std::fabs( cMid ) ) )
return(
nfu_Okay );
307 #if defined __cplusplus