| Geant4 Cross Reference |
1 /* 1
2 # <<BEGIN-copyright>>
3 # <<END-copyright>>
4 */
5
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <cmath>
9 #include <float.h>
10
11 #include "ptwXY.h"
12
13 #if defined __cplusplus
14 #include <cmath>
15 #include "G4Log.hh"
16 namespace GIDI {
17 using namespace GIDI;
18 #endif
19
20 static nfu_status ptwXY_createFromFunctionBise
21 void *argList, int level, int checkFor
22 static nfu_status ptwXY_createFromFunctionZero
23 ptwXY_createFromFunction_callback func
24 static nfu_status ptwXY_applyFunction2( ptwXYP
25 void *argList, int level, int checkForRoot
26 static nfu_status ptwXY_applyFunctionZeroCross
27 ptwXY_applyFunction_callback func, void *a
28 /*
29 **********************************************
30 */
31 void ptwXY_update_biSectionMax( ptwXYPoints *p
32
33 ptwXY1->biSectionMax = ptwXY1->biSectionMa
34 if( ptwXY1->biSectionMax < 0 ) ptwXY1->biS
35 if( ptwXY1->biSectionMax > ptwXY_maxBiSect
36 }
37 /*
38 **********************************************
39 */
40 ptwXYPoints *ptwXY_createFromFunction( int n,
41 int biSectionMax, nfu_status *status ) {
42
43 int64_t i;
44 double x1, y1, x2 = 0., y2, eps = ClosestA
45 ptwXYPoints *ptwXY;
46 ptwXYPoint *p1, *p2;
47
48 *status = nfu_Okay;
49 if( n < 2 ) { *status = nfu_tooFewPoints;
50 for( i = 1; i < n; i++ ) {
51 if( xs[i-1] >= xs[i] ) *status = nfu_X
52 }
53 if( *status == nfu_XNotAscending ) return(
54
55 x1 = xs[0];
56 if( ( *status = func( x1, &y1, argList ) )
57 if( ( ptwXY = ptwXY_new( ptwXY_interpolati
58 for( i = 1; i < n; i++ ) {
59 if( ( *status = ptwXY_setValueAtX_over
60 x2 = xs[i];
61 if( ( *status = func( x2, &y2, argList
62 if( ( *status = ptwXY_createFromFuncti
63 x1 = x2;
64 y1 = y2;
65 }
66 if( ( *status = ptwXY_setValueAtX_override
67
68 if( checkForRoots ) {
69 if( ( *status = ptwXY_simpleCoalescePo
70 for( i = ptwXY->length - 1, p2 = NULL;
71 p1 = &(ptwXY->points[i]);
72 if( p2 != NULL ) {
73 if( ( p1->y * p2->y ) < 0. ) {
74 if( ( *status = ptwXY_crea
75 }
76 }
77 }
78 }
79
80 return( ptwXY );
81
82 err:
83 ptwXY_free( ptwXY );
84 return( NULL );
85 }
86 /*
87 **********************************************
88 */
89 ptwXYPoints *ptwXY_createFromFunction2( ptwXPo
90 int biSectionMax, nfu_status *status ) {
91
92 return( ptwXY_createFromFunction( (int) xs
93 }
94 /*
95 **********************************************
96 */
97 static nfu_status ptwXY_createFromFunctionBise
98 void *argList, int level, int checkFor
99
100 nfu_status status;
101 double x, y, f;
102
103 if( ( x2 - x1 ) < ClosestAllowXFactor * DB
104 if( level >= ptwXY->biSectionMax ) return(
105 x = 0.5 * ( x1 + x2 );
106 if( ( status = ptwXY_interpolatePoint( ptw
107 if( ( status = func( x, &f, argList ) ) !=
108 if( std::fabs( f - y ) <= 0.8 * std::fabs(
109 if( ( status = ptwXY_createFromFunctionBis
110 if( ( status = ptwXY_setValueAtX_overrideI
111 return( ptwXY_createFromFunctionBisect( pt
112 }
113 /*
114 **********************************************
115 */
116 static nfu_status ptwXY_createFromFunctionZero
117 ptwXY_createFromFunction_callback func
118
119 //For coverity #63077
120 if ( y2 == y1 ) return ( nfu_badInput );
121
122 int i;
123 double x = 0, y = 0;
124 nfu_status status;
125
126 for( i = 0; i < 16; i++ ) {
127 if( y2 == y1 ) break;
128 x = ( y2 * x1 - y1 * x2 ) / ( y2 - y1
129 if( x <= x1 ) x = x1 + 0.1 * ( x2 - x1
130 if( x >= x2 ) x = x2 - 0.1 * ( x2 - x
131 if( ( status = func( x, &y, argList )
132 if( y == 0 ) break;
133 if( y1 * y < 0 ) {
134 x2 = x;
135 y2 = y; }
136 else {
137 x1 = x;
138 y1 = y;
139 }
140 }
141 return( ptwXY_setValueAtX_overrideIfClose(
142 }
143 /*
144 **********************************************
145 */
146 nfu_status ptwXY_applyFunction( ptwXYPoints *p
147
148 int64_t i, originalLength = ptwXY1->length
149 double y1, y2 = 0;
150 nfu_status status;
151 ptwXYPoint p1, p2;
152
153 checkForRoots = checkForRoots && ptwXY1->b
154 if( ptwXY1->status != nfu_Okay ) return( p
155 if( ptwXY1->interpolation == ptwXY_interpo
156 if( ptwXY1->interpolation == ptwXY_interpo
157 if( ( status = ptwXY_simpleCoalescePoints(
158 for( i = originalLength - 1; i >= 0; i-- )
159 y1 = ptwXY1->points[i].y;
160 if( ( status = func( &(ptwXY1->points[
161 p1 = ptwXY1->points[i];
162 if( notFirstPass ) {
163 if( ( status = ptwXY_applyFunction
164 }
165 notFirstPass = 1;
166 p2 = p1;
167 y2 = y1;
168 }
169 ptwXY_update_biSectionMax( ptwXY1, (double
170 return( status );
171 }
172 /*
173 **********************************************
174 */
175 static nfu_status ptwXY_applyFunction2( ptwXYP
176 void *argList, int level, int checkFor
177
178 double y;
179 ptwXYPoint p;
180 nfu_status status;
181
182 if( ( p2->x - p1->x ) < ClosestAllowXFacto
183 if( level >= ptwXY1->biSectionMax ) goto c
184 p.x = 0.5 * ( p1->x + p2->x );
185 if( ( status = ptwXY_interpolatePoint( ptw
186 p.y = y;
187 if( ( status = func( &p, argList ) ) != nf
188 if( std::fabs( ( p.x - p1->x ) * ( p2->y -
189 goto checkForZeroCrossing;
190 if( ( status = ptwXY_setValueAtX( ptwXY1,
191 if( ( status = ptwXY_applyFunction2( ptwXY
192 return( ptwXY_applyFunction2( ptwXY1, y, y
193
194 checkForZeroCrossing:
195 if( checkForRoots && ( ( p1->y * p2->y ) <
196 return( nfu_Okay );
197 }
198 /*
199 **********************************************
200 */
201 static nfu_status ptwXY_applyFunctionZeroCross
202 ptwXY_applyFunction_callback func, voi
203
204 int i;
205 double y, x1 = p1->x, x2 = p2->x, nY1 = p1
206 ptwXYPoint p;
207 nfu_status status;
208
209 //For coverity #63074
210 if ( nY2 == nY1 ) return ( nfu_badInput );
211
212 for( i = 0; i < 6; i++ ) {
213 if( nY2 == nY1 ) break;
214 p.x = ( nY2 * x1 - nY1 * x2 ) / ( nY2
215 if( p.x <= x1 ) p.x = 0.5 * ( x1 + x2
216 if( p.x >= x2 ) p.x = 0.5 * ( x1 + x2
217 if( ( status = ptwXY_interpolatePoint(
218 p.y = y;
219 if( ( status = func( &p, argList ) ) !
220 if( p.y == 0 ) break;
221 if( 0.5 * refY < std::fabs( p.y ) ) br
222 refY = std::fabs( p.y );
223 if( p1->y * p.y < 0 ) {
224 x2 = p.x;
225 nY2 = p.y; }
226 else {
227 x1 = p.x;
228 nY1 = p.y;
229 }
230 }
231 return( ptwXY_setValueAtX( ptwXY1, p.x, 0.
232 }
233 /*
234 **********************************************
235 */
236 ptwXYPoints *ptwXY_fromString( char const *str
237 double biSectionMax, double accuracy, char
238
239 int64_t numberConverted;
240 double *doublePtr;
241 ptwXYPoints *ptwXY = NULL;
242
243 if( ( *status = nfu_stringToListOfDoubles(
244 *status = nfu_oddNumberOfValues;
245 if( ( numberConverted % 2 ) == 0 )
246 ptwXY = ptwXY_create( interpolation, i
247 nfu_free( doublePtr );
248 return( ptwXY );
249 }
250 /*
251 **********************************************
252 */
253 void ptwXY_showInteralStructure( ptwXYPoints *
254
255 int64_t i, n = ptwXY_getNonOverflowLength(
256 ptwXYPoint *point = ptwXY->points;
257 ptwXYOverflowPoint *overflowPoint;
258
259 fprintf( f, "status = %d interpolation =
260 (int) ptwXY->status, (int) ptwXY->inte
261 fprintf( f, "userFlag = %d biSectionMax =
262 ptwXY->userFlag, ptwXY->biSectionMax,
263 fprintf( f, "interpolationString = %s\n",
264 fprintf( f, "getValueFunc is NULL = %d. ar
265 ptwXY->interpolationOtherInfo.getValue
266 fprintf( f, " overflowLength = %d overfl
267 (int) ptwXY->overflowLength, (int) ptw
268 fprintf( f, " Points data, points = %20p\
269 for( i = 0; i < n; i++, point++ ) fprintf
270 fprintf( f, " Overflow points data; %20p\
271 for( overflowPoint = ptwXY->overflowHeader
272 fprintf( f, " %14.7e %14.7e %8d %20
273 (void*) ( printPointersAsNull ? NULL : ove
274 (void*) ( printPointersAsNull ? NULL : ove
275 }
276 fprintf( f, " Points in order\n" );
277 for( i = 0; i < ptwXY->length; i++ ) {
278 point = ptwXY_getPointAtIndex( ptwXY,
279 fprintf( f, " %14.7e %14.7e\n", poi
280 }
281 }
282 /*
283 **********************************************
284 */
285 void ptwXY_simpleWrite( ptwXYPoints *ptwXY, FI
286
287 int64_t i;
288 ptwXYPoint *point;
289
290 for( i = 0; i < ptwXY->length; i++ ) {
291 point = ptwXY_getPointAtIndex( ptwXY,
292 fprintf( f, format, point->x, point->y
293 }
294 }
295 /*
296 **********************************************
297 */
298 void ptwXY_simplePrint( ptwXYPoints *ptwXY, ch
299
300 ptwXY_simpleWrite( ptwXY, stdout, format )
301 }
302
303 #if defined __cplusplus
304 }
305 #endif
306