| 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 #include <nf_Legendre.h> 10 #include <nf_Legendre.h>
11 #include <nf_integration.h> 11 #include <nf_integration.h>
12 12
13 #if defined __cplusplus 13 #if defined __cplusplus
14 #include "G4Log.hh" 14 #include "G4Log.hh"
15 #include "G4Pow.hh" 15 #include "G4Pow.hh"
16 namespace GIDI { 16 namespace GIDI {
17 using namespace GIDI; 17 using namespace GIDI;
18 #endif 18 #endif
19 19
20 typedef struct ptwXY_integrateWithFunctionInfo 20 typedef struct ptwXY_integrateWithFunctionInfo_s {
21 int degree; 21 int degree;
22 ptwXY_createFromFunction_callback func; 22 ptwXY_createFromFunction_callback func;
23 void *argList; 23 void *argList;
24 ptwXY_interpolation interpolation; 24 ptwXY_interpolation interpolation;
25 double x1, x2, y1, y2; 25 double x1, x2, y1, y2;
26 } ptwXY_integrateWithFunctionInfo; 26 } ptwXY_integrateWithFunctionInfo;
27 27
28 static nfu_status ptwXY_integrateWithFunction2 28 static nfu_status ptwXY_integrateWithFunction2( nf_Legendre_GaussianQuadrature_callback integrandFunction, void *argList, double x1,
29 double x2, double *integral ); 29 double x2, double *integral );
30 static nfu_status ptwXY_integrateWithFunction3 30 static nfu_status ptwXY_integrateWithFunction3( double x, double *y, void *argList );
31 /* 31 /*
32 ********************************************** 32 ************************************************************
33 */ 33 */
34 nfu_status ptwXY_f_integrate( ptwXY_interpolat 34 nfu_status ptwXY_f_integrate( ptwXY_interpolation interpolation, double x1, double y1, double x2, double y2, double *value ) {
35 35
36 nfu_status status = nfu_Okay; 36 nfu_status status = nfu_Okay;
37 double r; 37 double r;
38 38
39 *value = 0.; 39 *value = 0.;
40 switch( interpolation ) { 40 switch( interpolation ) {
41 case ptwXY_interpolationLinLin : 41 case ptwXY_interpolationLinLin : /* x linear, y linear */
42 *value = 0.5 * ( y1 + y2 ) * ( x2 - x1 42 *value = 0.5 * ( y1 + y2 ) * ( x2 - x1 );
43 break; 43 break;
44 case ptwXY_interpolationLinLog : 44 case ptwXY_interpolationLinLog : /* x linear, y log */
45 if( ( y1 <= 0. ) || ( y2 <= 0. ) ) { 45 if( ( y1 <= 0. ) || ( y2 <= 0. ) ) {
46 status = nfu_badIntegrationInput; 46 status = nfu_badIntegrationInput; }
47 else { 47 else {
48 r = y2 / y1; 48 r = y2 / y1;
49 if( std::fabs( r - 1. ) < 1e-4 ) { 49 if( std::fabs( r - 1. ) < 1e-4 ) {
50 r = r - 1.; 50 r = r - 1.;
51 *value = y1 * ( x2 - x1 ) / ( 51 *value = y1 * ( x2 - x1 ) / ( 1. + r * ( -0.5 + r * ( 1. / 3. + r * ( -0.25 + .2 * r ) ) ) ); }
52 else { 52 else {
53 *value = ( y2 - y1 ) * ( x2 - 53 *value = ( y2 - y1 ) * ( x2 - x1 ) / G4Log( r );
54 } 54 }
55 } 55 }
56 break; 56 break;
57 case ptwXY_interpolationLogLin : 57 case ptwXY_interpolationLogLin : /* x log, y linear */
58 if( ( x1 <= 0. ) || ( x2 <= 0. ) ) { 58 if( ( x1 <= 0. ) || ( x2 <= 0. ) ) {
59 status = nfu_badIntegrationInput; 59 status = nfu_badIntegrationInput; }
60 else { 60 else {
61 r = x2 / x1; 61 r = x2 / x1;
62 if( std::fabs( r - 1. ) < 1e-4 ) { 62 if( std::fabs( r - 1. ) < 1e-4 ) {
63 r = r - 1.; 63 r = r - 1.;
64 r = r * ( -0.5 + r * ( 1. / 3. 64 r = r * ( -0.5 + r * ( 1. / 3. + r * ( -0.25 + .2 * r ) ) );
65 *value = x1 * ( y2 - y1 ) * r 65 *value = x1 * ( y2 - y1 ) * r / ( 1. + r ) + y2 * ( x2 - x1 ); }
66 else { 66 else {
67 *value = ( y1 - y2 ) * ( x2 - 67 *value = ( y1 - y2 ) * ( x2 - x1 ) / G4Log( r ) + x2 * y2 - x1 * y1;
68 } 68 }
69 } 69 }
70 break; 70 break;
71 case ptwXY_interpolationLogLog : 71 case ptwXY_interpolationLogLog : /* x log, y log */
72 if( ( x1 <= 0. ) || ( x2 <= 0. ) || ( 72 if( ( x1 <= 0. ) || ( x2 <= 0. ) || ( y1 <= 0. ) || ( y2 <= 0. ) ) {
73 status = nfu_badIntegrationInput; 73 status = nfu_badIntegrationInput; }
74 else { 74 else {
75 int i, n; 75 int i, n;
76 double a, z, lx, ly, s, f; 76 double a, z, lx, ly, s, f;
77 77
78 r = y2 / y1; 78 r = y2 / y1;
79 if( std::fabs( r - 1. ) < 1e-4 ) { 79 if( std::fabs( r - 1. ) < 1e-4 ) {
80 ly = ( y2 - y1 ) / y1; 80 ly = ( y2 - y1 ) / y1;
81 ly = ly * ( 1. + ly * ( -0.5 + 81 ly = ly * ( 1. + ly * ( -0.5 + ly * ( 1. / 3. - 0.25 * ly ) ) ); }
82 else { 82 else {
83 ly = G4Log( r ); 83 ly = G4Log( r );
84 } 84 }
85 r = x2 / x1; 85 r = x2 / x1;
86 if( std::fabs( r - 1. ) < 1e-4 ) { 86 if( std::fabs( r - 1. ) < 1e-4 ) {
87 lx = ( x2 - x1 ) / x1; 87 lx = ( x2 - x1 ) / x1;
88 lx = lx * ( 1 + lx * ( -0.5 + 88 lx = lx * ( 1 + lx * ( -0.5 + lx * ( 1. / 3. - 0.25 * lx ) ) ); }
89 else { 89 else {
90 lx = G4Log( r ); 90 lx = G4Log( r );
91 } 91 }
92 a = ly / lx; 92 a = ly / lx;
93 if( std::fabs( r - 1. ) < 1e-3 ) { 93 if( std::fabs( r - 1. ) < 1e-3 ) {
94 z = ( x2 - x1 ) / x1; 94 z = ( x2 - x1 ) / x1;
95 n = (int) a; 95 n = (int) a;
96 if( n > 10 ) n = 12; 96 if( n > 10 ) n = 12;
97 if( n < 4 ) n = 6; 97 if( n < 4 ) n = 6;
98 a = a - n + 1; 98 a = a - n + 1;
99 f = n + 1.; 99 f = n + 1.;
100 for( i = 0, s = 0.; i < n; i++ 100 for( i = 0, s = 0.; i < n; i++, a++, f-- ) s = ( 1. + s ) * a * z / f;
101 *value = y1 * ( x2 - x1 ) * ( 101 *value = y1 * ( x2 - x1 ) * ( 1. + s ); }
102 else { 102 else {
103 *value = y1 * x1 * ( G4Pow::Ge 103 *value = y1 * x1 * ( G4Pow::GetInstance()->powA( r, a + 1. ) - 1. ) / ( a + 1. );
104 } 104 }
105 } 105 }
106 break; 106 break;
107 case ptwXY_interpolationFlat : 107 case ptwXY_interpolationFlat : /* x ?, y flat */
108 *value = y1 * ( x2 - x1 ); 108 *value = y1 * ( x2 - x1 );
109 break; 109 break;
110 case ptwXY_interpolationOther : 110 case ptwXY_interpolationOther :
111 status = nfu_otherInterpolation; 111 status = nfu_otherInterpolation;
112 } 112 }
113 return( status ); 113 return( status );
114 } 114 }
115 /* 115 /*
116 ********************************************** 116 ************************************************************
117 */ 117 */
118 double ptwXY_integrate( ptwXYPoints *ptwXY, do 118 double ptwXY_integrate( ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status ) {
119 119
120 int64_t i, n = ptwXY->length; 120 int64_t i, n = ptwXY->length;
121 double sum = 0., dSum, x, y, x1, x2, y1, y 121 double sum = 0., dSum, x, y, x1, x2, y1, y2, _sign = 1.;
122 ptwXYPoint *point; 122 ptwXYPoint *point;
123 123
124 if( ( *status = ptwXY->status ) != nfu_Oka 124 if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
125 *status = nfu_otherInterpolation; 125 *status = nfu_otherInterpolation;
126 if( ptwXY->interpolation == ptwXY_interpol 126 if( ptwXY->interpolation == ptwXY_interpolationOther ) return( 0. );
127 127
128 if( xMax < xMin ) { 128 if( xMax < xMin ) {
129 x = xMin; 129 x = xMin;
130 xMin = xMax; 130 xMin = xMax;
131 xMax = x; 131 xMax = x;
132 _sign = -1.; 132 _sign = -1.;
133 } 133 }
134 if( n < 2 ) return( 0. ); 134 if( n < 2 ) return( 0. );
135 135
136 if( ( *status = ptwXY_simpleCoalescePoints 136 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( 0. );
137 for( i = 0, point = ptwXY->points; i < n; 137 for( i = 0, point = ptwXY->points; i < n; i++, point++ ) {
138 if( point->x >= xMin ) break; 138 if( point->x >= xMin ) break;
139 } 139 }
140 if( i == n ) return( 0. ); 140 if( i == n ) return( 0. );
141 x2 = point->x; 141 x2 = point->x;
142 y2 = point->y; 142 y2 = point->y;
143 if( i > 0 ) { 143 if( i > 0 ) {
144 if( x2 > xMin ) { 144 if( x2 > xMin ) {
145 x1 = point[-1].x; 145 x1 = point[-1].x;
146 y1 = point[-1].y; 146 y1 = point[-1].y;
147 if( ( *status = ptwXY_interpolateP 147 if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMin, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( 0. );
148 if( x2 > xMax ) { 148 if( x2 > xMax ) {
149 double yMax; 149 double yMax;
150 150
151 if( ( *status = ptwXY_interpol 151 if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMax, &yMax, x1, y1, x2, y2 ) ) != nfu_Okay ) return( 0. );
152 if( ( *status = ptwXY_f_integr 152 if( ( *status = ptwXY_f_integrate( ptwXY->interpolation, xMin, y, xMax, yMax, &sum ) ) != nfu_Okay ) return( 0. );
153 return( sum ); } 153 return( sum ); }
154 else { 154 else {
155 if( ( *status = ptwXY_f_integr 155 if( ( *status = ptwXY_f_integrate( ptwXY->interpolation, xMin, y, x2, y2, &sum ) ) != nfu_Okay ) return( 0. );
156 } 156 }
157 } 157 }
158 } 158 }
159 i++; 159 i++;
160 point++; 160 point++;
161 for( ; i < n; i++, point++ ) { 161 for( ; i < n; i++, point++ ) {
162 x1 = x2; 162 x1 = x2;
163 y1 = y2; 163 y1 = y2;
164 x2 = point->x; 164 x2 = point->x;
165 y2 = point->y; 165 y2 = point->y;
166 if( x2 > xMax ) { 166 if( x2 > xMax ) {
167 if( ( *status = ptwXY_interpolateP 167 if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMax, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( 0. );
168 if( ( *status = ptwXY_f_integrate( 168 if( ( *status = ptwXY_f_integrate( ptwXY->interpolation, x1, y1, xMax, y, &dSum ) ) != nfu_Okay ) return( 0. );
169 sum += dSum; 169 sum += dSum;
170 break; 170 break;
171 } 171 }
172 if( ( *status = ptwXY_f_integrate( ptw 172 if( ( *status = ptwXY_f_integrate( ptwXY->interpolation, x1, y1, x2, y2, &dSum ) ) != nfu_Okay ) return( 0. );
173 sum += dSum; 173 sum += dSum;
174 } 174 }
175 175
176 return( _sign * sum ); 176 return( _sign * sum );
177 } 177 }
178 /* 178 /*
179 ********************************************** 179 ************************************************************
180 */ 180 */
181 double ptwXY_integrateDomain( ptwXYPoints *ptw 181 double ptwXY_integrateDomain( ptwXYPoints *ptwXY, nfu_status *status ) {
182 182
183 if( ( *status = ptwXY->status ) != nfu_Oka 183 if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
184 if( ptwXY->length > 0 ) return( ptwXY_inte 184 if( ptwXY->length > 0 ) return( ptwXY_integrate( ptwXY, ptwXY_getXMin( ptwXY ), ptwXY_getXMax( ptwXY ), status ) );
185 return( 0. ); 185 return( 0. );
186 } 186 }
187 /* 187 /*
188 ********************************************** 188 ************************************************************
189 */ 189 */
190 nfu_status ptwXY_normalize( ptwXYPoints *ptwXY 190 nfu_status ptwXY_normalize( ptwXYPoints *ptwXY ) {
191 /* 191 /*
192 * This function assumes ptwXY_integrateDomai 192 * This function assumes ptwXY_integrateDomain checks status and coalesces the points.
193 */ 193 */
194 194
195 int64_t i; 195 int64_t i;
196 nfu_status status; 196 nfu_status status;
197 double sum = ptwXY_integrateDomain( ptwXY, 197 double sum = ptwXY_integrateDomain( ptwXY, &status );
198 198
199 if( status != nfu_Okay ) return( status ); 199 if( status != nfu_Okay ) return( status );
200 if( sum == 0. ) { 200 if( sum == 0. ) {
201 status = nfu_badNorm; } 201 status = nfu_badNorm; }
202 else { 202 else {
203 for( i = 0; i < ptwXY->length; i++ ) p 203 for( i = 0; i < ptwXY->length; i++ ) ptwXY->points[i].y /= sum;
204 } 204 }
205 return( status ); 205 return( status );
206 } 206 }
207 /* 207 /*
208 ********************************************** 208 ************************************************************
209 */ 209 */
210 double ptwXY_integrateDomainWithWeight_x( ptwX 210 double ptwXY_integrateDomainWithWeight_x( ptwXYPoints *ptwXY, nfu_status *status ) {
211 211
212 if( ( *status = ptwXY->status ) != nfu_Oka 212 if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
213 if( ptwXY->length < 2 ) return( 0. ); 213 if( ptwXY->length < 2 ) return( 0. );
214 return( ptwXY_integrateWithWeight_x( ptwXY 214 return( ptwXY_integrateWithWeight_x( ptwXY, ptwXY_getXMin( ptwXY ), ptwXY_getXMax( ptwXY ), status ) );
215 } 215 }
216 /* 216 /*
217 ********************************************** 217 ************************************************************
218 */ 218 */
219 double ptwXY_integrateWithWeight_x( ptwXYPoint 219 double ptwXY_integrateWithWeight_x( ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status ) {
220 220
221 int64_t i, n = ptwXY->length; 221 int64_t i, n = ptwXY->length;
222 double sum = 0., x, y, x1, x2, y1, y2, _si 222 double sum = 0., x, y, x1, x2, y1, y2, _sign = 1.;
223 ptwXYPoint *point; 223 ptwXYPoint *point;
224 224
225 if( ( *status = ptwXY->status ) != nfu_Oka 225 if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
226 *status = nfu_unsupportedInterpolation; 226 *status = nfu_unsupportedInterpolation;
227 if( ( ptwXY->interpolation != ptwXY_interp 227 if( ( ptwXY->interpolation != ptwXY_interpolationLinLin ) &&
228 ( ptwXY->interpolation != ptwXY_interp 228 ( ptwXY->interpolation != ptwXY_interpolationFlat ) ) return( 0. );
229 229
230 if( n < 2 ) return( 0. ); 230 if( n < 2 ) return( 0. );
231 if( xMax < xMin ) { 231 if( xMax < xMin ) {
232 x = xMin; 232 x = xMin;
233 xMin = xMax; 233 xMin = xMax;
234 xMax = x; 234 xMax = x;
235 _sign = -1.; 235 _sign = -1.;
236 } 236 }
237 237
238 if( ( *status = ptwXY_simpleCoalescePoints 238 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( 0. );
239 for( i = 0, point = ptwXY->points; i < n; 239 for( i = 0, point = ptwXY->points; i < n; ++i, ++point ) {
240 if( point->x >= xMin ) break; 240 if( point->x >= xMin ) break;
241 } 241 }
242 if( i == n ) return( 0. ); 242 if( i == n ) return( 0. );
243 x2 = point->x; 243 x2 = point->x;
244 y2 = point->y; 244 y2 = point->y;
245 if( i > 0 ) { 245 if( i > 0 ) {
246 if( x2 > xMin ) { 246 if( x2 > xMin ) {
247 if( ( *status = ptwXY_interpolateP 247 if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMin, &y, point[-1].x, point[-1].y, x2, y2 ) ) != nfu_Okay ) return( 0. );
248 x2 = xMin; 248 x2 = xMin;
249 y2 = y; 249 y2 = y;
250 --i; 250 --i;
251 --point; 251 --point;
252 } 252 }
253 } 253 }
254 ++i; 254 ++i;
255 ++point; 255 ++point;
256 for( ; i < n; ++i, ++point ) { 256 for( ; i < n; ++i, ++point ) {
257 x1 = x2; 257 x1 = x2;
258 y1 = y2; 258 y1 = y2;
259 x2 = point->x; 259 x2 = point->x;
260 y2 = point->y; 260 y2 = point->y;
261 if( x2 > xMax ) { 261 if( x2 > xMax ) {
262 if( ( *status = ptwXY_interpolateP 262 if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMax, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( 0. );
263 x2 = xMax; 263 x2 = xMax;
264 y2 = y; 264 y2 = y;
265 } 265 }
266 switch( ptwXY->interpolation ) { 266 switch( ptwXY->interpolation ) {
267 case ptwXY_interpolationFlat : 267 case ptwXY_interpolationFlat :
268 sum += ( x2 - x1 ) * y1 * 3 * ( x1 268 sum += ( x2 - x1 ) * y1 * 3 * ( x1 + x2 );
269 break; 269 break;
270 case ptwXY_interpolationLinLin : 270 case ptwXY_interpolationLinLin :
271 sum += ( x2 - x1 ) * ( y1 * ( 2 * 271 sum += ( x2 - x1 ) * ( y1 * ( 2 * x1 + x2 ) + y2 * ( x1 + 2 * x2 ) );
272 break; 272 break;
273 default : /* Only to stop compil 273 default : /* Only to stop compilers from complaining. */
274 break; 274 break;
275 } 275 }
276 if( x2 == xMax ) break; 276 if( x2 == xMax ) break;
277 } 277 }
278 278
279 return( _sign * sum / 6 ); 279 return( _sign * sum / 6 );
280 } 280 }
281 /* 281 /*
282 ********************************************** 282 ************************************************************
283 */ 283 */
284 double ptwXY_integrateDomainWithWeight_sqrt_x( 284 double ptwXY_integrateDomainWithWeight_sqrt_x( ptwXYPoints *ptwXY, nfu_status *status ) {
285 285
286 if( ( *status = ptwXY->status ) != nfu_Oka 286 if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
287 if( ptwXY->length < 2 ) return( 0. ); 287 if( ptwXY->length < 2 ) return( 0. );
288 return( ptwXY_integrateWithWeight_sqrt_x( 288 return( ptwXY_integrateWithWeight_sqrt_x( ptwXY, ptwXY_getXMin( ptwXY ), ptwXY_getXMax( ptwXY ), status ) );
289 } 289 }
290 /* 290 /*
291 ********************************************** 291 ************************************************************
292 */ 292 */
293 double ptwXY_integrateWithWeight_sqrt_x( ptwXY 293 double ptwXY_integrateWithWeight_sqrt_x( ptwXYPoints *ptwXY, double xMin, double xMax, nfu_status *status ) {
294 294
295 int64_t i, n = ptwXY->length; 295 int64_t i, n = ptwXY->length;
296 double sum = 0., x, y, x1, x2, y1, y2, _si 296 double sum = 0., x, y, x1, x2, y1, y2, _sign = 1., sqrt_x1, sqrt_x2, inv_apb, c;
297 ptwXYPoint *point; 297 ptwXYPoint *point;
298 298
299 if( ( *status = ptwXY->status ) != nfu_Oka 299 if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
300 *status = nfu_unsupportedInterpolation; 300 *status = nfu_unsupportedInterpolation;
301 if( ( ptwXY->interpolation != ptwXY_interp 301 if( ( ptwXY->interpolation != ptwXY_interpolationLinLin ) &&
302 ( ptwXY->interpolation != ptwXY_interp 302 ( ptwXY->interpolation != ptwXY_interpolationFlat ) ) return( 0. );
303 303
304 if( n < 2 ) return( 0. ); 304 if( n < 2 ) return( 0. );
305 if( xMax < xMin ) { 305 if( xMax < xMin ) {
306 x = xMin; 306 x = xMin;
307 xMin = xMax; 307 xMin = xMax;
308 xMax = x; 308 xMax = x;
309 _sign = -1.; 309 _sign = -1.;
310 } 310 }
311 311
312 if( ( *status = ptwXY_simpleCoalescePoints 312 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( 0. );
313 for( i = 0, point = ptwXY->points; i < n; 313 for( i = 0, point = ptwXY->points; i < n; ++i, ++point ) {
314 if( point->x >= xMin ) break; 314 if( point->x >= xMin ) break;
315 } 315 }
316 if( i == n ) return( 0. ); 316 if( i == n ) return( 0. );
317 x2 = point->x; 317 x2 = point->x;
318 y2 = point->y; 318 y2 = point->y;
319 if( i > 0 ) { 319 if( i > 0 ) {
320 if( x2 > xMin ) { 320 if( x2 > xMin ) {
321 if( ( *status = ptwXY_interpolateP 321 if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMin, &y, point[-1].x, point[-1].y, x2, y2 ) ) != nfu_Okay ) return( 0. );
322 x2 = xMin; 322 x2 = xMin;
323 y2 = y; 323 y2 = y;
324 --i; 324 --i;
325 --point; 325 --point;
326 } 326 }
327 } 327 }
328 ++i; 328 ++i;
329 ++point; 329 ++point;
330 sqrt_x2 = std::sqrt( x2 ); 330 sqrt_x2 = std::sqrt( x2 );
331 for( ; i < n; ++i, ++point ) { 331 for( ; i < n; ++i, ++point ) {
332 x1 = x2; 332 x1 = x2;
333 y1 = y2; 333 y1 = y2;
334 sqrt_x1 = sqrt_x2; 334 sqrt_x1 = sqrt_x2;
335 x2 = point->x; 335 x2 = point->x;
336 y2 = point->y; 336 y2 = point->y;
337 if( x2 > xMax ) { 337 if( x2 > xMax ) {
338 if( ( *status = ptwXY_interpolateP 338 if( ( *status = ptwXY_interpolatePoint( ptwXY->interpolation, xMax, &y, x1, y1, x2, y2 ) ) != nfu_Okay ) return( 0. );
339 x2 = xMax; 339 x2 = xMax;
340 y2 = y; 340 y2 = y;
341 } 341 }
342 sqrt_x2 = std::sqrt( x2 ); 342 sqrt_x2 = std::sqrt( x2 );
343 inv_apb = sqrt_x1 + sqrt_x2; 343 inv_apb = sqrt_x1 + sqrt_x2;
344 c = 2. * ( sqrt_x1 * sqrt_x2 + x1 + x2 344 c = 2. * ( sqrt_x1 * sqrt_x2 + x1 + x2 );
345 switch( ptwXY->interpolation ) { 345 switch( ptwXY->interpolation ) {
346 case ptwXY_interpolationFlat : 346 case ptwXY_interpolationFlat :
347 sum += ( sqrt_x2 - sqrt_x1 ) * y1 347 sum += ( sqrt_x2 - sqrt_x1 ) * y1 * 2.5 * c;
348 break; 348 break;
349 case ptwXY_interpolationLinLin : 349 case ptwXY_interpolationLinLin :
350 sum += ( sqrt_x2 - sqrt_x1 ) * ( y 350 sum += ( sqrt_x2 - sqrt_x1 ) * ( y1 * ( c + x1 * ( 1. + sqrt_x2 / inv_apb ) ) + y2 * ( c + x2 * ( 1. + sqrt_x1 / inv_apb ) ) );
351 break; 351 break;
352 default : /* Only to stop compil 352 default : /* Only to stop compilers from complaining. */
353 break; 353 break;
354 } 354 }
355 if( x2 == xMax ) break; 355 if( x2 == xMax ) break;
356 } 356 }
357 357
358 return( 2. / 15. * _sign * sum ); 358 return( 2. / 15. * _sign * sum );
359 } 359 }
360 /* 360 /*
361 ********************************************** 361 ************************************************************
362 */ 362 */
363 ptwXPoints *ptwXY_groupOneFunction( ptwXYPoint 363 ptwXPoints *ptwXY_groupOneFunction( ptwXYPoints *ptwXY, ptwXPoints *groupBoundaries, ptwXY_group_normType normType, ptwXPoints *ptwX_norm, nfu_status *status ) {
364 364
365 int64_t i, igs, ngs; 365 int64_t i, igs, ngs;
366 double x1, y1, x2, y2, y2p, xg1, xg2, sum; 366 double x1, y1, x2, y2, y2p, xg1, xg2, sum;
367 ptwXYPoints *f; 367 ptwXYPoints *f;
368 ptwXPoints *groupedData = NULL; 368 ptwXPoints *groupedData = NULL;
369 369
370 if( ( *status = ptwXY_simpleCoalescePoints 370 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( NULL );
371 if( ( *status = groupBoundaries->status ) 371 if( ( *status = groupBoundaries->status ) != nfu_Okay ) return( NULL );
372 *status = nfu_otherInterpolation; 372 *status = nfu_otherInterpolation;
373 if( ptwXY->interpolation == ptwXY_interpol 373 if( ptwXY->interpolation == ptwXY_interpolationOther ) return( NULL );
374 374
375 ngs = ptwX_length( groupBoundaries ) - 1; 375 ngs = ptwX_length( groupBoundaries ) - 1;
376 if( normType == ptwXY_group_normType_norm 376 if( normType == ptwXY_group_normType_norm ) {
377 if( ptwX_norm == NULL ) { 377 if( ptwX_norm == NULL ) {
378 *status = nfu_badNorm; 378 *status = nfu_badNorm;
379 return( NULL ); 379 return( NULL );
380 } 380 }
381 *status = ptwX_norm->status; 381 *status = ptwX_norm->status;
382 if( ptwX_norm->status != nfu_Okay ) re 382 if( ptwX_norm->status != nfu_Okay ) return( NULL );
383 if( ptwX_length( ptwX_norm ) != ngs ) 383 if( ptwX_length( ptwX_norm ) != ngs ) {
384 *status = nfu_badNorm; 384 *status = nfu_badNorm;
385 return( NULL ); 385 return( NULL );
386 } 386 }
387 } 387 }
388 388
389 if( ( f = ptwXY_intersectionWith_ptwX( ptw 389 if( ( f = ptwXY_intersectionWith_ptwX( ptwXY, groupBoundaries, status ) ) == NULL ) return( NULL );
390 if( f->length == 0 ) return( ptwX_createLi 390 if( f->length == 0 ) return( ptwX_createLine( ngs, ngs, 0, 0, status ) );
391 391
392 if( ( groupedData = ptwX_new( ngs, status 392 if( ( groupedData = ptwX_new( ngs, status ) ) == NULL ) goto err;
393 xg1 = groupBoundaries->points[0]; 393 xg1 = groupBoundaries->points[0];
394 x1 = f->points[0].x; 394 x1 = f->points[0].x;
395 y1 = f->points[0].y; 395 y1 = f->points[0].y;
396 for( igs = 0, i = 1; igs < ngs; igs++ ) { 396 for( igs = 0, i = 1; igs < ngs; igs++ ) {
397 xg2 = groupBoundaries->points[igs+1]; 397 xg2 = groupBoundaries->points[igs+1];
398 sum = 0; 398 sum = 0;
399 if( xg2 > x1 ) { 399 if( xg2 > x1 ) {
400 for( ; i < f->length; i++, x1 = x2 400 for( ; i < f->length; i++, x1 = x2, y1 = y2 ) {
401 x2 = f->points[i].x; 401 x2 = f->points[i].x;
402 if( x2 > xg2 ) break; 402 if( x2 > xg2 ) break;
403 y2p = y2 = f->points[i].y; 403 y2p = y2 = f->points[i].y;
404 if( f->interpolation == ptwXY_ 404 if( f->interpolation == ptwXY_interpolationFlat ) y2p = y1;
405 sum += ( y1 + y2p ) * ( x2 - x 405 sum += ( y1 + y2p ) * ( x2 - x1 );
406 } 406 }
407 } 407 }
408 if( sum != 0. ) { 408 if( sum != 0. ) {
409 if( normType == ptwXY_group_normTy 409 if( normType == ptwXY_group_normType_dx ) {
410 sum /= ( xg2 - xg1 ); } 410 sum /= ( xg2 - xg1 ); }
411 else if( normType == ptwXY_group_n 411 else if( normType == ptwXY_group_normType_norm ) {
412 if( ptwX_norm->points[igs] == 412 if( ptwX_norm->points[igs] == 0. ) {
413 *status = nfu_divByZero; 413 *status = nfu_divByZero;
414 goto err; 414 goto err;
415 } 415 }
416 sum /= ptwX_norm->points[igs]; 416 sum /= ptwX_norm->points[igs];
417 } 417 }
418 } 418 }
419 groupedData->points[igs] = 0.5 * sum; 419 groupedData->points[igs] = 0.5 * sum;
420 groupedData->length++; 420 groupedData->length++;
421 xg1 = xg2; 421 xg1 = xg2;
422 } 422 }
423 423
424 ptwXY_free( f ); 424 ptwXY_free( f );
425 return( groupedData ); 425 return( groupedData );
426 426
427 err: 427 err:
428 ptwXY_free( f ); 428 ptwXY_free( f );
429 if( groupedData != NULL ) ptwX_free( group 429 if( groupedData != NULL ) ptwX_free( groupedData );
430 return( NULL ); 430 return( NULL );
431 } 431 }
432 /* 432 /*
433 ********************************************** 433 ************************************************************
434 */ 434 */
435 ptwXPoints *ptwXY_groupTwoFunctions( ptwXYPoin 435 ptwXPoints *ptwXY_groupTwoFunctions( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, ptwXPoints *groupBoundaries, ptwXY_group_normType normType,
436 ptwXPoints *ptwX_norm, nfu_status *sta 436 ptwXPoints *ptwX_norm, nfu_status *status ) {
437 437
438 int64_t i, igs, ngs; 438 int64_t i, igs, ngs;
439 double x1, fy1, gy1, x2, fy2, gy2, fy2p, g 439 double x1, fy1, gy1, x2, fy2, gy2, fy2p, gy2p, xg1, xg2, sum;
440 ptwXYPoints *f = NULL, *ff, *g = NULL, *gg 440 ptwXYPoints *f = NULL, *ff, *g = NULL, *gg = NULL;
441 ptwXPoints *groupedData = NULL; 441 ptwXPoints *groupedData = NULL;
442 442
443 if( ( *status = ptwXY_simpleCoalescePoints 443 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
444 if( ( *status = ptwXY_simpleCoalescePoints 444 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY2 ) ) != nfu_Okay ) return( NULL );
445 if( ( *status = groupBoundaries->status ) 445 if( ( *status = groupBoundaries->status ) != nfu_Okay ) return( NULL );
446 *status = nfu_otherInterpolation; 446 *status = nfu_otherInterpolation;
447 if( ptwXY1->interpolation == ptwXY_interpo 447 if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
448 if( ptwXY2->interpolation == ptwXY_interpo 448 if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( NULL );
449 449
450 ngs = ptwX_length( groupBoundaries ) - 1; 450 ngs = ptwX_length( groupBoundaries ) - 1;
451 if( normType == ptwXY_group_normType_norm 451 if( normType == ptwXY_group_normType_norm ) {
452 if( ptwX_norm == NULL ) { 452 if( ptwX_norm == NULL ) {
453 *status = nfu_badNorm; 453 *status = nfu_badNorm;
454 return( NULL ); 454 return( NULL );
455 } 455 }
456 if( ( *status = ptwX_norm->status ) != 456 if( ( *status = ptwX_norm->status ) != nfu_Okay ) return( NULL );
457 if( ptwX_length( ptwX_norm ) != ngs ) 457 if( ptwX_length( ptwX_norm ) != ngs ) {
458 *status = nfu_badNorm; 458 *status = nfu_badNorm;
459 return( NULL ); 459 return( NULL );
460 } 460 }
461 } 461 }
462 462
463 if( ( ff = ptwXY_intersectionWith_ptwX( pt 463 if( ( ff = ptwXY_intersectionWith_ptwX( ptwXY1, groupBoundaries, status ) ) == NULL ) return( NULL );
464 if( ( gg = ptwXY_intersectionWith_ptwX( pt 464 if( ( gg = ptwXY_intersectionWith_ptwX( ptwXY2, groupBoundaries, status ) ) == NULL ) goto err;
465 if( ( ff->length == 0 ) || ( gg->length == 465 if( ( ff->length == 0 ) || ( gg->length == 0 ) ) {
466 ptwXY_free( ff ); 466 ptwXY_free( ff );
467 ptwXY_free( gg ); 467 ptwXY_free( gg );
468 return( ptwX_createLine( ngs, ngs, 0, 468 return( ptwX_createLine( ngs, ngs, 0, 0, status ) );
469 } 469 }
470 470
471 if( ( *status = ptwXY_tweakDomainsToMutual 471 if( ( *status = ptwXY_tweakDomainsToMutualify( ff, gg, 4, 0 ) ) != nfu_Okay ) goto err;
472 if( ( f = ptwXY_union( ff, gg, status, ptw 472 if( ( f = ptwXY_union( ff, gg, status, ptwXY_union_fill ) ) == NULL ) goto err;
473 if( ( g = ptwXY_union( gg, f, status, ptwX 473 if( ( g = ptwXY_union( gg, f, status, ptwXY_union_fill ) ) == NULL ) goto err;
474 474
475 if( ( groupedData = ptwX_new( ngs, status 475 if( ( groupedData = ptwX_new( ngs, status ) ) == NULL ) goto err;
476 xg1 = groupBoundaries->points[0]; 476 xg1 = groupBoundaries->points[0];
477 x1 = f->points[0].x; 477 x1 = f->points[0].x;
478 fy1 = f->points[0].y; 478 fy1 = f->points[0].y;
479 gy1 = g->points[0].y; 479 gy1 = g->points[0].y;
480 for( igs = 0, i = 1; igs < ngs; igs++ ) { 480 for( igs = 0, i = 1; igs < ngs; igs++ ) {
481 xg2 = groupBoundaries->points[igs+1]; 481 xg2 = groupBoundaries->points[igs+1];
482 sum = 0; 482 sum = 0;
483 if( xg2 > x1 ) { 483 if( xg2 > x1 ) {
484 for( ; i < f->length; i++, x1 = x2 484 for( ; i < f->length; i++, x1 = x2, fy1 = fy2, gy1 = gy2 ) {
485 x2 = f->points[i].x; 485 x2 = f->points[i].x;
486 if( x2 > xg2 ) break; 486 if( x2 > xg2 ) break;
487 fy2p = fy2 = f->points[i].y; 487 fy2p = fy2 = f->points[i].y;
488 if( f->interpolation == ptwXY_ 488 if( f->interpolation == ptwXY_interpolationFlat ) fy2p = fy1;
489 gy2p = gy2 = g->points[i].y; 489 gy2p = gy2 = g->points[i].y;
490 if( g->interpolation == ptwXY_ 490 if( g->interpolation == ptwXY_interpolationFlat ) gy2p = gy1;
491 sum += ( ( fy1 + fy2p ) * ( gy 491 sum += ( ( fy1 + fy2p ) * ( gy1 + gy2p ) + fy1 * gy1 + fy2p * gy2p ) * ( x2 - x1 );
492 } 492 }
493 } 493 }
494 if( sum != 0. ) { 494 if( sum != 0. ) {
495 if( normType == ptwXY_group_normTy 495 if( normType == ptwXY_group_normType_dx ) {
496 sum /= ( xg2 - xg1 ); } 496 sum /= ( xg2 - xg1 ); }
497 else if( normType == ptwXY_group_n 497 else if( normType == ptwXY_group_normType_norm ) {
498 if( ptwX_norm->points[igs] == 498 if( ptwX_norm->points[igs] == 0. ) {
499 *status = nfu_divByZero; 499 *status = nfu_divByZero;
500 goto err; 500 goto err;
501 } 501 }
502 sum /= ptwX_norm->points[igs]; 502 sum /= ptwX_norm->points[igs];
503 } 503 }
504 } 504 }
505 groupedData->points[igs] = sum / 6.; 505 groupedData->points[igs] = sum / 6.;
506 groupedData->length++; 506 groupedData->length++;
507 xg1 = xg2; 507 xg1 = xg2;
508 } 508 }
509 509
510 ptwXY_free( f ); 510 ptwXY_free( f );
511 ptwXY_free( g ); 511 ptwXY_free( g );
512 ptwXY_free( ff ); 512 ptwXY_free( ff );
513 ptwXY_free( gg ); 513 ptwXY_free( gg );
514 return( groupedData ); 514 return( groupedData );
515 515
516 err: 516 err:
517 ptwXY_free( ff ); 517 ptwXY_free( ff );
518 if( gg != NULL ) ptwXY_free( gg ); 518 if( gg != NULL ) ptwXY_free( gg );
519 // Coverity #63063 519 // Coverity #63063
520 if( f != NULL ) ptwXY_free( f ); 520 if( f != NULL ) ptwXY_free( f );
521 if( g != NULL ) ptwXY_free( g ); 521 if( g != NULL ) ptwXY_free( g );
522 if( groupedData != NULL ) ptwX_free( group 522 if( groupedData != NULL ) ptwX_free( groupedData );
523 return( NULL ); 523 return( NULL );
524 } 524 }
525 /* 525 /*
526 ********************************************** 526 ************************************************************
527 */ 527 */
528 ptwXPoints *ptwXY_groupThreeFunctions( ptwXYPo 528 ptwXPoints *ptwXY_groupThreeFunctions( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, ptwXYPoints *ptwXY3, ptwXPoints *groupBoundaries,
529 ptwXY_group_normType normType, ptwXPoi 529 ptwXY_group_normType normType, ptwXPoints *ptwX_norm, nfu_status *status ) {
530 530
531 int64_t i, igs, ngs; 531 int64_t i, igs, ngs;
532 double x1, fy1, gy1, hy1, x2, fy2, gy2, hy 532 double x1, fy1, gy1, hy1, x2, fy2, gy2, hy2, fy2p, gy2p, hy2p, xg1, xg2, sum;
533 ptwXYPoints *f = NULL, *ff, *fff = NULL, * 533 ptwXYPoints *f = NULL, *ff, *fff = NULL, *g = NULL, *gg = NULL, *h = NULL, *hh = NULL;
534 ptwXPoints *groupedData = NULL; 534 ptwXPoints *groupedData = NULL;
535 535
536 if( ( *status = ptwXY_simpleCoalescePoints 536 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
537 if( ( *status = ptwXY_simpleCoalescePoints 537 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY2 ) ) != nfu_Okay ) return( NULL );
538 if( ( *status = ptwXY_simpleCoalescePoints 538 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY3 ) ) != nfu_Okay ) return( NULL );
539 if( ( *status = groupBoundaries->status ) 539 if( ( *status = groupBoundaries->status ) != nfu_Okay ) return( NULL );
540 *status = nfu_otherInterpolation; 540 *status = nfu_otherInterpolation;
541 if( ptwXY1->interpolation == ptwXY_interpo 541 if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
542 if( ptwXY2->interpolation == ptwXY_interpo 542 if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( NULL );
543 if( ptwXY3->interpolation == ptwXY_interpo 543 if( ptwXY3->interpolation == ptwXY_interpolationOther ) return( NULL );
544 544
545 ngs = ptwX_length( groupBoundaries ) - 1; 545 ngs = ptwX_length( groupBoundaries ) - 1;
546 if( normType == ptwXY_group_normType_norm 546 if( normType == ptwXY_group_normType_norm ) {
547 if( ptwX_norm == NULL ) { 547 if( ptwX_norm == NULL ) {
548 *status = nfu_badNorm; 548 *status = nfu_badNorm;
549 return( NULL ); 549 return( NULL );
550 } 550 }
551 if( ( *status = ptwX_norm->status ) != 551 if( ( *status = ptwX_norm->status ) != nfu_Okay ) return( NULL );
552 if( ptwX_length( ptwX_norm ) != ngs ) 552 if( ptwX_length( ptwX_norm ) != ngs ) {
553 *status = nfu_badNorm; 553 *status = nfu_badNorm;
554 return( NULL ); 554 return( NULL );
555 } 555 }
556 } 556 }
557 557
558 if( ( ff = ptwXY_intersectionWith_ptwX( pt 558 if( ( ff = ptwXY_intersectionWith_ptwX( ptwXY1, groupBoundaries, status ) ) == NULL ) return( NULL );
559 if( ( gg = ptwXY_intersectionWith_ptwX( pt 559 if( ( gg = ptwXY_intersectionWith_ptwX( ptwXY2, groupBoundaries, status ) ) == NULL ) goto err;
560 if( ( hh = ptwXY_intersectionWith_ptwX( pt 560 if( ( hh = ptwXY_intersectionWith_ptwX( ptwXY3, groupBoundaries, status ) ) == NULL ) goto err;
561 if( ( ff->length == 0 ) || ( gg->length == 561 if( ( ff->length == 0 ) || ( gg->length == 0 ) || ( hh->length == 0 ) ) return( ptwX_createLine( ngs, ngs, 0, 0, status ) );
562 562
563 if( ( *status = ptwXY_tweakDomainsToMutual 563 if( ( *status = ptwXY_tweakDomainsToMutualify( ff, gg, 4, 0 ) ) != nfu_Okay ) goto err;
564 if( ( *status = ptwXY_tweakDomainsToMutual 564 if( ( *status = ptwXY_tweakDomainsToMutualify( ff, hh, 4, 0 ) ) != nfu_Okay ) goto err;
565 if( ( *status = ptwXY_tweakDomainsToMutual 565 if( ( *status = ptwXY_tweakDomainsToMutualify( gg, hh, 4, 0 ) ) != nfu_Okay ) goto err;
566 if( ( fff = ptwXY_union( ff, gg, status, 566 if( ( fff = ptwXY_union( ff, gg, status, ptwXY_union_fill ) ) == NULL ) goto err;
567 if( ( h = ptwXY_union( hh, fff, status, 567 if( ( h = ptwXY_union( hh, fff, status, ptwXY_union_fill ) ) == NULL ) goto err;
568 if( ( f = ptwXY_union( fff, h, status, 568 if( ( f = ptwXY_union( fff, h, status, ptwXY_union_fill ) ) == NULL ) goto err;
569 if( ( g = ptwXY_union( gg, h, status, 569 if( ( g = ptwXY_union( gg, h, status, ptwXY_union_fill ) ) == NULL ) goto err;
570 570
571 if( ( groupedData = ptwX_new( ngs, status 571 if( ( groupedData = ptwX_new( ngs, status ) ) == NULL ) goto err;
572 xg1 = groupBoundaries->points[0]; 572 xg1 = groupBoundaries->points[0];
573 x1 = f->points[0].x; 573 x1 = f->points[0].x;
574 fy1 = f->points[0].y; 574 fy1 = f->points[0].y;
575 gy1 = g->points[0].y; 575 gy1 = g->points[0].y;
576 hy1 = h->points[0].y; 576 hy1 = h->points[0].y;
577 for( igs = 0, i = 1; igs < ngs; igs++ ) { 577 for( igs = 0, i = 1; igs < ngs; igs++ ) {
578 xg2 = groupBoundaries->points[igs+1]; 578 xg2 = groupBoundaries->points[igs+1];
579 sum = 0; 579 sum = 0;
580 if( xg2 > x1 ) { 580 if( xg2 > x1 ) {
581 for( ; i < f->length; i++, x1 = x2 581 for( ; i < f->length; i++, x1 = x2, fy1 = fy2, gy1 = gy2, hy1 = hy2 ) {
582 x2 = f->points[i].x; 582 x2 = f->points[i].x;
583 if( x2 > xg2 ) break; 583 if( x2 > xg2 ) break;
584 fy2p = fy2 = f->points[i].y; 584 fy2p = fy2 = f->points[i].y;
585 if( f->interpolation == ptwXY_ 585 if( f->interpolation == ptwXY_interpolationFlat ) fy2p = fy1;
586 gy2p = gy2 = g->points[i].y; 586 gy2p = gy2 = g->points[i].y;
587 if( g->interpolation == ptwXY_ 587 if( g->interpolation == ptwXY_interpolationFlat ) gy2p = gy1;
588 hy2p = hy2 = h->points[i].y; 588 hy2p = hy2 = h->points[i].y;
589 if( h->interpolation == ptwXY_ 589 if( h->interpolation == ptwXY_interpolationFlat ) hy2p = hy1;
590 sum += ( ( fy1 + fy2p ) * ( gy 590 sum += ( ( fy1 + fy2p ) * ( gy1 + gy2p ) * ( hy1 + hy2p ) + 2 * fy1 * gy1 * hy1 + 2 * fy2p * gy2p * hy2p ) * ( x2 - x1 );
591 } 591 }
592 } 592 }
593 if( sum != 0. ) { 593 if( sum != 0. ) {
594 if( normType == ptwXY_group_normTy 594 if( normType == ptwXY_group_normType_dx ) {
595 sum /= ( xg2 - xg1 ); } 595 sum /= ( xg2 - xg1 ); }
596 else if( normType == ptwXY_group_n 596 else if( normType == ptwXY_group_normType_norm ) {
597 if( ptwX_norm->points[igs] == 597 if( ptwX_norm->points[igs] == 0. ) {
598 *status = nfu_divByZero; 598 *status = nfu_divByZero;
599 goto err; 599 goto err;
600 } 600 }
601 sum /= ptwX_norm->points[igs]; 601 sum /= ptwX_norm->points[igs];
602 } 602 }
603 } 603 }
604 groupedData->points[igs] = sum / 12.; 604 groupedData->points[igs] = sum / 12.;
605 groupedData->length++; 605 groupedData->length++;
606 xg1 = xg2; 606 xg1 = xg2;
607 } 607 }
608 608
609 ptwXY_free( f ); 609 ptwXY_free( f );
610 ptwXY_free( g ); 610 ptwXY_free( g );
611 ptwXY_free( h ); 611 ptwXY_free( h );
612 ptwXY_free( ff ); 612 ptwXY_free( ff );
613 ptwXY_free( gg ); 613 ptwXY_free( gg );
614 ptwXY_free( hh ); 614 ptwXY_free( hh );
615 ptwXY_free( fff ); 615 ptwXY_free( fff );
616 return( groupedData ); 616 return( groupedData );
617 617
618 err: 618 err:
619 ptwXY_free( ff ); 619 ptwXY_free( ff );
620 if( fff != NULL ) ptwXY_free( fff ); 620 if( fff != NULL ) ptwXY_free( fff );
621 if( gg != NULL ) ptwXY_free( gg ); 621 if( gg != NULL ) ptwXY_free( gg );
622 if( hh != NULL ) ptwXY_free( hh ); 622 if( hh != NULL ) ptwXY_free( hh );
623 if( f != NULL ) ptwXY_free( f ); 623 if( f != NULL ) ptwXY_free( f );
624 if( g != NULL ) ptwXY_free( g ); 624 if( g != NULL ) ptwXY_free( g );
625 if( h != NULL ) ptwXY_free( h ); 625 if( h != NULL ) ptwXY_free( h );
626 if( groupedData != NULL ) ptwX_free( group 626 if( groupedData != NULL ) ptwX_free( groupedData );
627 return( NULL ); 627 return( NULL );
628 } 628 }
629 /* 629 /*
630 ********************************************** 630 ************************************************************
631 */ 631 */
632 ptwXPoints *ptwXY_runningIntegral( ptwXYPoints 632 ptwXPoints *ptwXY_runningIntegral( ptwXYPoints *ptwXY, nfu_status *status ) {
633 633
634 int i; 634 int i;
635 ptwXPoints *runningIntegral = NULL; 635 ptwXPoints *runningIntegral = NULL;
636 double integral = 0., sum; 636 double integral = 0., sum;
637 637
638 if( ( *status = ptwXY_simpleCoalescePoints 638 if( ( *status = ptwXY_simpleCoalescePoints( ptwXY ) ) != nfu_Okay ) return( NULL );
639 if( ( runningIntegral = ptwX_new( ptwXY->l 639 if( ( runningIntegral = ptwX_new( ptwXY->length, status ) ) == NULL ) goto err;
640 640
641 if( ( *status = ptwX_setPointAtIndex( runn 641 if( ( *status = ptwX_setPointAtIndex( runningIntegral, 0, 0. ) ) != nfu_Okay ) goto err;
642 for( i = 1; i < ptwXY->length; i++ ) { 642 for( i = 1; i < ptwXY->length; i++ ) {
643 if( ( *status = ptwXY_f_integrate( ptw 643 if( ( *status = ptwXY_f_integrate( ptwXY->interpolation, ptwXY->points[i-1].x, ptwXY->points[i-1].y,
644 ptwXY->points[i].x, ptwXY->points[ 644 ptwXY->points[i].x, ptwXY->points[i].y, &sum ) ) != nfu_Okay ) goto err;
645 integral += sum; 645 integral += sum;
646 if( ( *status = ptwX_setPointAtIndex( 646 if( ( *status = ptwX_setPointAtIndex( runningIntegral, i, integral ) ) != nfu_Okay ) goto err;
647 } 647 }
648 return( runningIntegral ); 648 return( runningIntegral );
649 649
650 err: 650 err:
651 if( runningIntegral != NULL ) ptwX_free( r 651 if( runningIntegral != NULL ) ptwX_free( runningIntegral );
652 return( NULL ); 652 return( NULL );
653 } 653 }
654 /* 654 /*
655 ********************************************** 655 ************************************************************
656 */ 656 */
657 double ptwXY_integrateWithFunction( ptwXYPoint 657 double ptwXY_integrateWithFunction( ptwXYPoints *ptwXY, ptwXY_createFromFunction_callback func, void *argList,
658 double xMin, double xMax, int degree, 658 double xMin, double xMax, int degree, int recursionLimit, double tolerance, nfu_status *status ) {
659 659
660 int64_t i1, i2, n1 = ptwXY->length; 660 int64_t i1, i2, n1 = ptwXY->length;
661 long evaluations; 661 long evaluations;
662 double integral = 0., integral_, sign = -1 662 double integral = 0., integral_, sign = -1., xa, xb;
663 ptwXY_integrateWithFunctionInfo integrateW 663 ptwXY_integrateWithFunctionInfo integrateWithFunctionInfo;
664 ptwXYPoint *point; 664 ptwXYPoint *point;
665 665
666 if( ( *status = ptwXY->status ) != nfu_Oka 666 if( ( *status = ptwXY->status ) != nfu_Okay ) return( 0. );
667 667
668 if( xMin == xMax ) return( 0. ); 668 if( xMin == xMax ) return( 0. );
669 if( n1 < 2 ) return( 0. ); 669 if( n1 < 2 ) return( 0. );
670 670
671 ptwXY_simpleCoalescePoints( ptwXY ); 671 ptwXY_simpleCoalescePoints( ptwXY );
672 672
673 if( xMin > xMax ) { 673 if( xMin > xMax ) {
674 sign = xMin; 674 sign = xMin;
675 xMin = xMax; 675 xMin = xMax;
676 xMax = sign; 676 xMax = sign;
677 sign = -1.; 677 sign = -1.;
678 } 678 }
679 if( xMin >= ptwXY->points[n1-1].x ) return 679 if( xMin >= ptwXY->points[n1-1].x ) return( 0. );
680 if( xMax <= ptwXY->points[0].x ) return( 0 680 if( xMax <= ptwXY->points[0].x ) return( 0. );
681 681
682 for( i1 = 0; i1 < ( n1 - 1 ); i1++ ) { 682 for( i1 = 0; i1 < ( n1 - 1 ); i1++ ) {
683 if( ptwXY->points[i1+1].x > xMin ) bre 683 if( ptwXY->points[i1+1].x > xMin ) break;
684 } 684 }
685 for( i2 = n1 - 1; i2 > i1; i2-- ) { 685 for( i2 = n1 - 1; i2 > i1; i2-- ) {
686 if( ptwXY->points[i2-1].x < xMax ) bre 686 if( ptwXY->points[i2-1].x < xMax ) break;
687 } 687 }
688 point = &(ptwXY->points[i1]); 688 point = &(ptwXY->points[i1]);
689 689
690 integrateWithFunctionInfo.degree = degree; 690 integrateWithFunctionInfo.degree = degree;
691 integrateWithFunctionInfo.func = func; 691 integrateWithFunctionInfo.func = func;
692 integrateWithFunctionInfo.argList = argLis 692 integrateWithFunctionInfo.argList = argList;
693 integrateWithFunctionInfo.interpolation = 693 integrateWithFunctionInfo.interpolation = ptwXY->interpolation;
694 integrateWithFunctionInfo.x2 = point->x; 694 integrateWithFunctionInfo.x2 = point->x;
695 integrateWithFunctionInfo.y2 = point->y; 695 integrateWithFunctionInfo.y2 = point->y;
696 696
697 xa = xMin; 697 xa = xMin;
698 for( ; i1 < i2; i1++ ) { 698 for( ; i1 < i2; i1++ ) {
699 integrateWithFunctionInfo.x1 = integra 699 integrateWithFunctionInfo.x1 = integrateWithFunctionInfo.x2;
700 integrateWithFunctionInfo.y1 = integra 700 integrateWithFunctionInfo.y1 = integrateWithFunctionInfo.y2;
701 ++point; 701 ++point;
702 integrateWithFunctionInfo.x2 = point-> 702 integrateWithFunctionInfo.x2 = point->x;
703 integrateWithFunctionInfo.y2 = point-> 703 integrateWithFunctionInfo.y2 = point->y;
704 xb = point->x; 704 xb = point->x;
705 if( xb > xMax ) xb = xMax; 705 if( xb > xMax ) xb = xMax;
706 *status = nf_GnG_adaptiveQuadrature( p 706 *status = nf_GnG_adaptiveQuadrature( ptwXY_integrateWithFunction2, ptwXY_integrateWithFunction3, &integrateWithFunctionInfo,
707 xa, xb, recursionLimit, tolerance, 707 xa, xb, recursionLimit, tolerance, &integral_, &evaluations );
708 if( *status != nfu_Okay ) return( 0. ) 708 if( *status != nfu_Okay ) return( 0. );
709 integral += integral_; 709 integral += integral_;
710 xa = xb; 710 xa = xb;
711 } 711 }
712 712
713 return( integral ); 713 return( integral );
714 } 714 }
715 /* 715 /*
716 ********************************************** 716 ************************************************************
717 */ 717 */
718 static nfu_status ptwXY_integrateWithFunction2 718 static nfu_status ptwXY_integrateWithFunction2( nf_Legendre_GaussianQuadrature_callback integrandFunction, void *argList, double x1,
719 double x2, double *integral ) { 719 double x2, double *integral ) {
720 720
721 ptwXY_integrateWithFunctionInfo *integrate 721 ptwXY_integrateWithFunctionInfo *integrateWithFunctionInfo = (ptwXY_integrateWithFunctionInfo *) argList;
722 nfu_status status; 722 nfu_status status;
723 723
724 status = nf_Legendre_GaussianQuadrature( i 724 status = nf_Legendre_GaussianQuadrature( integrateWithFunctionInfo->degree, x1, x2, integrandFunction, argList, integral );
725 return( status ); 725 return( status );
726 } 726 }
727 /* 727 /*
728 ********************************************** 728 ************************************************************
729 */ 729 */
730 static nfu_status ptwXY_integrateWithFunction3 730 static nfu_status ptwXY_integrateWithFunction3( double x, double *y, void *argList ) {
731 731
732 double yf; 732 double yf;
733 ptwXY_integrateWithFunctionInfo *integrate 733 ptwXY_integrateWithFunctionInfo *integrateWithFunctionInfo = (ptwXY_integrateWithFunctionInfo *) argList;
734 nfu_status status; 734 nfu_status status;
735 735
736 if( ( status = ptwXY_interpolatePoint( int 736 if( ( status = ptwXY_interpolatePoint( integrateWithFunctionInfo->interpolation, x, &yf,
737 integrateWithFunctionInfo->x1, int 737 integrateWithFunctionInfo->x1, integrateWithFunctionInfo->y1,
738 integrateWithFunctionInfo->x2, int 738 integrateWithFunctionInfo->x2, integrateWithFunctionInfo->y2 ) ) == nfu_Okay ) {
739 status = integrateWithFunctionInfo->fu 739 status = integrateWithFunctionInfo->func( x, y, integrateWithFunctionInfo->argList );
740 *y *= yf; 740 *y *= yf;
741 } 741 }
742 return( status ); 742 return( status );
743 } 743 }
744 744
745 #if defined __cplusplus 745 #if defined __cplusplus
746 } 746 }
747 #endif 747 #endif
748 748