Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lend/src/ptwXY_functions.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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