;***********************************************************************;
; Function : hist_columns                                               ;
;                   wks: workstation object                             ;
;                    xy: graphic                                        ;
;             binvalues: numeric                                        ;
;               barlocs: numeric                                        ;
;              barwidth: numeric                                        ;
;                colors                                                 ;
;               compare: logical                                        ;
;                 gsres: logical                                        ;
;                                                                       ;
;  xy        - xy plot id to draw columns on                            ;
;  bins      - the center of each bin range.                            ;
;  binvalues - the number of values in the corresponding bin            ;
;  barlocs   - the start of the first bar in each bin                   ;
;  width     - the width of the bar                                     ;
;  colors    - array of colors to use                                   ;
;  gsres     - optional primitive resources                             ;
;                                                                       ;
; This function creates the columns for a histogram plot. The Y axis    ;
; will represent the number of values in a bin and a percentage.        ;
;                                                                       ;
;***********************************************************************;
undef("hist_columns")
function hist_columns(wks[1]:graphic,xy[1]:graphic,binvalues:numeric, \
                      barlocs[*][*]:numeric, barwidth[*]:numeric, \
                      colors, compare:logical, gsres:logical)
local i, nbins, dims, nbinvalues, gsres, xpoints, ypoints, multibars,\
col_dims, col_rank
begin
  dims       = dimsizes(barlocs)
  nbarsinbin = dims(0)
  nbars      = dims(1)
  delete(dims)

  dims = dimsizes(binvalues)
  if(dimsizes(dims).eq.1) then
    nbinvalues = dims(0)
  else
    nbinvalues = dims(1)
  end if
  delete(dims)
  if(nbars.ne.nbinvalues) then
    print("Error: hist_columns: Dimension sizes of bins (" + nbars+ ") and binvalues (" + nbinvalues + ") must be the same")
    return
  end if

  if(nbarsinbin.ge.2)
    multibars   = True
  else
    multibars   = False
  end if

;
; Set up arrays to hold polygon points.
;
  if(multibars) then
    xpoints   = new((/nbarsinbin,nbars,5/),float)
    ypoints   = new((/nbarsinbin,nbars,5/),float)
    polygons  = new((/nbarsinbin,nbars/),graphic)
  else
    xpoints  = new((/nbars,5/),float)
    ypoints  = new((/nbars,5/),float)
    polygons = new(nbars,graphic)
  end if
;
; Set up variable to hold resources.
;
  gsres           = True
  set_attr(gsres,"gsEdgesOn",True)
;
; Begin assigning polygon points.
;
  if(multibars)
    do i=0,nbarsinbin-1
      ypoints(i,:,0) = (/0/)
      ypoints(i,:,1) = (/binvalues(i,:)/)
      ypoints(i,:,2) = (/ypoints(i,:,1)/)
      ypoints(i,:,3) = (/0/)
      ypoints(i,:,4) = (/0/)

      xpoints(i,:,0) = (/barlocs(i,:)/)
      xpoints(i,:,1) = (/xpoints(i,:,0)/)
      xpoints(i,:,2) = (/barlocs(i,:) + barwidth/)
      xpoints(i,:,3) = (/xpoints(i,:,2)/)
      xpoints(i,:,4) = (/xpoints(i,:,0)/)
    end do
  else
    ypoints(:,0) = (/0/)
    ypoints(:,1) = (/binvalues/)
    ypoints(:,2) = (/ypoints(:,1)/)
    ypoints(:,3) = (/0/)
    ypoints(:,4) = (/0/)

    xpoints(:,0) = (/barlocs(0,:)/)
    xpoints(:,1) = (/xpoints(:,0)/)
    xpoints(:,2) = (/barlocs(0,:) + barwidth/)
    xpoints(:,3) = (/xpoints(:,2)/)
    xpoints(:,4) = (/xpoints(:,0)/)
  end if

  if(compare)
    fillindex = get_res_value(gsres,"gsFillIndex",(/0,6/))
  else
    if(multibars) then
      fillindex = get_res_value(gsres,"gsFillIndex",0)
      if(dimsizes(dimsizes(fillindex)).eq.1) then
         ftmp = new(nbarsinbin,typeof(fillindex))
         ftmp = fillindex
         delete(fillindex)
         fillindex = ftmp
         delete(ftmp)
      end if
    else
      fillindex = get_res_value(gsres,"gsFillIndex",0)
    end if
  end if
