| Geant4 Cross Reference |
1 /* 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( ptwXYPoi
18 static nfu_status ptwXY_exp_s( ptwXYPoints *pt
19 static nfu_status ptwXY_convolution2( ptwXYPoi
20 static nfu_status ptwXY_convolution3( ptwXYPoi
21 /*
22 **********************************************
23 */
24 nfu_status ptwXY_pow( ptwXYPoints *ptwXY, doub
25
26 return( ptwXY_applyFunction( ptwXY, ptwXY_
27 }
28 /*
29 **********************************************
30 */
31 static nfu_status ptwXY_pow_callback( ptwXYPoi
32
33 nfu_status status = nfu_Okay;
34 double *v = (double *) argList;
35
36 point->y = G4Pow::GetInstance()->powA( poi
37 /* ???? Need to test for valid y-value. */
38 return( status );
39 }
40 /*
41 **********************************************
42 */
43 nfu_status ptwXY_exp( ptwXYPoints *ptwXY, doub
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_interpol
52 if( ptwXY->interpolation == ptwXY_interpol
53
54 if( ( status = ptwXY_simpleCoalescePoints(
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,
63 x2 = x1;
64 y2 = y1;
65 }
66 return( status );
67 }
68 /*
69 **********************************************
70 */
71 static nfu_status ptwXY_exp_s( ptwXYPoints *pt
72
73 nfu_status status;
74 double x, y, dx, dy, z, zp, s;
75
76 if( ( x1 == x2 ) || ( y1 == y2 ) ) return(
77 if( level >= ptwXY->biSectionMax ) return(
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 ) )
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 * p
88 if( ( status = ptwXY_setValueAtX( ptwXY, x
89 if( ( status = ptwXY_exp_s( ptwXY, x, y, z
90 if( ( status = ptwXY_exp_s( ptwXY, x1, y1,
91 return( status );
92 }
93 /*
94 **********************************************
95 */
96 ptwXYPoints *ptwXY_convolution( ptwXYPoints *p
97 /*
98 * Currently, only supports linear-linear int
99 *
100 * This function calculates c(y) = int
101 *
102 */
103 int64_t i1, i2, n1, n2, n;
104 ptwXYPoints *f1 = ptwXY1, *f2 = ptwXY2, *c
105 double accuracy = ptwXY1->accuracy, yMin,
106
107 if( ( *status = ptwXY_simpleCoalescePoints
108 if( ( *status = ptwXY_simpleCoalescePoints
109
110 *status = nfu_unsupportedInterpolation;
111 if( ( ptwXY1->interpolation != ptwXY_inter
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_interpola
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
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_interpo
135
136 yMin = f1->points[0].x + f2->points[0].x;
137 yMax = f1->points[n1 - 1].x + f2->points[n
138
139 if( ( *status = ptwXY_setValueAtX( convolu
140
141 if( mode < 0 ) {
142 dy = ( yMax - yMin ) / 400;
143 for( y = yMin + dy; y < yMax; y += dy
144 if( ( *status = ptwXY_convolution2
145 if( ( *status = ptwXY_setValueAtX(
146 } }
147 else {
148 for( i1 = 0; i1 < n1; i1++ ) {
149 for( i2 = 0; i2 < n2; i2++ ) {
150 y = yMin + ( f1->points[i1].x
151 if( y <= yMin ) continue;
152 if( y >= yMax ) continue;
153 if( ( *status = ptwXY_convolut
154 if( ( *status = ptwXY_setValue
155 }
156 }
157 }
158 if( ( *status = ptwXY_setValueAtX( convolu
159 if( ( *status = ptwXY_simpleCoalescePoints
160 for( i1 = convolute->length - 1; i1 > 0; i
161 if( ( *status = ptwXY_convolution3( co
162 convolute->points[i1].x, convolute
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( ptwXYPoi
175
176 int64_t i1 = 0, i2 = 0, n1 = f1->length, n
177 double dx1, dx2, x1MinP, x1Min, x2Max;
178 double f1x1 = 0, f1y1 = 0, f1x2 = 0, f1
179 double f1x1p, f1y1p, f1x2p, f1y2p, f2x1p,
180 ptwXY_lessEqualGreaterX legx;
181 ptwXYOverflowPoint lessThanEqualXPoint, gr
182 nfu_status status;
183
184 x2Max = f2->points[0].x + ( y - yMin );
185 if( x2Max > f2->points[n2 - 1].x ) x2Max =
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,
192 case ptwXY_lessEqualGreaterX_empty :
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_be
206 if( ( status = ptwXY_interpolatePo
207 }
208 break;
209 }
210
211 switch( legx = ptwXY_getPointsAroundX( f2,
212 case ptwXY_lessEqualGreaterX_empty :
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_be
226 if( ( status = ptwXY_interpolatePo
227 }
228 break;
229 }
230
231 f1x1p = x1Min;
232 f2x2p = x2Max;
233 f1y2p = f1y2;
234 f2y1p = f2y1;
235 while( ( i1 < n1 ) && ( i2 >= 0 ) ) { // L
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 ) {
243 f2x1p = f2x2p - dx1;
244 if( f2x1p < f2->points[i2].x ) {
245 f2x1p = f2x2;
246 f2y1p = f2y2; }
247 else {
248 if( ( status = ptwXY_interpola
249 }
250 *c += ( ( f1y1p + f1y2p ) * ( f2y1
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 ) |
263 f1x2p = f1x2;
264 f1y2p = f1y2; }
265 else {
266 if( ( status = ptwXY_interpola
267 }
268 *c += ( ( f1y1p + f1y2p ) * ( f2y1
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].
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( ptwXYPoi
295
296 nfu_status status;
297 double yMid = 0.5 * ( y1 + y2 ), cMid = 0.
298
299 if( ( y2 - yMid ) <= 1e-5 * ( ptwXY_getXMa
300 if( ( status = ptwXY_convolution2( f1, f2,
301 if( std::fabs( c - cMid ) <= convolute->ac
302 if( ( status = ptwXY_setValueAtX( convolut
303 if( ( status = ptwXY_convolution3( convolu
304 return( ptwXY_convolution3( convolute, f1,
305 }
306
307 #if defined __cplusplus
308 }
309 #endif
310