From 0f82c4fc97437ca623ec3e8f4646ded126a55293 Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (o)" <j.wuttke@fz-juelich.de>
Date: Mon, 20 Mar 2017 12:59:00 +0100
Subject: [PATCH] 2D plot with special treatment for values <0 or >10

---
 pub/lib/plot.cpp          |  13 ++++-
 pub/plot/axis.cpp         |   2 +-
 pub/share/setup1D.ps      |   2 +-
 pub/share/setup2D.ps      |   2 +-
 pub/share/wups17a.ps      | 113 +++++++++++++++++++++++---------------
 pub/trivia/vector_ops.cpp |  14 +++++
 pub/trivia/vector_ops.hpp |   2 +
 7 files changed, 100 insertions(+), 48 deletions(-)

diff --git a/pub/lib/plot.cpp b/pub/lib/plot.cpp
index 2d2e5055..24af67fe 100644
--- a/pub/lib/plot.cpp
+++ b/pub/lib/plot.cpp
@@ -7,6 +7,7 @@
 //! \file  plot.cpp
 //! \brief NPlot: plot data and curves.
 
+#include <cmath>
 #include <algorithm>
 #include <boost/format.hpp>
 
@@ -524,18 +525,26 @@ namespace {
         double xinf = xlim.front();
         double xsup = xlim.back();
         double yinf, ysup;
-        triv::getMinMax(fd->VS(0)->y, yinf, ysup);
+        if (plot->Y.logflag)
+            triv::getPosminMax(fd->VS(0)->y, yinf, ysup);
+        else
+            triv::getMinMax(fd->VS(0)->y, yinf, ysup);
         for (size_t j=1; j<m; ++j) {
             xlim = triv::histogram_limits(fd->VS(j)->x);
             xinf = std::min(xinf, xlim.front());
             xsup = std::max(xsup, xlim.back());
             double yyinf, yysup;
-            triv::getMinMax(fd->VS(j)->y, 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);
         }
         plot->X.set_limits(xinf, xsup);
         plot->Z.set_limits(zlim.front(), zlim.back());
+        if (plot->Y.logflag && !(std::isfinite(yinf) && std::isfinite(ysup)))
+            throw "No nonzero y values";
         plot->Y.set_limits(yinf, ysup);
         // plot
         plot->start_frame2D("Frida version " VERSION);
diff --git a/pub/plot/axis.cpp b/pub/plot/axis.cpp
index 8cb3840f..1b0c2fe9 100644
--- a/pub/plot/axis.cpp
+++ b/pub/plot/axis.cpp
@@ -369,7 +369,7 @@ void CAxis::calc_ticks(vector<double>& Tacks, int* ntpt, double* ticklim) const
         double rlgmax = log10(sup);
         double rlgrel = rlgmax - rlgmin;
         if (!std::isfinite(rlgrel) || rlgrel > 100)
-            throw S("Excessive log range");
+            throw "Excessive log range " + S(inf) + " .. " + S(sup);
 
         // incr := #decades per tack
         int incr;
diff --git a/pub/share/setup1D.ps b/pub/share/setup1D.ps
index 03429028..3778174b 100644
--- a/pub/share/setup1D.ps
+++ b/pub/share/setup1D.ps
@@ -18,7 +18,7 @@ WuGdict17a begin
 % /setboxbackgroundcolor { 0.93 setgray } def % default is white
 % setPalatino
 
-{ 8 aCol5 iColA } /ipCol x bind def % number of colours and colour style
+{ 7 aCol5 iColA } /ipCol x bind def % number of colours and colour style
 /pStyles [
    { 11 0 0 1. 1. pset  0 ipCol }
    { 12 1 0 1. 1. pset  1 ipCol }
diff --git a/pub/share/setup2D.ps b/pub/share/setup2D.ps
index 28d0d58d..569eabc7 100644
--- a/pub/share/setup2D.ps
+++ b/pub/share/setup2D.ps
@@ -18,7 +18,7 @@ WuGdict17a begin
 % /setboxbackgroundcolor { 0.93 setgray } def % default is white
 % setPalatino
 
-{ 10 aCol6 iColA } /icCol x bind def % number of colours and colour style
+{ 10 aCol6 {black} {white} iColAA } /icCol x bind def % number of colours and colour style
 
 /hxlow 11.3 def
 /hxhig 12.3 def
diff --git a/pub/share/wups17a.ps b/pub/share/wups17a.ps
index 6684cdc7..ac134036 100644
--- a/pub/share/wups17a.ps
+++ b/pub/share/wups17a.ps
@@ -157,7 +157,6 @@ WuGdict17a begin
    /wxmax x def
    /wxmin x def
    /wxlog x 0 eq not def
-   % prepare conversion world coord -> frame coord
    /wxdel wxmax wxmin wxlog { div log } { sub } ifelse def
    /wxd { % dx(world) | dx(frame)
       wxlog { log } if wxdel div 10 mul
@@ -180,19 +179,29 @@ WuGdict17a begin
       wyd
       } bind def
    } def
-/hSetCoord { % log min max | -
+/hSetCoord { % log min max | - % for use in 2D plot
    /whmax x def
    /whmin x def
    /whlog x 0 eq not def
-   % prepare conversion world coord -> frame coord
-   /whdel whmax whmin whlog { div log } { sub } ifelse def
-   /whd { % dx(world) | dx(frame)
-      whlog { log } if whdel div 10 mul
-      } bind def
-   /wh { % x(world) | x(frame)
-      whmin whlog { div } { sub } ifelse
-      whd
+   whlog {
+      /whdel whmax whmin div log def
+      /wh { % h(world) | h(col)
+         dup whmin lt {
+            pop -1
+         } {
+            dup whmax gt {
+	       pop 11
+	    } {
+	       whmin div log whdel div 10 mul
+	    } ifelse
+         } ifelse
       } bind def
+   } {
+      /whdel whmax whmin sub def
+      /wh { % h(world) | h(col)
+         whmin sub whdel div 10 mul
+         } bind def
+      } ifelse
    } def
 
 % pair conversion
@@ -218,8 +227,8 @@ WuGdict17a begin
    2 -1 roll weightA mul x weightB mul add 3 1 roll
    } def
 
-/relcol { % i_col n_col | rel(0..1) : for one-dimensional choices
-   1 sub div 0 max 1 min
+/relcol { % i_col max_i | rel(0..1) % where i_col=0,..,max_i
+   div 0 max 1 min
    } def
 
 %%  Named colors:
@@ -283,6 +292,23 @@ WuGdict17a begin
    4 3 roll aCol x get exec colormix setRGBcolor
    } def
 
+%% ditto with discontinuous values for h<0 and h>1
+
+/iColAA { % i i_max arr {low_col} {hig_col} | -
+   5 3 roll
+   div
+   dup 0 lt {
+      pop pop x pop exec
+   } {
+      dup 1 gt {
+         pop 3 1 roll pop pop exec
+      } {
+         1 % arr {l} {h} i i_max
+	 5 2 roll pop pop
+	 iColA
+         } ifelse
+      } ifelse
+   } def
 
 %% Color arrays for non-linear one-dimensional choices:
 
@@ -302,29 +328,29 @@ WuGdict17a begin
    {   0   0 255 } % 13
    ] def
 /aCol2 [ % orange-red-blue-darkblue
-   { 255 180   0 } %  1
-   { 255 160   0 } %  1
-   { 255 120   0 } %  2
-   { 255  70   0 } %  3
-   { 255   0   0 } %  4
-   { 220  30  30 } %  5
-   { 220  70  60 } %  6
-   { 220 100 110 } %  7
-   { 200 130 130 } %  8
-   { 200 130 160 } %  9
-   { 180 110 180 } % 10
-   { 165 110 185 } % 11
-   { 150 130 190 } % 12
-   { 130 150 210 } % 13
-   { 100 120 220 } % 14
-   {  85 105 230 } % 15
-   {  70  90 255 } % 16
-   {   0   0 255 } % 17
-   {   0   0 180 } % 18
-   {  10  10 150 } % 19
-   {  30  30 130 } % 20
+   { 255 180   0 }
+   { 255 160   0 }
+   { 255 120   0 }
+   { 255  70   0 }
+   { 255   0   0 }
+   { 220  30  30 }
+   { 220  70  60 }
+   { 220 100 110 }
+   { 200 130 130 }
+   { 200 130 160 }
+   { 180 110 180 }
+   { 165 110 185 }
+   { 150 130 190 }
+   { 130 150 210 }
+   { 100 120 220 }
+   {  85 105 230 }
+   {  70  90 255 }
+   {   0   0 255 }
+   {   0   0 180 }
+   {  10  10 150 }
+   {  30  30 130 }
    ] def
-/aCol3 [ % [fixed size: 9] siemenscolors
+/aCol3 [ % [fixed size: max_i=8] siemenscolors
    { 165   0  33 } % siemensred
    {  33 153 102 } % siemensgreen
    { 0   102 153 } % siemensblue
@@ -358,7 +384,7 @@ WuGdict17a begin
    { 120  80  30 }
    { 100  60  20 }
    ] def
-/aCol5 [ % [fixed size: 8] old gnuplot default (see man gnuplot and rgb.txt)
+/aCol5 [ % [fixed size: max_i=7] old gnuplot default (see man gnuplot and rgb.txt)
    { 255   0   0 } % red
    {   0 255   0 } % green
    {   0   0 255 } % blue
@@ -1120,10 +1146,10 @@ WuGdict17a begin
       } {
          1 eq {
             fill
-	 } {
-	    gsave fill grestore
-	    gsave black st grestore
-	 } ifelse
+         } {
+            gsave fill grestore
+            gsave black st grestore
+         } ifelse
       } ifelse
    } def
 
@@ -1249,10 +1275,11 @@ WuGdict17a begin
 } def
 
 /ColorLegend { % - | - [plot color legend to rectangle hxlow hxhig 0 10]
-   /eps .05 def
-   eps eps 10.001 {
-      dup icCol
-      dup eps sub x
+   /ncol 100 def
+   1 1 ncol {
+      dup ncol div 10 mul icCol
+      dup ncol div 10 mul
+      x 1 sub ncol div 10 mul
       hxlow hxhig 4 2 roll rect fill
       } for
    } def
diff --git a/pub/trivia/vector_ops.cpp b/pub/trivia/vector_ops.cpp
index 36f30708..49a0d6b6 100644
--- a/pub/trivia/vector_ops.cpp
+++ b/pub/trivia/vector_ops.cpp
@@ -167,3 +167,17 @@ std::string triv::indices_to_s(const vector<int>& v)
     }
     return ret;
 }
