diff --git a/pub/lib/plot.cpp b/pub/lib/plot.cpp
index 25657418173c00927927ca7432a4a7221386b522..4722b156889d09193f4b9af3b5db47e2ee3c2132 100644
--- a/pub/lib/plot.cpp
+++ b/pub/lib/plot.cpp
@@ -11,6 +11,7 @@
 #include <boost/format.hpp>
 #include <cfloat>
 #include <cmath>
+#include <limits>
 
 #include "defs.hpp"
 
@@ -48,6 +49,18 @@ void plot2D(class CPlot* plot);
 //*  Subroutines: determine data ranges
 //**************************************************************************************************
 
+//! Adjusts limits inf, sup to include points from dat.
+void update_minmax(const vector<double>& dat, bool logflag, double& inf, double& sup)
+{
+    double i, s;
+    if (logflag)
+	triv::getPosminMax(dat, i, s);
+    else
+	triv::getMinMax(dat, i, s);
+    inf = std::min(inf, i);
+    sup = std::max(sup, s);
+}
+
 //! Loop over selected data files to determine x range.
 
 void determine_Xrange(CPlot* plot, vector<int>& JSel)
@@ -536,30 +549,18 @@ void plot2D(class CPlot* plot)
 
     // Determine plot limits.
     // Intensity limits:
-    double yinf, ysup;
-    for (size_t j = 1; j < m; ++j) {
-        double yyinf, yysup;
-        if (plot->Y.logflag)
-            triv::getPosminMax(fd->VS(j)->y, yyinf, yysup);
-        else
-            triv::getMinMax(fd->VS(j)->y, yyinf, yysup);
-        yinf = std::min(yinf, yyinf);
-        ysup = std::max(ysup, yysup);
-    }
-    if (plot->Y.logflag)
-        triv::getPosminMax(fd->VS(0)->y, yinf, ysup);
-    else
-        triv::getMinMax(fd->VS(0)->y, yinf, ysup);
+    double yinf = +std::numeric_limits<double>::infinity();
+    double ysup = -std::numeric_limits<double>::infinity();
+    for (size_t j = 0; j < m; ++j)
+	update_minmax(fd->VS(j)->y, plot->Y.logflag, yinf, ysup);
     // Limits of coordinates x,z:
     vector<double> xlim, zlim;
-    double xinf, xsup;
-    double zinf, zsup;
+    double xinf = +std::numeric_limits<double>::infinity();
+    double xsup = -std::numeric_limits<double>::infinity();
+    double zinf = +std::numeric_limits<double>::infinity();
+    double zsup = -std::numeric_limits<double>::infinity();
     if (plot->mode2D==0) {
-	zlim = triv::histogram_limits(zval);
-	xlim = triv::histogram_limits(fd->VS(0)->x);
-	xinf = xlim.front();
-	xsup = xlim.back();
-	for (size_t j = 1; j < m; ++j) {
+	for (size_t j = 0; j < m; ++j) {
 	    xlim = triv::histogram_limits(fd->VS(j)->x);
 	    xinf = std::min(xinf, xlim.front());
 	    xsup = std::max(xsup, xlim.back());\
@@ -567,19 +568,8 @@ void plot2D(class CPlot* plot)
 	zinf = zlim.front();
 	zsup = zlim.back();
     } else if (plot->mode2D==1) {
-	if (plot->X.logflag)
-	    triv::getPosminMax(fd->VS(0)->x, xinf, xsup);
-	else
-	    triv::getMinMax(fd->VS(0)->x, xinf, xsup);
-	for (size_t j = 1; j < m; ++j) {
-	    double xi, xs;
-	    if (plot->X.logflag)
-		triv::getPosminMax(fd->VS(0)->x, xi, xs);
-	    else
-		triv::getMinMax(fd->VS(0)->x, xi, xs);
-	    xinf = std::min(xinf, xi);
-	    xsup = std::max(xsup, xs);
-	}
+	for (size_t j = 0; j < m; ++j)
+	    update_minmax(fd->VS(j)->x, plot->X.logflag, xinf, xsup);
 	if (plot->Z.logflag)
 	    triv::getPosminMax(zval, zinf, zsup);
 	else