diff --git a/src/dualplot.cpp b/src/dualplot.cpp
index 8e3abf106ef726d0c691be18caffe5c02627f776..c8dca908b7217171c599e1eaef742a3efb353a00 100644
--- a/src/dualplot.cpp
+++ b/src/dualplot.cpp
@@ -370,58 +370,25 @@ void NPlot::Box()
 }
 
 void NPlot::Line( const int lstyle, 
-		  const vector<double> x, const vector<double> y,
+		  const vector<double> xp, const vector<double> yp,
 		  const vector<double>* z,
 		  const string xco, const string yco,
                   const string info )
 {
-    uint np=0, nxl=0, nxh=0, nyl=0, nyh=0;
-    uint n=x.size();
-    if (n!=y.size())
+    // Checks:
+    uint np=xp.size();
+    if ( np<0 )
+        throw string( "PROG ERR NPLot::Line x.size=0" );
+    if ( np!=yp.size() )
         throw string( "PROG ERR NPLot::Line x.size<>y.size" );
-    vector<double> xp, yp;
 
-    for (uint i=0; i<n; i++) {
-        if(!X.force && !X.R.contained(x[i])) {
-            if (x[i]<=X.R.inf) nxl++;
-            if (x[i]>=X.R.sup) nxh++;
-            continue;
-        }
-        if(!Y.force && !Y.R.contained(y[i])) {
-            if (y[i]<=Y.R.inf) nyl++;
-            if (y[i]>=Y.R.sup) nyh++;
-            continue;
-        }
-        if( np==maxpoints && np<n) {
-            printf("reached maxpoints at %g\n", x[i]);
-        }
-        xp.push_back( x[i] );
-        yp.push_back( y[i] );
-        np = xp.size();
-        if( np > maxpoints ) // continue in large steps
-            i += maxpoints;
-    }
-    if ( np==0 ) {
-        printf( "%s : all %d points outside plot range:\n",
-                info.c_str(), n );
-        if( nxl )
-            printf( "  %d points < xmin\n", nxl );
-        if( nxh )
-            printf( "  %d points > xmax\n", nxh );
-        if( nyl )
-            printf( "  %d points < ymin\n", nyl );
-        if( nyh )
-            printf( "  %d points > ymax\n", nyh );
-        return;
-    }
-
-    // save to file:
+    // Save to file:
     char gp_fnam[40];
     sprintf( gp_fnam, "/tmp/%s-%03d.gnu", getenv("LOGNAME"), gp_fno++);
     if (gp_fnames!="") gp_fnames += ", ";
     gp_fnames += string("\"") + gp_fnam + "\"" + 
-//		string(" title \"") + xco + string (" -> ") 
-//		                    + yco + string ("\"");
+        //		string(" title \"") + xco + string (" -> ") 
+        //		                    + yco + string ("\"");
         " notitle";
     if (lstyle<0)
         gp_fnames += " with lines";
@@ -432,7 +399,7 @@ void NPlot::Line( const int lstyle,
         fprintf(gp_fd, "%16.8g %16.8g\n", xp[i], yp[i]);
     fclose(gp_fd);
 	
-    // plot command:
+    // Plot command:
     char aux[80];
     sprintf( aux, "[%12.8g:%12.8g] [%12.8g:%12.8g] ", 
              X.R.inf, X.R.sup, Y.R.inf, Y.R.sup);
diff --git a/src/dualplot.h b/src/dualplot.h
index f357fe308645283235af61e69da51f3c5dca0e5b..bd7bea8a425ccf091778eea3f9a5cc543225ef74 100644
--- a/src/dualplot.h
+++ b/src/dualplot.h
@@ -37,6 +37,7 @@ class CAxis {
 
 namespace NPlot { // additional declarations in dualplot.cpp
     extern class CAxis X, Y;
+    extern uint maxpoints;
 
     void Multiask();
 
diff --git a/src/edif.cpp b/src/edif.cpp
index 53f062b0ac32ddc3971f4b81ff2bf350102fdabf..032a05fdbdecdc0f7e6308dae695f09429de8b00 100644
--- a/src/edif.cpp
+++ b/src/edif.cpp
@@ -1068,21 +1068,63 @@ void NEdif::Plot( bool add )
         for ( int iv=0; iv<JSel.size(); ++iv){
             j = JSel.V[iv];
             if ( fd ) {
-                lstyle = 0;
-                NPlot::Line( lstyle, fd->VS[j].x, fd->VS[j].y, e->Z(j),
-                             e->P()->xco.str(), e->P()->yco.str(),
-                             "data file "+strg(k)+" scan "+strg(j) );
+                const CScan *s;
+                s = &( fd->VS[j] );
+                uint n=s->x.size();
+                if ( n!=s->y.size() )
+                    throw string( "BUG: plot: x.size<>y.size" );
+                // Filter: x,y -> xp,yp if inside frame
+                uint np=0, nxl=0, nxh=0, nyl=0, nyh=0;
+                vector<double> xp, yp;
+                for ( uint i=0; i<n; i++ ) {
+                    if( !NPlot::X.force && !NPlot::X.R.contained(s->x[i]) ) {
+                        if (s->x[i]<=NPlot::X.R.inf) nxl++;
+                        if (s->x[i]>=NPlot::X.R.sup) nxh++;
+                        continue;
+                    }
+                    if( !NPlot::Y.force && !NPlot::Y.R.contained(s->y[i]) ) {
+                        if (s->y[i]<=NPlot::Y.R.inf) nyl++;
+                        if (s->y[i]>=NPlot::Y.R.sup) nyh++;
+                        continue;
+                    }
+                    if( np==NPlot::maxpoints && np<n ) {
+                        printf("reached maxpoints at %g\n", s->x[i]);
+                    }
+                    xp.push_back( s->x[i] );
+                    yp.push_back( s->y[i] );
+                    np = xp.size();
+                    if( np > NPlot::maxpoints ) // continue in large steps
+                        i += NPlot::maxpoints;
+                }
+                if ( np==0 ) {
+                    printf( "file %i spec %i:" 
+                            " all %d points outside plot range:\n",
+                            k, j, n );
+                    if( nxl )
+                        printf( "  %d points < xmin\n", nxl );
+                    if( nxh )
+                        printf( "  %d points > xmax\n", nxh );
+                    if( nyl )
+                        printf( "  %d points < ymin\n", nyl );
+                    if( nyh )
+                        printf( "  %d points > ymax\n", nyh );
+                } else {
+                    lstyle = 0;
+                    NPlot::Line( lstyle, fd->VS[j].x, fd->VS[j].y, e->Z(j),
+                                 e->P()->xco.str(), e->P()->yco.str(),
+                                 "data file "+strg(k)+" scan "+strg(j) );
+                }
             } else if ( fc && false ) {
                 // keep old vectorial version;
                 // might be useful for accelerated plotting
                 int npts = 987;
-                S.x.resize( npts );
+                vector<double> xp(npts), yp;
                 for( uint i=0; i<npts; ++i )
-                    S.x[i] = NPlot::X.R.plotcoord2value(
+                    xp[i] = NPlot::X.R.plotcoord2value(
                         NPlot::X.log, i/(npts-1.0) );
-                fc->T->tree_curve_val_vec( &(S.y), S.x, k, j );
+                fc->T->tree_curve_val_vec( &yp, xp, k, j );
                 lstyle = -1;
-                NPlot::Line( lstyle, S.x, S.y, e->Z(j),
+                NPlot::Line( lstyle, xp, yp, e->Z(j),
                              e->P()->xco.str(), e->P()->yco.str(),
                              "curve file "+strg(k)+" scan "+strg(j) );
             } else if ( fc ) {