;
; Make sure fill indices are between 0 and 17.
; No more, because you have have added your own fill index.
;
;  fillindex = min((/max((/max(fillindex),0/)),17/))

  rgba_arr = convert_color_to_rgba(wks,colors)
  col_dims = dimsizes(rgba_arr)   ; better be 2D or 3D
  col_rank = dimsizes(col_dims)
  if(col_rank.gt.3.or.(col_rank.eq.3.and.col_dims(0).ne.nbarsinbin)) then
    print("Error: hist_columns: Invalid dimensionality of array of colors passed in.")
    print("   Was expecting a 1D or 2D array of colors")
    exit
  end if
  if(col_rank.eq.3) then
    ncolors = col_dims(1)
  else
    ncolors = col_dims(0)
  end if
  do i = 0, nbars - 1
    if(col_rank.eq.1) then
      gsres@gsFillColor = rgba_arr
    else if(col_rank.eq.2) then
      gsres@gsFillColor = rgba_arr(i % ncolors,:)
    end if
    end if
    if(multibars) 
;
; Add the remaining histogram bars first, so that they will be drawn
; second. The first histogram bars will thus be drawn on top, if
; they are being stacked.
;
      do j=nbarsinbin-1,0,1
        if(col_rank.eq.3) then
          gsres@gsFillColor = rgba_arr(j,i % ncolors,:)
        end if
        gsres@gsFillIndex = fillindex(j)
        if(binvalues@horizontal) then
          polygons(j,i) = gsn_add_polygon(wks,xy,ypoints(j,i,:), \
                                         xpoints(j,i,:), gsres)
        else
          polygons(j,i) = gsn_add_polygon(wks,xy,xpoints(j,i,:), \
                                          ypoints(j,i,:), gsres)
        end if
      end do
    else
      gsres@gsFillIndex = fillindex
      if(binvalues@horizontal) then
        polygons(i) = gsn_add_polygon(wks,xy,ypoints(i,:),xpoints(i,:),gsres)
      else
        polygons(i) = gsn_add_polygon(wks,xy,xpoints(i,:),ypoints(i,:),gsres)
      end if
    end if
  end do
;
; Return the polygons created as attributes of the XY plot. This is
; necessary, b/c otherwise the polygons will go away when you exit the
; function.
;
  var_string      = unique_string("hpolygons")
  xy@$var_string$ = polygons
  return(xy)
end

