ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ptwXY_interpolation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ptwXY_interpolation.cc
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)( ptwXYPoints *desc, double x1, double y1, double x2, double y2, int depth );
20 
21 static double ptwXY_flatInterpolationToLinear_eps( double px, double eps );
23 static nfu_status ptwXY_LogLogToLinLin( ptwXYPoints *desc, double x1, double y1, double x2, double y2, int depth );
24 static nfu_status ptwXY_LinLogToLinLin( ptwXYPoints *desc, double x1, double y1, double x2, double y2, int depth );
25 static nfu_status ptwXY_LogLinToLinLin( ptwXYPoints *desc, double x1, double y1, double x2, double y2, int depth );
26 static nfu_status ptwXY_otherToLinLin( ptwXYPoints *desc, double x1, double y1, double x2, double y2, int depth );
27 /*
28 ************************************************************
29 */
30 nfu_status ptwXY_interpolatePoint( ptwXY_interpolation interpolation, double x, double *y, double x1, double y1, double x2, double y2 ) {
31 
32  nfu_status status = nfu_Okay;
33 
34  if( interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
35  if( ( x1 > x2 ) || ( x < x1 ) || ( x > x2 ) ) return( nfu_invalidInterpolation );
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 ) {
47  *y = ( y1 * ( x2 - x ) + y2 * ( x - x1 ) ) / ( x2 - x1 );
48  break;
50  if( ( x <= 0. ) || ( x1 <= 0. ) || ( x2 <= 0. ) ) return( nfu_invalidInterpolation );
51  *y = ( y1 * G4Log( x2 / x ) + y2 * G4Log( x / x1 ) ) / G4Log( x2 / x1 );
52  break;
54  if( ( y1 <= 0. ) || ( y2 <= 0. ) ) return( nfu_invalidInterpolation );
55  *y = G4Exp( ( G4Log( y1 ) * ( x2 - x ) + G4Log( y2 ) * ( x - x1 ) ) / ( x2 - x1 ) );
56  break;
58  if( ( x <= 0. ) || ( x1 <= 0. ) || ( x2 <= 0. ) ) return( nfu_invalidInterpolation );
59  if( ( y1 <= 0. ) || ( y2 <= 0. ) ) return( nfu_invalidInterpolation );
60  *y = G4Exp( ( G4Log( y1 ) * G4Log( x2 / x ) + G4Log( y2 ) * G4Log( x / x1 ) ) / G4Log( x2 / x1 ) );
61  break;
63  *y = y1;
64  break;
65  default :
66  status = nfu_invalidInterpolation;
67  }
68  }
69  return( status );
70 }
71 /*
72 ************************************************************
73 */
74 ptwXYPoints *ptwXY_flatInterpolationToLinear( ptwXYPoints *ptwXY, double lowerEps, double upperEps, nfu_status *status ) {
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( ptwXY ) ) != nfu_Okay ) return( NULL );
84  *status = nfu_invalidInterpolation;
85  if( ptwXY->interpolation != ptwXY_interpolationFlat ) return( NULL );
86  *status = nfu_badInput;
87  if( ( lowerEps < 0 ) || ( upperEps < 0 ) || ( ( lowerEps == 0 ) && ( upperEps == 0 ) ) ) return( NULL );
88  if( ( lowerEps != 0 ) && ( lowerEps < minEps ) ) lowerEps = minEps;
89  if( ( upperEps != 0 ) && ( upperEps < minEps ) ) upperEps = minEps;
90 
91  length = ptwXY->length * ( 1 + ( lowerEps == 0 ? 0 : 1 ) + ( lowerEps == 0 ? 0 : 1 ) );
92  if( ( n = ptwXY_new( ptwXY_interpolationLinLin, NULL, ptwXY->biSectionMax, ptwXY->accuracy, length, ptwXY->overflowLength, status, ptwXY->userFlag ) ) == NULL ) return( NULL );
93 
94  p3 = ptwXY->points;
95  if( ptwXY->length > 0 ) ptwXY_setValueAtX( n, p3->x, p3->y );
96  for( i = 0; i < ptwXY->length; i++, p3++ ) {
97  if( i > 1 ) {
98  if( lowerEps > 0 ) {
99  x = ptwXY_flatInterpolationToLinear_eps( p2->x, -lowerEps );
100  if( x > p1->x ) {
101  if( ( *status = ptwXY_setValueAtX( n, x, p1->y ) ) != nfu_Okay ) goto Err;
102  }
103  }
104  if( lowerEps == 0 ) if( ( *status = ptwXY_setValueAtX( n, p2->x, p1->y ) ) != nfu_Okay ) goto Err;
105  if( upperEps == 0 ) if( ( *status = ptwXY_setValueAtX( n, p2->x, p2->y ) ) != nfu_Okay ) goto Err;
106  if( upperEps > 0 ) {
107  x = ptwXY_flatInterpolationToLinear_eps( p2->x, upperEps );
108  if( x < p3->x ) {
109  if( ( *status = ptwXY_setValueAtX( n, x, p2->y ) ) != nfu_Okay ) goto Err;
110  }
111  }
112  }
113  p1 = p2;
114  p2 = p3;
115  }
116  if( ptwXY->length > 1 ) {
117  if( ( lowerEps != 0 ) && ( p1->y != p2->y ) ) {
118  x = ptwXY_flatInterpolationToLinear_eps( p2->x, -lowerEps );
119  if( x > p1->x ) {
120  if( ( *status = ptwXY_setValueAtX( n, x, p1->y ) ) != nfu_Okay ) goto Err;
121  }
122  }
123  if( ( *status = ptwXY_setValueAtX( n, p2->x, p2->y ) ) != nfu_Okay ) goto Err;
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_eps( double px, double eps ) {
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( ptwXYPoints *ptwXY, ptwXY_interpolation interpolationTo, double accuracy, nfu_status *status ) {
154 /*
155 * This function only works when 'ptwXY->interpolation == interpolationTo' or when interpolationTo is ptwXY_interpolationLinLin.
156 */
157  ptwXYPoints *n1;
158  interpolation_func func = NULL;
159 
160  if( ( *status = ptwXY->status ) != nfu_Okay ) return( NULL );
161  if( ptwXY->interpolation == interpolationTo ) {
162  *status = nfu_Okay;
163  return( ptwXY_clone( ptwXY, status ) ); }
164  else {
165  if( interpolationTo == ptwXY_interpolationLinLin ) {
166  switch( ptwXY->interpolation ) {
168  func = ptwXY_LogLogToLinLin; break;
170  func = ptwXY_LinLogToLinLin; break;
172  func = ptwXY_LogLinToLinLin; break;
174  if( ptwXY->interpolationOtherInfo.getValueFunc != NULL ) func = ptwXY_otherToLinLin;
175  break;
176  case ptwXY_interpolationLinLin : /* Stops compilers from complaining. */
178  break;
179  }
180  }
181  }
183  if( func == NULL ) return( NULL );
184 
185  *status = nfu_Okay;
186  if( ( n1 = ptwXY_cloneToInterpolation( ptwXY, interpolationTo, status ) ) == NULL ) return( NULL );
187  if( accuracy < ptwXY->accuracy ) accuracy = ptwXY->accuracy;
188  n1->accuracy = accuracy;
189 
192  *status = ptwXY_toOtherInterpolation2( n1, ptwXY, func );
194  n1->interpolationOtherInfo.argList = NULL;
195  if( *status != nfu_Okay ) n1 = ptwXY_free( n1 );
196  return( n1 );
197 }
198 /*
199 ************************************************************
200 */
202 
203  nfu_status status;
204  int64_t i;
205  double x1, y1, x2, y2;
206 
207  if( ( status = ptwXY_simpleCoalescePoints( src ) ) != nfu_Okay ) return( status );
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, x2, y2, 0 ) ) != nfu_Okay ) break;
216  }
217  x1 = x2;
218  y1 = y2;
219  }
220  return( status );
221 }
222 /*
223 ************************************************************
224 */
225 static nfu_status ptwXY_LogLogToLinLin( ptwXYPoints *desc, double x1, double y1, double x2, double y2, int depth ) {
226 
227  nfu_status status = nfu_Okay;
228  double x, y, u, u2 = x2 / x1, v2 = y2 / y1, logYXs, logXs = G4Log( u2 ), logYs = G4Log( v2 ), vLin, vLog, w;
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 - logYXs ) * ( v2 - 1 ) );
239  vLog = G4Pow::GetInstance()->powA( u, logYXs );
240  }
241  vLin = ( u2 - u + v2 * ( u - 1 ) ) / ( u2 - 1 );
242  if( std::fabs( vLog - vLin ) <= ( vLog * desc->accuracy ) ) return( status );
243  x = x1 * u;
244  y = y1 * vLog;
245  if( ( status = ptwXY_setValueAtX( desc, x, y ) ) != nfu_Okay ) return( status );
246  if( ( status = ptwXY_LogLogToLinLin( desc, x1, y1, x, y, depth + 1 ) ) != nfu_Okay ) return( status );
247  return( ptwXY_LogLogToLinLin( desc, x, y, x2, y2, depth + 1 ) );
248 }
249 /*
250 ************************************************************
251 */
252 static nfu_status ptwXY_LinLogToLinLin( ptwXYPoints *desc, double x1, double y1, double x2, double y2, int depth ) {
253 
254  nfu_status status = nfu_Okay;
255  double x, y, logYs = G4Log( y2 / y1 ), yLinLin;
256 
257  if( depth > 16 ) return( nfu_Okay );
258  x = ( x2 - x1 ) / ( y2 - y1 ) * ( ( y2 - y1 ) / logYs - y1 ) + x1;
259  y = y1 * G4Exp( logYs / ( x2 - x1 ) * ( x - x1 ) );
260  yLinLin = ( y1 * ( x2 - x ) + y2 * ( x - x1 ) ) / ( x2 - x1 );
261  if( std::fabs( y - yLinLin ) <= ( y * desc->accuracy ) ) return( status );
262  if( ( status = ptwXY_setValueAtX( desc, x, y ) ) != nfu_Okay ) return( status );
263  if( ( status = ptwXY_LinLogToLinLin( desc, x1, y1, x, y, depth + 1 ) ) != nfu_Okay ) return( status );
264  return( ptwXY_LinLogToLinLin( desc, x, y, x2, y2, depth + 1 ) );
265 }
266 /*
267 ************************************************************
268 */
269 static nfu_status ptwXY_LogLinToLinLin( ptwXYPoints *desc, double x1, double y1, double x2, double y2, int depth ) {
270 
271  nfu_status status = nfu_Okay;
272  double x = std::sqrt( x2 * x1 ), y, logXs = G4Log( x2 / x1 ), yLinLin;
273 
274  if( depth > 16 ) return( nfu_Okay );
275 #if 0 /* The next line is very unstable at determineing x. Initial x must be chosen better. */
276  x = ( y1 * x2 - y2 * x1 ) / ( y1 * logXs + ( y2 - y1 ) * ( std::log( x / x1 ) - 1 ) );
277 #endif
278  y = ( y2 - y1 ) * G4Log( x / x1 ) / logXs + y1;
279  yLinLin = ( y1 * ( x2 - x ) + y2 * ( x - x1 ) ) / ( x2 - x1 );
280  if( std::fabs( y - yLinLin ) <= ( y * desc->accuracy ) ) return( status );
281  if( ( status = ptwXY_setValueAtX( desc, x, y ) ) != nfu_Okay ) return( status );
282  if( ( status = ptwXY_LogLinToLinLin( desc, x1, y1, x, y, depth + 1 ) ) != nfu_Okay ) return( status );
283  return( ptwXY_LogLinToLinLin( desc, x, y, x2, y2, depth + 1 ) );
284 }
285 /*
286 ************************************************************
287 */
288 static nfu_status ptwXY_otherToLinLin( ptwXYPoints *desc, double x1, double y1, double x2, double y2, int depth ) {
289 
290  nfu_status status;
291  double x = 0.5 * ( x1 + x2 ), y, yLinLin;
293  void *argList = desc->interpolationOtherInfo.argList;
294 
295  if( depth > 16 ) return( nfu_Okay );
296  if( ( status = getValueFunc( argList, x, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( status );
297  yLinLin = ( y1 * ( x2 - x ) + y2 * ( x - x1 ) ) / ( x2 - x1 );
298  if( std::fabs( y - yLinLin ) <= ( y * desc->accuracy ) ) return( status );
299  if( ( status = ptwXY_setValueAtX( desc, x, y ) ) != nfu_Okay ) return( status );
300  if( ( status = ptwXY_otherToLinLin( desc, x1, y1, x, y, depth + 1 ) ) != nfu_Okay ) return( status );
301  return( ptwXY_otherToLinLin( desc, x, y, x2, y2, depth + 1 ) );
302 }
303 /*
304 ************************************************************
305 */
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 ) ) == NULL ) return( NULL );
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; i++, p++ ) {
322  p->x = ( p->x - xMin ) * inverseDx;
323  p->y = p->y * dx;
324  }
325  n->points[n->length-1].x = 1.; /* Make sure last point is realy 1. */
326  return( n );
327 }
328 /*
329 ************************************************************
330 */
331 ptwXYPoints *ptwXY_fromUnitbase( ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status ) {
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 ) ) == NULL ) return( NULL );
341 
342  dx = xMax - xMin;
343  inverseDx = 1. / dx;
344  length = n->length;
345  for( i = 0, p2 = p = n->points; i < length; ++i, ++p ) {
346  p2->x = p->x * dx + xMin;
347  if( i > 0 ) {
348  if( std::fabs( p2->x - xLast ) <= 10. * DBL_EPSILON * ( std::fabs( p2->x ) + std::fabs( 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; /* Make sure last point is realy xMax. */
358  return( n );
359 }
360 /*
361 ************************************************************
362 */
363 ptwXYPoints *ptwXY_unitbaseInterpolate( double w, double w1, ptwXYPoints *ptwXY1, double w2, ptwXYPoints *ptwXY2, nfu_status *status ) {
364 /*
365 * Should we not be checking the interpolation members???????
366 */
367  int64_t i;
368  ptwXYPoints *n1 = NULL, *n2 = NULL, *a = NULL, *r = NULL;
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, status ) ) == NULL ) return( NULL );
382  if( ( n2 = ptwXY_toUnitbase( ptwXY2, status ) ) == NULL ) goto Err;
383  f = ( w - w1 ) / ( w2 - w1 );
384  g = 1. - f;
385  for( i = 0, p = n1->points; i < n1->length; i++, p++ ) p->y *= g;
386  for( i = 0, p = n2->points; i < n2->length; i++, p++ ) p->y *= f;
387  if( ( a = ptwXY_add_ptwXY( n1, n2, status ) ) == NULL ) goto Err;
388 
389  xMin = g * ptwXY1->points[0].x + f * ptwXY2->points[0].x;
390  xMax = g * ptwXY1->points[ptwXY1->length-1].x + f * ptwXY2->points[ptwXY2->length-1].x;
391  if( ( r = ptwXY_fromUnitbase( a, xMin, xMax, status ) ) == NULL ) goto Err;
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