| Geant4 Cross Reference |
1 /* 1
2 # <<BEGIN-copyright>>
3 # <<END-copyright>>
4 */
5
6 #include <cmath>
7 #include <float.h>
8
9 #include "ptwXY.h"
10
11 #if defined __cplusplus
12 #include "G4Exp.hh"
13 #include "G4Log.hh"
14 #include "G4Pow.hh"
15 namespace GIDI {
16 using namespace GIDI;
17 #endif
18
19 typedef nfu_status (*interpolation_func)( ptwX
20
21 static double ptwXY_flatInterpolationToLinear_
22 static nfu_status ptwXY_toOtherInterpolation2(
23 static nfu_status ptwXY_LogLogToLinLin( ptwXYP
24 static nfu_status ptwXY_LinLogToLinLin( ptwXYP
25 static nfu_status ptwXY_LogLinToLinLin( ptwXYP
26 static nfu_status ptwXY_otherToLinLin( ptwXYPo
27 /*
28 **********************************************
29 */
30 nfu_status ptwXY_interpolatePoint( ptwXY_inter
31
32 nfu_status status = nfu_Okay;
33
34 if( interpolation == ptwXY_interpolationOt
35 if( ( x1 > x2 ) || ( x < x1 ) || ( x > x2
36 if( y1 == y2 ) {
37 *y = y1; }
38 else if( x1 == x2 ) {
39 *y = 0.5 * ( y1 + y2 ); }
40 else if( x == x1 ) {
41 *y = y1; }
42 else if( x == x2 ) {
43 *y = y2; }
44 else {
45 switch( interpolation ) {
46 case ptwXY_interpolationLinLin :
47 *y = ( y1 * ( x2 - x ) + y2 * ( x
48 break;
49 case ptwXY_interpolationLogLin :
50 if( ( x <= 0. ) || ( x1 <= 0. ) ||
51 *y = ( y1 * G4Log( x2 / x ) + y2 *
52 break;
53 case ptwXY_interpolationLinLog :
54 if( ( y1 <= 0. ) || ( y2 <= 0. ) )
55 *y = G4Exp( ( G4Log( y1 ) * ( x2 -
56 break;
57 case ptwXY_interpolationLogLog :
58 if( ( x <= 0. ) || ( x1 <= 0. ) ||
59 if( ( y1 <= 0. ) || ( y2 <= 0. ) )
60 *y = G4Exp( ( G4Log( y1 ) * G4Log(
61 break;
62 case ptwXY_interpolationFlat :
63 *y = y1;
64 break;
65 default :
66 status = nfu_invalidInterpolation;
67 }
68 }
69 return( status );
70 }
71 /*
72 **********************************************
73 */
74 ptwXYPoints *ptwXY_flatInterpolationToLinear(
75
76 int64_t i, length;
77 double x;
78 ptwXYPoints *n;
79 ptwXYPoint *p1 = NULL, *p2 = NULL, *p3;
80
81 #define minEps 5e-16
82
83 if( ( *status = ptwXY_simpleCoalescePoints
84 *status = nfu_invalidInterpolation;
85 if( ptwXY->interpolation != ptwXY_interpol
86 *status = nfu_badInput;
87 if( ( lowerEps < 0 ) || ( upperEps < 0 ) |
88 if( ( lowerEps != 0 ) && ( lowerEps < minE
89 if( ( upperEps != 0 ) && ( upperEps < minE
90
91 length = ptwXY->length * ( 1 + ( lowerEps
92 if( ( n = ptwXY_new( ptwXY_interpolationLi
93
94 p3 = ptwXY->points;
95 if( ptwXY->length > 0 ) ptwXY_setValueAtX(
96 for( i = 0; i < ptwXY->length; i++, p3++ )
97 if( i > 1 ) {
98 if( lowerEps > 0 ) {
99 x = ptwXY_flatInterpolationToL
100 if( x > p1->x ) {
101 if( ( *status = ptwXY_setV
102 }
103 }
104 if( lowerEps == 0 ) if( ( *status
105 if( upperEps == 0 ) if( ( *status
106 if( upperEps > 0 ) {
107 x = ptwXY_flatInterpolationToL
108 if( x < p3->x ) {
109 if( ( *status = ptwXY_setV
110 }
111 }
112 }
113 p1 = p2;
114 p2 = p3;
115 }
116 if( ptwXY->length > 1 ) {
117 if( ( lowerEps != 0 ) && ( p1->y != p2
118 x = ptwXY_flatInterpolationToLinea
119 if( x > p1->x ) {
120 if( ( *status = ptwXY_setValue
121 }
122 }
123 if( ( *status = ptwXY_setValueAtX( n,
124 }
125
126 return( n );
127
128 Err:
129 ptwXY_free( n );
130 return( NULL );
131
132 #undef minEps
133 }
134 /*
135 **********************************************
136 */
137 static double ptwXY_flatInterpolationToLinear_
138
139 double x;
140
141 if( px < 0 ) {
142 x = ( 1 - eps ) * px; }
143 else if( px > 0 ) {
144 x = ( 1 + eps ) * px; }
145 else {
146 x = eps;
147 }
148 return( x );
149 }
150 /*
151 **********************************************
152 */
153 ptwXYPoints *ptwXY_toOtherInterpolation( ptwXY
154 /*
155 * This function only works when 'ptwXY->inte
156 */
157 ptwXYPoints *n1;
158 interpolation_func func = NULL;
159
160 if( ( *status = ptwXY->status ) != nfu_Oka
161 if( ptwXY->interpolation == interpolationT
162 *status = nfu_Okay;
163 return( ptwXY_clone( ptwXY, status ) )
164 else {
165 if( interpolationTo == ptwXY_interpola
166 switch( ptwXY->interpolation ) {
167 case ptwXY_interpolationLogLog :
168 func = ptwXY_LogLogToLinLin; b
169 case ptwXY_interpolationLinLog :
170 func = ptwXY_LinLogToLinLin; b
171 case ptwXY_interpolationLogLin :
172 func = ptwXY_LogLinToLinLin; b
173 case ptwXY_interpolationOther :
174 if( ptwXY->interpolationOtherI
175 break;
176 case ptwXY_interpolationLinLin :
177 case ptwXY_interpolationFlat :
178 break;
179 }
180 }
181 }
182 *status = nfu_unsupportedInterpolationConv
183 if( func == NULL ) return( NULL );
184
185 *status = nfu_Okay;
186 if( ( n1 = ptwXY_cloneToInterpolation( ptw
187 if( accuracy < ptwXY->accuracy ) accuracy
188 n1->accuracy = accuracy;
189
190 n1->interpolationOtherInfo.getValueFunc =
191 n1->interpolationOtherInfo.argList = ptwXY
192 *status = ptwXY_toOtherInterpolation2( n1,
193 n1->interpolationOtherInfo.getValueFunc =
194 n1->interpolationOtherInfo.argList = NULL;
195 if( *status != nfu_Okay ) n1 = ptwXY_free(
196 return( n1 );
197 }
198 /*
199 **********************************************
200 */
201 static nfu_status ptwXY_toOtherInterpolation2(
202
203 nfu_status status;
204 int64_t i;
205 double x1, y1, x2, y2;
206
207 if( ( status = ptwXY_simpleCoalescePoints(
208
209 x1 = src->points[0].x;
210 y1 = src->points[0].y;
211 for( i = 1; i < src->length; i++ ) {
212 x2 = src->points[i].x;
213 y2 = src->points[i].y;
214 if( ( x1 != x2 ) && ( y1 != y2 ) ) {
215 if( ( status = func( desc, x1, y1,
216 }
217 x1 = x2;
218 y1 = y2;
219 }
220 return( status );
221 }
222 /*
223 **********************************************
224 */
225 static nfu_status ptwXY_LogLogToLinLin( ptwXYP
226
227 nfu_status status = nfu_Okay;
228 double x, y, u, u2 = x2 / x1, v2 = y2 / y1
229
230 logYXs = logYs / logXs;
231
232 if( depth > 16 ) return( nfu_Okay );
233 if( std::fabs( logYXs - 1 ) < 1e-5 ) {
234 u = 0.5 * ( 1 + u2 );
235 w = ( logYXs - 1 ) * logXs;
236 vLog = u * ( 1. + w * ( 1 + 0.5 * w )
237 else {
238 u = logYXs * ( u2 - v2 ) / ( ( 1 - log
239 vLog = G4Pow::GetInstance()->powA( u,
240 }
241 vLin = ( u2 - u + v2 * ( u - 1 ) ) / ( u2
242 if( std::fabs( vLog - vLin ) <= ( vLog * d
243 x = x1 * u;
244 y = y1 * vLog;
245 if( ( status = ptwXY_setValueAtX( desc, x,
246 if( ( status = ptwXY_LogLogToLinLin( desc,
247 return( ptwXY_LogLogToLinLin( desc, x, y,
248 }
249 /*
250 **********************************************
251 */
252 static nfu_status ptwXY_LinLogToLinLin( ptwXYP
253
254 nfu_status status = nfu_Okay;
255 double x, y, logYs = G4Log( y2 / y1 ), yLi
256
257 if( depth > 16 ) return( nfu_Okay );
258 x = ( x2 - x1 ) / ( y2 - y1 ) * ( ( y2 - y
259 y = y1 * G4Exp( logYs / ( x2 - x1 ) * ( x
260 yLinLin = ( y1 * ( x2 - x ) + y2 * ( x - x
261 if( std::fabs( y - yLinLin ) <= ( y * desc
262 if( ( status = ptwXY_setValueAtX( desc, x,
263 if( ( status = ptwXY_LinLogToLinLin( desc,
264 return( ptwXY_LinLogToLinLin( desc, x, y,
265 }
266 /*
267 **********************************************
268 */
269 static nfu_status ptwXY_LogLinToLinLin( ptwXYP
270
271 nfu_status status = nfu_Okay;
272 double x = std::sqrt( x2 * x1 ), y, logXs
273
274 if( depth > 16 ) return( nfu_Okay );
275 #if 0 /* The next line is very unstable at det
276 x = ( y1 * x2 - y2 * x1 ) / ( y1 * logXs +
277 #endif
278 y = ( y2 - y1 ) * G4Log( x / x1 ) / logXs
279 yLinLin = ( y1 * ( x2 - x ) + y2 * ( x - x
280 if( std::fabs( y - yLinLin ) <= ( y * desc
281 if( ( status = ptwXY_setValueAtX( desc, x,
282 if( ( status = ptwXY_LogLinToLinLin( desc,
283 return( ptwXY_LogLinToLinLin( desc, x, y,
284 }
285 /*
286 **********************************************
287 */
288 static nfu_status ptwXY_otherToLinLin( ptwXYPo
289
290 nfu_status status;
291 double x = 0.5 * ( x1 + x2 ), y, yLinLin;
292 ptwXY_getValue_callback getValueFunc = des
293 void *argList = desc->interpolationOtherIn
294
295 if( depth > 16 ) return( nfu_Okay );
296 if( ( status = getValueFunc( argList, x, &
297 yLinLin = ( y1 * ( x2 - x ) + y2 * ( x - x
298 if( std::fabs( y - yLinLin ) <= ( y * desc
299 if( ( status = ptwXY_setValueAtX( desc, x,
300 if( ( status = ptwXY_otherToLinLin( desc,
301 return( ptwXY_otherToLinLin( desc, x, y, x
302 }
303 /*
304 **********************************************
305 */
306 ptwXYPoints *ptwXY_toUnitbase( ptwXYPoints *pt
307
308 int64_t i;
309 ptwXYPoints *n;
310 ptwXYPoint *p;
311 double xMin, xMax, dx, inverseDx;
312
313 *status = nfu_tooFewPoints;
314 if( ptwXY->length < 2 ) return( NULL );
315 if( ( n = ptwXY_clone( ptwXY, status ) ) =
316
317 xMin = n->points[0].x;
318 xMax = n->points[n->length-1].x;
319 dx = xMax - xMin;
320 inverseDx = 1. / dx;
321 for( i = 0, p = n->points; i < n->length;
322 p->x = ( p->x - xMin ) * inverseDx;
323 p->y = p->y * dx;
324 }
325 n->points[n->length-1].x = 1.;
326 return( n );
327 }
328 /*
329 **********************************************
330 */
331 ptwXYPoints *ptwXY_fromUnitbase( ptwXYPoints *
332
333 int64_t i, length;
334 ptwXYPoints *n;
335 ptwXYPoint *p, *p2;
336 double dx, inverseDx, xLast = 0.;
337
338 *status = nfu_tooFewPoints;
339 if( ptwXY->length < 2 ) return( NULL );
340 if( ( n = ptwXY_clone( ptwXY, status ) ) =
341
342 dx = xMax - xMin;
343 inverseDx = 1. / dx;
344 length = n->length;
345 for( i = 0, p2 = p = n->points; i < length
346 p2->x = p->x * dx + xMin;
347 if( i > 0 ) {
348 if( std::fabs( p2->x - xLast ) <=
349 --(n->length);
350 continue;
351 }
352 }
353 p2->y = p->y * inverseDx;
354 xLast = p2->x;
355 ++p2;
356 }
357 n->points[n->length-1].x = xMax;
358 return( n );
359 }
360 /*
361 **********************************************
362 */
363 ptwXYPoints *ptwXY_unitbaseInterpolate( double
364 /*
365 * Should we not be checking the interpolatio
366 */
367 int64_t i;
368 ptwXYPoints *n1 = NULL, *n2 = NULL, *a = N
369 ptwXYPoint *p;
370 double f, g, xMin, xMax;
371
372 *status = nfu_XOutsideDomain;
373 if( w <= w1 ) {
374 if( w < w1 ) return( NULL );
375 return( ptwXY_clone( ptwXY1, status )
376 }
377 if( w >= w2 ) {
378 if( w > w2 ) return( NULL );
379 return( ptwXY_clone( ptwXY2, status )
380 }
381 if( ( n1 = ptwXY_toUnitbase( ptwXY1, statu
382 if( ( n2 = ptwXY_toUnitbase( ptwXY2, statu
383 f = ( w - w1 ) / ( w2 - w1 );
384 g = 1. - f;
385 for( i = 0, p = n1->points; i < n1->length
386 for( i = 0, p = n2->points; i < n2->length
387 if( ( a = ptwXY_add_ptwXY( n1, n2, status
388
389 xMin = g * ptwXY1->points[0].x + f * ptwXY
390 xMax = g * ptwXY1->points[ptwXY1->length-1
391 if( ( r = ptwXY_fromUnitbase( a, xMin, xMa
392 ptwXY_free( n1 );
393 ptwXY_free( n2 );
394 ptwXY_free( a );
395 return( r );
396
397 Err:
398 if( n1 != NULL ) ptwXY_free( n1 );
399 if( n2 != NULL ) ptwXY_free( n2 );
400 if( a != NULL ) ptwXY_free( a );
401 return( NULL );
402 }
403
404 #if defined __cplusplus
405 }
406 #endif
407