;***********************************************************************;
; Function : compute_hist_vals                                          ;
;                     x: numeric or string                              ;
;               binlocs: numeric or string                              ;
;              nbinlocs: integer                                        ;
;             bin_width: numeric                                        ;
;           setinterval: logical                                        ;
;           setdiscrete: integer                                        ;
;            minmaxbins: logical                                        ;
;             count_msg: logical                                        ;
;                isnice: integer                                        ;
;               compare: integer                                        ;
;            sidebyside: integer                                        ;
;                                                                       ;
; By default, this routine calculates a nice set of ranges for "binning";
; the given data. The following cases are possible:                     ;
;                                                                       ;
; 1. If setinterval is True, then the user has set their own bin        ;
;    intervals via either the gsnHistogramBinIntervals or the           ;
;    gsnHistogramClassIntervals resource (stored in binlocs array).     ;
;                                                                       ;
; 2. If setdiscrete is True, then the user has set discrete bin values  ;
;    via the gsnHistogramDiscreteBinValues resource (stored in          ;
;    "binlocs" array).                                                  ;
;                                                                       ;
; 3. If neither setinterval or setdiscrete is True, and if the resource ;
;    gsnHistogramBinWidth is set, then its value will be used as a bin  ;
;    width (bin_width). By default, bin_width will only be used as an   ;
;    approximate value, because it attempts to select "nice" values     ;
;    based on a width close to bin_width. If the resource               ;
;    gsnHistogramSelectNiceIntervals (isnice) is set to False, then     ;
;    you will get a bin width exactly equal to bin_width.               ;
;                                                                       ;
; 4. If neither setinterval or setdiscrete is True, and if the resource ;
;    gsnHistogramNumberOfBins is set, then its value (nbinlocs) will be ;
;    used to determine the number of bins. By default, nbinlocs is only ;
;    used as an approximate value, because it attempts to select "nice" ;
;    values based on a number of bins close to nbinlocs. If the         ;
;    resource gsnHistogramSelectNiceIntervals (isnice) is set to False, ;
;    then you will get a number of bins exactly equal to nbinlocs.      ;
;                                                                       ;
; 5. If no special resources are set, then this routine defaults to     ;
;    calculating approximately 10 "nice" bin intervals, based on the    ;
;    range of the data.                                                 ;
;                                                                       ;
;***********************************************************************;
undef("compute_hist_vals")
function compute_hist_vals(xx,binlocs,nbinlocs[1]:integer, \
                           bin_width:numeric,setinterval:logical, \
                           setdiscrete:logical, minmaxbins:logical, \
                           count_msg:logical,isnice:logical, \
                           compare:logical,sidebyside:logical)
local xmin, xmax, new_binlocs, nbars, buckets, x, i
begin
  if(isnumeric(xx)) then
    xmin = tofloat(min(xx))
    xmax = tofloat(max(xx))
  end if

;
; If the bins are set, need to determine if you want to have the bins
; represent ranges of values or discrete values.
;
  if(setdiscrete) then
    new_binlocs = binlocs
    nbars       = dimsizes(new_binlocs)   ; # of bars equals # of binlocs
    nsets       = dimsizes(new_binlocs)   ; # of bars equals # of binlocs
  else
;
; Check if range values have been set by user, or if we need to
; calculate them.
;
    if(setinterval) then
      new_binlocs = binlocs
    else
      if(nbinlocs.lt.0) then
        if(bin_width.lt.0.) then
          nbinlocs = 10                   ; Default to 10 bin locations.
        else 
          nbinlocs = floattointeger(((xmax - xmin)/bin_width))
          if(nbinlocs.le.0) then
            print("Warning: compute_hist_vals: cannot use given bin width. Defaulting...")
            nbinlocs = 10 
          end if
        end if
      end if
      if(.not.setdiscrete) then
        if(isnice) then
;
; Based on min and max of data, compute a new min/max/step that will
; give us "nice" bin values.
;
           nicevals = nice_mnmxintvl(xmin,xmax,nbinlocs,True)
           nvals    = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
           new_binlocs = fspan(nicevals(0),nicevals(1),nvals)
         else
;
; Don't bother with "nice" values; just span the data. 
;
           new_binlocs = fspan(xmin,xmax,nbinlocs+1)
         end if
      end if
    end if
    nbars = dimsizes(new_binlocs)-1
  end if
;
; Count number of values in a particular bin range, or exactly
; equal to a bin value if discrete.
;
  if(compare.or.sidebyside) then
    dims  = dimsizes(xx)
    nsets = dims(0)
    npts  = dims(1)
    x     = xx
  else
    nsets  = 1
    npts   = dimsizes(xx)
    x      = new((/1,npts/),typeof(xx))
    x(0,:) = xx
  end if

;
; Set up variable to hold binned values. Binned values can have
; unequal spacing.
;
  num_in_bins = new((/nsets,nbars/),integer)
