Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lend/src/ptwXY_misc.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/lend/src/ptwXY_misc.cc (Version 11.3.0) and /processes/hadronic/models/lend/src/ptwXY_misc.cc (Version 11.0.p1)


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