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