;
; Count the values in each discrete bin.
;
  if(setdiscrete) then
    do j = 0,nsets-1
      do i = 0, nbars-1
        num_in_bins(j,i) = num(x(j,:).eq.new_binlocs(i))
      end do
    end do
  else
;
; Count the values in each bin interval. Bin intervals can be
; of length 0, meaning an exact count is done.
;
    do j = 0,nsets-1
      do i = 0, nbars-1
        if(new_binlocs(i).eq.new_binlocs(i+1)) then
          num_in_bins(j,i) = num(x(j,:).eq.new_binlocs(i))
        else
;
; Special tests for last interval are required.
;
          if(i.eq.(nbars-1)) then
            if(nbars.gt.1.and.new_binlocs(i).eq.new_binlocs(i-1)) then
              num_in_bins(j,i) = num(x(j,:).gt.new_binlocs(i).and. \
                                     x(j,:).le.new_binlocs(i+1))
            else
              num_in_bins(j,i) = num(x(j,:).ge.new_binlocs(i).and. \
                                     x(j,:).le.new_binlocs(i+1))
            end if
          else
;
; If the previous interval was not really an interval, but an exact
; bin value, then be careful not to count those values in the current
; interval.
;
            if(i.gt.0.and.new_binlocs(i).eq.new_binlocs(i-1)) then
              num_in_bins(j,i) = num(x(j,:).gt.new_binlocs(i).and. \
                                     x(j,:).lt.new_binlocs(i+1))
            else
              num_in_bins(j,i) = num(x(j,:).ge.new_binlocs(i).and. \
                                     x(j,:).lt.new_binlocs(i+1))
            end if
          end if
        end if
      end do
    end do
  end if
;
; If minmaxbins is True, then we need to also count the values
; outside the range of our new_binlocs.
;
  if(minmaxbins) then
    new_num_in_bins            = new((/nsets,nbars+2/),integer)
    new_num_in_bins(:,1:nbars) = num_in_bins
    do j = 0,nsets-1
      new_num_in_bins(j,0)       = num(x(j,:).lt.new_binlocs(0))
      new_num_in_bins(j,nbars+1) = num(x(j,:).gt.new_binlocs(nbars))
    end do
    delete(num_in_bins)
    num_in_bins = new_num_in_bins
    delete(new_num_in_bins)
    nbars = nbars + 2
  end if
;
; Count number of missing values.
;
  num_missing = new(nsets,integer)
  if(isatt(x,"_FillValue")) then
    do j = 0,nsets-1
      num_missing(j) = num(ismissing(x(j,:)))
    end do
  else
    num_missing = 0
  end if
;
; If count_msg is True, then we need to bin the number of missing values.
;
  if(count_msg) then
    new_num_in_bins              = new((/nsets,nbars+1/),integer)
    new_num_in_bins(:,0:nbars-1) = num_in_bins
    new_num_in_bins(:,nbars)     = num_missing

    delete(num_in_bins)
    num_in_bins = new_num_in_bins
    delete(new_num_in_bins)
    nbars = nbars + 1
  end if

;
; Calculate percentages, both with and without missing points included.
;
  npts_nomiss = npts - num_missing

  if(compare.or.sidebyside) then
    percs          = (100.*num_in_bins)/tofloat(npts)
    percs_nomiss   = new((/nsets,nbars/),float)
    do i=0,nsets-1
      percs_nomiss(i,:) = (100.*num_in_bins(i,:))/tofloat(npts_nomiss(i))
    end do
  else
    percs        = (100.*num_in_bins(0,:))/tofloat(npts)
    percs_nomiss = (100.*num_in_bins(0,:))/tofloat(npts_nomiss)
  end if
;
; Return information.
; 
  num_in_bins@NumMissing     = num_missing
  num_in_bins@binlocs        = new_binlocs
  num_in_bins@percentages    = percs
  num_in_bins@percentages_nm = percs_nomiss

  delete(x)
  if(compare.or.sidebyside) then
    return(num_in_bins)
  else
    return(num_in_bins(0,:))
  end if
end