13 #if defined __cplusplus
20 static nfu_status
ptwXY_createFromFunctionBisect( ptwXYPoints *ptwXY,
double x1,
double y1,
double x2,
double y2, ptwXY_createFromFunction_callback func,
21 void *argList,
int level,
int checkForRoots,
double eps );
23 ptwXY_createFromFunction_callback func,
void *argList,
double eps );
24 static nfu_status
ptwXY_applyFunction2( ptwXYPoints *ptwXY1,
double y1,
double y2, ptwXYPoint *p1, ptwXYPoint *p2, ptwXY_applyFunction_callback func,
25 void *argList,
int level,
int checkForRoots );
27 ptwXY_applyFunction_callback func,
void *argList );
33 ptwXY1->biSectionMax = ptwXY1->biSectionMax - 1.442695 *
G4Log( ptwXY1->length / oldLength );
34 if( ptwXY1->biSectionMax < 0 ) ptwXY1->biSectionMax = 0;
35 if( ptwXY1->biSectionMax > ptwXY_maxBiSectionMax ) ptwXY1->biSectionMax = ptwXY_maxBiSectionMax;
40 ptwXYPoints *
ptwXY_createFromFunction(
int n,
double *xs, ptwXY_createFromFunction_callback func,
void *argList,
double accuracy,
int checkForRoots,
41 int biSectionMax, nfu_status *status ) {
49 if( n < 2 ) { *status = nfu_tooFewPoints;
return( NULL ); }
50 for( i = 1; i <
n; i++ ) {
51 if( xs[i-1] >= xs[i] ) *status = nfu_XNotAscending;
53 if( *status == nfu_XNotAscending )
return( NULL );
56 if( ( *status = func( x1, &y1, argList ) ) != nfu_Okay )
return( NULL );
57 if( ( ptwXY =
ptwXY_new( ptwXY_interpolationLinLin, NULL, biSectionMax, accuracy, 500, 50, status, 0 ) ) == NULL )
return( NULL );
58 for( i = 1; i <
n; i++ ) {
61 if( ( *status = func( x2, &y2, argList ) ) != nfu_Okay )
goto err;
70 for( i = ptwXY->length - 1, p2 = NULL; i >= 0; i--, p2 = p1 ) {
71 p1 = &(ptwXY->points[i]);
73 if( ( p1->y * p2->y ) < 0. ) {
89 ptwXYPoints *
ptwXY_createFromFunction2( ptwXPoints *xs, ptwXY_createFromFunction_callback func,
void *argList,
double accuracy,
int checkForRoots,
90 int biSectionMax, nfu_status *status ) {
92 return(
ptwXY_createFromFunction( (
int) xs->length, xs->points, func, argList, accuracy, checkForRoots, biSectionMax, status ) );
98 void *argList,
int level,
int checkForRoots,
double eps ) {
103 if( ( x2 - x1 ) < ClosestAllowXFactor *
DBL_EPSILON * ( std::fabs( x1 ) + std::fabs( x2 ) ) )
return( nfu_Okay );
104 if( level >= ptwXY->biSectionMax )
return( nfu_Okay );
105 x = 0.5 * ( x1 + x2 );
106 if( ( status =
ptwXY_interpolatePoint( ptwXY->interpolation, x, &y, x1, y1, x2, y2 ) ) != nfu_Okay )
return( status );
107 if( ( status = func( x, &f, argList ) ) != nfu_Okay )
return( status );
108 if( std::fabs( f - y ) <= 0.8 * std::fabs( f * ptwXY->accuracy ) )
return( nfu_Okay );
117 ptwXY_createFromFunction_callback func,
void *argList,
double eps ) {
123 for( i = 0; i < 16; i++ ) {
124 if( y2 == y1 )
break;
125 x = ( y2 * x1 - y1 * x2 ) / ( y2 - y1 );
126 if( x <= x1 ) x = x1 + 0.1 * ( x2 - x1 );
127 if( x >= x2 ) x = x2 - 0.1 * ( x2 - x1 );
128 if( ( status = func( x, &y, argList ) ) != nfu_Okay )
return( status );
143 nfu_status
ptwXY_applyFunction( ptwXYPoints *ptwXY1, ptwXY_applyFunction_callback func,
void *argList,
int checkForRoots ) {
145 int64_t i, originalLength = ptwXY1->length, notFirstPass = 0;
150 checkForRoots = checkForRoots && ptwXY1->biSectionMax;
151 if( ptwXY1->status != nfu_Okay )
return( ptwXY1->status );
152 if( ptwXY1->interpolation == ptwXY_interpolationOther )
return( nfu_otherInterpolation );
153 if( ptwXY1->interpolation == ptwXY_interpolationFlat )
return( nfu_invalidInterpolation );
155 for( i = originalLength - 1; i >= 0; i-- ) {
156 y1 = ptwXY1->points[i].y;
157 if( ( status = func( &(ptwXY1->points[i]), argList ) ) != nfu_Okay )
return( status );
158 p1 = ptwXY1->points[i];
160 if( ( status =
ptwXY_applyFunction2( ptwXY1, y1, y2, &p1, &p2, func, argList, 0, checkForRoots ) ) != nfu_Okay )
return( status );
172 static nfu_status
ptwXY_applyFunction2( ptwXYPoints *ptwXY1,
double y1,
double y2, ptwXYPoint *p1, ptwXYPoint *p2, ptwXY_applyFunction_callback func,
173 void *argList,
int level,
int checkForRoots ) {
179 if( ( p2->x - p1->x ) < ClosestAllowXFactor *
DBL_EPSILON * ( std::fabs( p1->x ) + std::fabs( p2->x ) ) )
return( nfu_Okay );
180 if( level >= ptwXY1->biSectionMax )
goto checkForZeroCrossing;
181 p.x = 0.5 * ( p1->x + p2->x );
182 if( ( status =
ptwXY_interpolatePoint( ptwXY1->interpolation, p.x, &y, p1->x, y1, p2->x, y2 ) ) != nfu_Okay )
return( status );
184 if( ( status = func( &p, argList ) ) != nfu_Okay )
return( status );
185 if( std::fabs( ( p.x - p1->x ) * ( p2->y - p1->y ) + ( p2->x - p1->x ) * ( p1->y - p.y ) ) <= 0.8 * std::fabs( ( p2->x - p1->x ) * p.y * ptwXY1->accuracy ) )
186 goto checkForZeroCrossing;
187 if( ( status =
ptwXY_setValueAtX( ptwXY1, p.x, p.y ) ) != nfu_Okay )
return( status );
188 if( ( status =
ptwXY_applyFunction2( ptwXY1, y1, y, p1, &p, func, argList, level + 1, checkForRoots ) ) )
return( status );
189 return(
ptwXY_applyFunction2( ptwXY1, y, y2, &p, p2, func, argList, level + 1, checkForRoots ) );
191 checkForZeroCrossing:
199 ptwXY_applyFunction_callback func,
void *argList ) {
202 double y, x1 = p1->x, x2 = p2->x, nY1 = p1->y, nY2 = p2->y, refY = 0.5 * ( std::fabs( p1->y ) + std::fabs( p2->y ) );
206 for( i = 0; i < 6; i++ ) {
207 if( nY2 == nY1 )
break;
208 p.x = ( nY2 * x1 - nY1 * x2 ) / ( nY2 - nY1 );
209 if( p.x <= x1 ) p.x = 0.5 * ( x1 + x2 );
210 if( p.x >= x2 ) p.x = 0.5 * ( x1 + x2 );
211 if( ( status =
ptwXY_interpolatePoint( ptwXY1->interpolation, p.x, &y, p1->x, y1, p2->x, y2 ) ) != nfu_Okay )
return( status );
213 if( ( status = func( &p, argList ) ) != nfu_Okay )
return( status );
214 if( p.y == 0 )
break;
215 if( 0.5 * refY < std::fabs( p.y ) )
break;
216 refY = std::fabs( p.y );
217 if( p1->y * p.y < 0 ) {
230 ptwXYPoints *
ptwXY_fromString(
char const *str, ptwXY_interpolation interpolation, ptwXY_interpolationOtherInfo
const *interpolationOtherInfo,
231 double biSectionMax,
double accuracy,
char **endCharacter, nfu_status *status ) {
233 int64_t numberConverted;
235 ptwXYPoints *ptwXY = NULL;
237 if( ( *status =
nfu_stringToListOfDoubles( str, &numberConverted, &doublePtr, endCharacter ) ) != nfu_Okay )
return( NULL );
238 *status = nfu_oddNumberOfValues;
239 if( ( numberConverted % 2 ) == 0 )
240 ptwXY =
ptwXY_create( interpolation, interpolationOtherInfo, biSectionMax, accuracy, numberConverted, 10, numberConverted / 2, doublePtr, status, 0 );
250 ptwXYPoint *point = ptwXY->points;
251 ptwXYOverflowPoint *overflowPoint;
253 fprintf( f,
"status = %d interpolation = %d length = %d allocatedSize = %d\n",
254 (
int) ptwXY->status, (
int) ptwXY->interpolation, (
int) ptwXY->length, (
int) ptwXY->allocatedSize );
255 fprintf( f,
"userFlag = %d biSectionMax = %.8e accuracy = %.2e minFractional_dx = %.6e\n",
256 ptwXY->userFlag, ptwXY->biSectionMax, ptwXY->accuracy, ptwXY->minFractional_dx );
257 fprintf( f,
"interpolationString = %s\n", ptwXY->interpolationOtherInfo.interpolationString );
258 fprintf( f,
"getValueFunc is NULL = %d. argList is NULL = %d.\n",
259 ptwXY->interpolationOtherInfo.getValueFunc == NULL, ptwXY->interpolationOtherInfo.argList == NULL );
260 fprintf( f,
" overflowLength = %d overflowAllocatedSize = %d mallocFailedSize = %d\n",
261 (
int) ptwXY->overflowLength, (
int) ptwXY->overflowAllocatedSize, (
int) ptwXY->mallocFailedSize );
262 fprintf( f,
" Points data, points = %20p\n", ( printPointersAsNull ? NULL : (
void*)ptwXY->points ) );
263 for( i = 0; i <
n; i++, point++ ) fprintf( f,
" %14.7e %14.7e\n", point->x, point->y );
264 fprintf( f,
" Overflow points data; %20p\n", ( printPointersAsNull ? NULL : (
void*)&(ptwXY->overflowHeader) ) );
265 for( overflowPoint = ptwXY->overflowHeader.next; overflowPoint != &(ptwXY->overflowHeader); overflowPoint = overflowPoint->next ) {
266 fprintf( f,
" %14.7e %14.7e %8d %20p %20p %20p\n", overflowPoint->point.x, overflowPoint->point.y, (
int) overflowPoint->index,
267 (
void*) ( printPointersAsNull ? NULL : overflowPoint ), (
void*) ( printPointersAsNull ? NULL : overflowPoint->prior ),
268 (
void*) ( printPointersAsNull ? NULL : overflowPoint->next ) );
270 fprintf( f,
" Points in order\n" );
271 for( i = 0; i < ptwXY->length; i++ ) {
273 fprintf( f,
" %14.7e %14.7e\n", point->x, point->y );
284 for( i = 0; i < ptwXY->length; i++ ) {
286 fprintf( f, format, point->x, point->y );
297 #if defined __cplusplus
void ptwXY_simpleWrite(ptwXYPoints *ptwXY, FILE *f, char *format)
static nfu_status ptwXY_applyFunction2(ptwXYPoints *ptwXY1, double y1, double y2, ptwXYPoint *p1, ptwXYPoint *p2, ptwXY_applyFunction_callback func, void *argList, int level, int checkForRoots)
ptwXYPoints * ptwXY_createFromFunction2(ptwXPoints *xs, ptwXY_createFromFunction_callback func, void *argList, double accuracy, int checkForRoots, int biSectionMax, nfu_status *status)
static nfu_status ptwXY_createFromFunctionZeroCrossing(ptwXYPoints *ptwXY, double x1, double y1, double x2, double y2, ptwXY_createFromFunction_callback func, void *argList, double eps)
static nfu_status ptwXY_createFromFunctionBisect(ptwXYPoints *ptwXY, double x1, double y1, double x2, double y2, ptwXY_createFromFunction_callback func, void *argList, int level, int checkForRoots, double eps)
ptwXYPoints * ptwXY_create(ptwXY_interpolation interpolation, ptwXY_interpolationOtherInfo const *interpolationOtherInfo, double biSectionMax, double accuracy, int64_t primarySize, int64_t secondarySize, int64_t length, double const *xy, nfu_status *status, int userFlag)
int64_t ptwXY_getNonOverflowLength(ptwXYPoints const *ptwXY)
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
static const G4double eps
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
nfu_status ptwXY_applyFunction(ptwXYPoints *ptwXY1, ptwXY_applyFunction_callback func, void *argList, int checkForRoots)
nfu_status nfu_stringToListOfDoubles(char const *str, int64_t *numberConverted, double **doublePtr, char **endCharacter)
static nfu_status ptwXY_applyFunctionZeroCrossing(ptwXYPoints *ptwXY1, double y1, double y2, ptwXYPoint *p1, ptwXYPoint *p2, ptwXY_applyFunction_callback func, void *argList)
nfu_status ptwXY_setValueAtX_overrideIfClose(ptwXYPoints *ptwXY, double x, double y, double eps, int override)
G4double G4Log(G4double x)
void ptwXY_update_biSectionMax(ptwXYPoints *ptwXY1, double oldLength)
const G4double x[NPOINTSGL]
void ptwXY_showInteralStructure(ptwXYPoints *ptwXY, FILE *f, int printPointersAsNull)
ptwXYPoints * ptwXY_createFromFunction(int n, double *xs, ptwXY_createFromFunction_callback func, void *argList, double accuracy, int checkForRoots, int biSectionMax, nfu_status *status)
ptwXYPoint * ptwXY_getPointAtIndex(ptwXYPoints *ptwXY, int64_t index)
void ptwXY_simplePrint(ptwXYPoints *ptwXY, char *format)
ptwXYPoints * ptwXY_new(ptwXY_interpolation interpolation, ptwXY_interpolationOtherInfo const *interpolationOtherInfo, double biSectionMax, double accuracy, int64_t primarySize, int64_t secondarySize, nfu_status *status, int userFlag)
nfu_status ptwXY_interpolatePoint(ptwXY_interpolation interpolation, double x, double *y, double x1, double y1, double x2, double y2)
ptwXYPoints * ptwXY_fromString(char const *str, ptwXY_interpolation interpolation, ptwXY_interpolationOtherInfo const *interpolationOtherInfo, double biSectionMax, double accuracy, char **endCharacter, nfu_status *status)