simplify and accelerate determination of x grid, and make it more precise

pub/src/axis.cpp

@@ -234,16 +234,26 @@ double CAxis::value2ploterror( double v, double dv ) const
 }
 
 
-//! Map linear plot scale (0..1) to application scale (inf..sup).
+//! Map linear plot scale (0..n-1) to application scale (inf..sup).
 
-double CAxis::plotcoord2value( double c ) const
+void CAxis::set_xgrid( vector<double>& x, int n ) const
 {
     if ( !finite() )
-        throw string( "undefined plot range" );
-    if ( logflag ) {
-        return inf * exp( c*log(sup/inf) );
-    } else {
-        return inf + c*(sup-inf);
+        throw "undefined plot range";
+    if ( n<2 )
+        throw "call to set_xgrid with less than two points";
+    if ( inf>=sup )
+        throw "call to set_xgrid with invalid limits";   
+    x.resize( n );
+    x[0] = inf;
+    x[n-1] = sup;
+    for( int i=1; i<n-1; ++i ){
+        if ( logflag ) {
+            double c = ((double)(i)) / ((double)(n-1));
+            x[i] = pow( inf, 1-c ) * pow( sup, c );
+        } else {
+            x[i] = ( (n-1-i)*inf + i*sup ) / (n-1);
+        }
     }
 }
 

pub/src/axis.h

@@ -20,7 +20,7 @@ class CAxis {
     bool contains( double val ) const;
     double value2plotcoord( double v ) const;
     double value2ploterror( double v, double dv ) const;
-    double plotcoord2value( double c ) const;
+    void set_xgrid( vector<double>& x, int n ) const;
     double pc( double v ) const;
     double pcerr( double v, double dv ) const;
     double inf_pos() const;

pub/src/dualplot.cpp

@@ -193,9 +193,9 @@ void CPlot::addSpec( bool as_line, int style_no,
     // Checks:
     uint np=xp.size();
     if ( !np )
-        throw string( "BUG: NPLot::Line x.size=0" );
+        throw "invalid call to CPLot::addSpec: no data points";
     if ( np!=yp.size() )
-        throw string( "BUG: NPLot::Line x.size<>y.size" );
+        throw "invalid call to CPLot::addSpec: x.size<>y.size";
 
     // Prepare for live display, to be shown by showSpecs():
     string gp_fnam = str( format( "/tmp/%s-%d-%03d.gnu" )
@@ -216,9 +216,9 @@ void CPlot::addSpec( bool as_line, int style_no,
             throw "Data point number " + strg(i) + " is invalid: x=" +
                 strg(xp[i]) + ", y=" + strg(yp[i]);
         if( xp[i]<X.inf || xp[i]>X.sup )
-            throw "Plot::Line: x["+strg(i)+"]="+strg(xp[i])+" out of range";
+            throw "CPlot::addSpec: x["+strg(i)+"]="+strg(xp[i])+" out of range";
         if( yp[i]<Y.inf || yp[i]>Y.sup )
-            throw "Plot::Line: y["+strg(i)+"]="+strg(yp[i])+" out of range";
+            throw "CPlot::addSpec: y["+strg(i)+"]="+strg(yp[i])+" out of range";
         if( with_errors && dyp.size() )
             fprintf(gp_fd, "%16.8g %16.8g %16.8g\n", xp[i], yp[i], dyp[i] );
         else
@@ -243,8 +243,7 @@ void CPlot::addSpec( bool as_line, int style_no,
     else
         snprintf( outlin, mLin, "%2d pstyle", style_no+1 );
     ps_accu.push_back( outlin );
-    snprintf( outlin, mLin-2, " %% (%s -> %s)",
-              xco.c_str(), yco.c_str() );
+    snprintf( outlin, mLin-2, " %% (%s -> %s)", xco.c_str(), yco.c_str() );
     strncat( outlin, "\n", mLin );
     ps_accu.push_back( outlin );
     for (uint i=0; i<np; i++) {

pub/src/plot.cpp

@@ -92,9 +92,7 @@ void NPlot::Plot( class CPlot *plot, bool add )
                     if ( fd ){
                         s = fd->VS(j);
                     } else {
-                        S->x.resize( npts );
-                        for ( uint i=0; i<npts; i++ )
-                            S->x[i] = plot->X.plotcoord2value( i/(npts-1.0) );
+                        plot->X.set_xgrid( S->x, npts );
                         fc->curve_val_vec( &(S->y), S->x, k, j );
                         s = S;
                     }
@@ -240,9 +238,8 @@ void NPlot::Plot( class CPlot *plot, bool add )
                 // equidistant grid
 
                 uint npts = plot->equipoints;
-                vector<double> xc(npts), yc, xp, yp;
-                for ( uint i=0; i<npts; ++i )
-                    xc[i] = plot->X.plotcoord2value( i/(npts-1.0) );
+                vector<double> xc, yc, xp, yp;
+                plot->X.set_xgrid( xc, npts );
                 fc->curve_val_vec( &yc, xc, k, j );
                 for ( uint i=0; i<npts; ++i ) {
                     if ( plot->X.contains(xc[i]) && plot->Y.contains(yc[i]) ) {
@@ -260,14 +257,11 @@ void NPlot::Plot( class CPlot *plot, bool add )
                 // - TODO interpolate when crossing the horizontal border
                 // - interpolate at discontinuities     
 
-                // ??? bool   logx = plot->X.logflag;
-
                 vector<double> xc, yc, xn, yn;
                 
                 // start with equidistant grid:
                 uint npts = plot->equipoints;
-                for ( uint i=0; i<npts; ++i )
-                    xc.push_back( plot->X.plotcoord2value( i/(npts-1.0) ) );
+                plot->X.set_xgrid( xc, npts );
 
                 // refinement loop:
                 for ( int iref=0; ; ++iref ) {