+
+//! Determines minimal positive value and maximal value in vector v.
+//! Results must be tested for finiteness.
+void triv::getPosminMax(const std::vector<double>& v, double& minval, double&maxval)
+{
+    minval = +INFINITY;
+    maxval = -INFINITY;
+    for (auto t = v.begin(); t<v.end(); ++t) {
+        if (*t<=0)
+            continue;
+        minval = std::min(minval, *t);
+        maxval = std::max(maxval, *t);
+    }
+}
diff --git a/pub/trivia/vector_ops.hpp b/pub/trivia/vector_ops.hpp
index 22edd3e4..b47e28bb 100644
--- a/pub/trivia/vector_ops.hpp
+++ b/pub/trivia/vector_ops.hpp
@@ -8,6 +8,7 @@
 #ifndef VECTOR_OPS_H
 #define VECTOR_OPS_H
 
+#include <cmath>
 #include <vector>
 #include <functional>
 #include <algorithm>
@@ -25,6 +26,7 @@ namespace triv
     std::string indices_to_s(const std::vector<int>& v);
 
     template <class T> void getMinMax(const std::vector<T>& v, T& minval, T&maxval);
+    void getPosminMax(const std::vector<double>& v, double& minval, double&maxval);
     template <class T> bool contains(const std::vector<T>& v, const T e);
     template <class T, class Pred> std::vector<T> merge_sorted(
         const std::vector<T>& a, const std::vector<T>& b, Pred a_before_b);
-- 
GitLab