histogram.rb in gphys-1.4.3

- old
+ new

@@ -32,29 +32,37 @@
       else
         nbins = opts["nbins"] || gphys0.total/500
         nbins = 10 if nbins < 10
         min = opts["min"] || gphys0.min.val
         max = opts["max"] || gphys0.max.val
+        if log_bins = (opts["log_bins"] && (min > 0))
+          min = Math.log10(min)
+          max = Math.log10(max)
+        end
         hist = GSL::Histogram.alloc(nbins,[min,max])
       end
       val = gphys0.val
       val = val.get_array![val.get_mask!] if NArrayMiss === val
+      val = NMath.log10(val) if log_bins
       hist.increment(val)
 
       bounds = hist.range.to_na
+      bounds = 10 ** bounds if log_bins
       center = (bounds[0..-2]+bounds[1..-1])/2
+      cell_width = (bounds[1..-1] - bounds[0..-2]) / 2
       name = gphys0.name
       attr = gphys0.data.attr_copy
       bounds = VArray.new(bounds, attr, name)
       center = VArray.new(center, attr, name)
       axis = Axis.new(true)
       axis.set_cell(center, bounds, name)
       axis.set_pos_to_center
 
       bin = hist.bin.to_na
+      bin /= cell_width if opts["log_bins"]
       bin = VArray.new(bin,
-                       {"long_name"=>"number in bins", "units"=>"1"},
+                       {"long_name" => (log_bins ? "number per unit bin width" : "number in bins"), "units"=>"1"},
                        "bin")
       new_gphys = GPhys.new(Grid.new(axis), bin)
       new_gphys.set_att("mean",[hist.mean])
       new_gphys.set_att("standard_deviation",[hist.sigma])
       return new_gphys
@@ -155,11 +163,12 @@
 
     @@histogram_options = Misc::KeywordOptAutoHelp.new(
             ['window', [nil,nil,nil,nil], "window bounds"],
             ['title', "histogram", "window title"],
             ['exchange', false, "exchange x and y"],
-            ['fill', false, "fill bars"]
+            ['fill', false, "fill bars"],
+            ['fill_pattern', nil, "fill pattern"]
                                                        )
     def histogram(gphys, newframe=true, options=nil)
       gropn_1_if_not_yet
       if newframe!=true && newframe!=false
         raise ArgumentError, "2nd arg (newframe) must be true of false"
@@ -171,11 +180,14 @@
         raise ArgumentError, "axis must be cell type"
       end
       # if window is specified via GGraph#fig or GGraph#set_fig, use it.
       opts = @@histogram_options.interpret(options)
       exchange = opts["exchange"]
-      window = opts["window"] || @@fig['window']
+      raise "Option 'window' must be an Array of length == 4" unless opts["window"].is_a?(Array) && opts["window"].length == 4
+      window = opts["window"].dup
+      4.times{|i| window[i] ||= @@fig['window'][i]} if @@fig['window']
+      4.times{|i| window[i] ||= DCL::sgqwnd[i]} unless newframe
       unless exchange
         x = gphys.axis(0).cell_bounds
         y = gphys
         window[2] ||= 0
       else
@@ -199,11 +211,13 @@
           tmp.set_mask(mask * tmp.get_mask)
         else
           tmp = NArrayMiss.to_nam_no_dup(tmp, mask)
         end
         y.replace_val(tmp)
-        window[2] = tmp.abs.min * (tmp.min < 0 ? -1 : 1)
+        if(window[2].nil? || window[2] == 0)
+          window[2] = tmp.abs.min * (tmp.min < 0 ? -1 : 1)
+        end
       end
       if (itr==3 || itr==4)
         tmp = x.val
         if tmp.min * tmp.max < 0
           if tmp.min.abs < tmp.max
@@ -218,42 +232,39 @@
           tmp.set_mask(mask * tmp.get_mask)
         else
           tmp = NArrayMiss.to_nam_no_dup(tmp, mask)
         end
         x.replace_val(tmp)
-        window[0] = tmp.abs.min * (tmp.min < 0 ? -1 : 1)
+        if(window[0].nil? || window[0] == 0)
+          window[0] = tmp.abs.min * (tmp.min < 0 ? -1 : 1)
+        end
       end
       opts["window"] = window
 
-      if newframe
-        fig(x, y, "window"=>opts["window"])
-        axes(x, y, "title"=>opts["title"])
-      end
+      fig(x, y, "window" => window) if newframe
       lmiss = DCL::gllget("lmiss")
       DCL::gllset("lmiss", true)
+      if opts["fill"]
+        itps = DCL::uuqarp
+        itps[0] = itps[1] = opts["fill_pattern"] if opts["fill_pattern"]
+        DCL::uusarp(*itps)
+      end
       unless exchange
         if opts["fill"]
           DCL::uvbxa(x.val, [window[2]] * y.length, y.val)
         end
-        if (itr == 2) || (itr == 4)
-          bottom = [y.val.min] * y.length
-        else
-          bottom = [window[2]] * y.length
-        end
+        bottom = [window[2]] * y.length
         DCL::uvbxf(x.val, bottom, y.val)
       else
         if opts["fill"]
           DCL::uhbxa([window[0]] * x.length, x.val, y.val)
         end
-        if (itr == 3) || (itr == 4)
-          bottom = [x.val.min] * x.length
-        else
-          bottom = [window[0]] * x.length
-        end
+        bottom = [window[0]] * x.length
         DCL::uhbxf(bottom, x.val, y.val)
       end
       DCL::gllset("lmiss", lmiss)
+      axes(x, y, "title" => opts["title"]) if newframe
       return nil
     end
     alias :histogram1D :histogram
 
   end
@@ -307,11 +318,17 @@
 
   GGraph.set_fig("itr"=>4)
   GGraph.histogram(hist)
 
   GGraph.set_fig("itr"=>1)
+  GGraph.histogram(hist, true, "fill" => true, "fill_pattern" => 15999)
+
+  GGraph.set_fig("itr"=>1)
   hist = gphys1D.histogram("nbins"=>10)
   GGraph.histogram(hist, true, "title"=>"histogram 1D")
+
+  hist = GAnalysis.histogram(gphys1D, 'log_bins' => true)
+  GGraph.histogram(hist)
 
   hist = GAnalysis.histogram2D(gphys2D_0, gphys2D_1, "nbins0"=>50)
   GGraph.tone(hist, true, "tonc"=>true)
 
 