;----------------------------------------------------------------------
; This tests NCL-2161, a possible bug in bin_sum
;----------------------------------------------------------------------

APPLY_THE_FIX = True    ; Whether to attempt to fix the NCL version of "bin_sum"

;----------------------------------------------------------------------
; This calls the built-in version of bin_sum and returns
; both gcnt and gbin.
;----------------------------------------------------------------------
function bin_sum_builtin(x, y, z, binMin, binMax, binWidth )
local nBin, glat, glon, gbin
begin
  nBin = 1 + round( ( binMax - binMin ) / binWidth, 3 )
  glat = fspan( binMin, binMax, nBin )
  glon = fspan( binMin, binMax, nBin )

;---Do the binning using "bin_sum"
  gcnt = new( (/ nBin, nBin /), integer )
  gbin = new( (/ nBin, nBin /), float )
  gcnt = 0
  gbin = 0
  bin_sum( gbin, gcnt, glon, glat, x, y, z )

  gcnt!0 = "y"            ;---For plotting purposes
  gcnt!1 = "x"
  gcnt&y = glat
  gcnt&x = glon

  return([/gbin,gcnt/])
end

;----------------------------------------------------------------------
; This is an NCL version of the Fortran "bin_sum" function.
;
; It should produce identical (i.e. incorrect) results as bin_sum if
; APPLY_THE_FIX is set to False. Otherwise, it will produce 
; identical results as bin_sum_other (see below).
;
; Note: I believe the bug is in the code that calculates "nl"
; and "ml" below, which is also in the Fortran version of this
; routine. Note how it is binning values that fall outside the
; range of glat/glon.
;----------------------------------------------------------------------
function bin_sum_ncl(x, y, z, binMin, binMax, binWidth )
local binWidth, nBin, glat, glon, gbin, dx, dy, fmt, \
      dlat, dlon, glatbnd, blonbnd, k, nl, ml, delta
begin
  fmt = "%9.2f"

  nBin     = 1 + round( ( binMax - binMin ) / binWidth, 3 )

  glat = fspan( binMin, binMax, nBin )
  glon = fspan( binMin, binMax, nBin )

  min_lat = min(y)
  max_lat = max(y)
  min_lon = min(x)
  max_lon = max(x)

  min_glat = min(glat)
  max_glat = max(glat)
  min_glon = min(glon)
  max_glon = max(glon)

  org_range = "Original range: [" + sprintf(fmt,min_lat) + "," +  sprintf(fmt,max_lat) + "] [" + \
                                    sprintf(fmt,min_lon) + "," +  sprintf(fmt,max_lon) + "]"
  new_range = "New range:      [" + sprintf(fmt,min_glat) + "," +  sprintf(fmt,max_glat) + "] [" + \
                                    sprintf(fmt,min_glon) + "," +  sprintf(fmt,max_glon) + "]"

  dx = (glon(1)-glon(0))/2.
  dy = (glat(1)-glat(0))/2.

  dlat = abs(glat(1)-glat(0))
  dlon = abs(glon(1)-glon(0))

  glatbnd = glat(0) - dlat/2.
  glonbnd = glon(0) - dlon/2.

  delta = (glat(1)-glat(0))/2.

  gbin = new((/ nBin, nBin /), typeof(z))
  gcnt = new((/ nBin, nBin /), integer)
  gbin = 0
  gcnt = 0
  
  do k = 0, dimsizes(z)-1
    nl = toint(abs((y(k)-glatbnd)/dlat))
    ml = toint(abs((x(k)-glonbnd)/dlon))
    if(ml.ge.0.and.ml.lt.nBin.and.nl.ge.0.and.nl.lt.nBin) then
      glat1 = glat(nl)-delta
      glat2 = glat(nl)+delta
      glon1 = glon(ml)-delta
      glon2 = glon(ml)+delta

      print("============================================================================")
      print(""+org_range)
      print(""+new_range)
      print("Value            (" + sprintf(fmt,y(k)) + ", " + sprintf(fmt,x(k)) +")")

      if(APPLY_THE_FIX.and.\
         (y(k).lt.glat1.or.y(k).gt.glat2.or.\
          x(k).lt.glon1.or.x(k).gt.glon2)) then
          print("OUT-OF-RANGE VALUE FOUND: fix being applied: skipping...")
          continue
      else
        gcnt(nl,ml) = gcnt(nl,ml) + 1
        gbin(nl,ml) = gbin(nl,ml) + z(k)
      end if

      print("got binned into: [" + sprintf(fmt,glat1) + "," + \
                                   sprintf(fmt,glat2) + "] [" + \
                                   sprintf(fmt,glon1) + "," + \
                                   sprintf(fmt,glon2) + "]")
      print("(ilt,iln)        (" + nl + "," + ml + ")")

      if(y(k).lt.glat1.or.y(k).gt.glat2) then
        print("ERROR: THIS FALLS OUTSIDE THE LAT RANGE")
      end if
      if(x(k).lt.glon1.or.x(k).gt.glon2) then
        print("ERROR: THIS FALLS OUTSIDE THE LON RANGE")
      end if
    end if
  end do

  gcnt!0 = "y"            ;---For plotting purposes
  gcnt!1 = "x"
  gcnt&y = glat
  gcnt&x = glon

  return([/gbin,gcnt/])
end

;----------------------------------------------------------------------
; This is a different binning function, which
; uses the actual bounds of the lat/lon data for
; the binning.
;----------------------------------------------------------------------
function bin_sum_other(x, y, z, binMin, binMax, binWidth )
local binWidth, nBin, glat, glon, gbin, dx, dy, \
      dlat, dlon, glatbnd, blonbnd, k, ilt, iln, delta 
begin
  nBin = 1 + round( ( binMax - binMin ) / binWidth, 3 )

  glat = fspan( binMin, binMax, nBin )
  glon = fspan( binMin, binMax, nBin )

  dx = (glon(1)-glon(0))/2.
  dy = (glat(1)-glat(0))/2.

  dlat = abs(glat(1)-glat(0))
  dlon = abs(glon(1)-glon(0))

  glatbnd = glat(0) - dlat/2.
  glonbnd = glon(0) - dlon/2.

  delta = (glat(1)-glat(0))/2.

  gbin = new( (/ nBin, nBin /), float )
  gcnt = new((/ nBin, nBin /), integer)
  gcnt = 0
  gbin = 0
  do k = 0,dimsizes(z)-1
    ilt = ind((glat-delta).le.y(k).and.y(k).lt.(glat+delta))
    iln = ind((glon-delta).le.x(k).and.x(k).lt.(glon+delta))
    if (.not.ismissing(ilt).and..not.ismissing(iln).and.\
        ilt.ge.0 .and. ilt.lt.nBin.and.iln.ge.0 .and.iln.lt.nBin) then
;;      print("other bin method: point " + y(k) + "/" + x(k) + " " + \
;;            "got binned into lat: " + (glat(ilt)-delta) + " to " + \
;;                                      (glat(ilt)+delta) + ", lon: " + \
;;                                      (glon(iln)-delta) + " to " + \
;;                                      (glon(iln)+delta) + \
;;                                     " (" + ilt + "," + iln + ")")
      gcnt(ilt,iln) = gcnt(ilt,iln)+1
      gbin(ilt,iln) = gbin(ilt,iln)+z(k)
    else
;---Make sure we get any data equal to upper bounds of lat/lon
      ilt = ind(y(k).eq.(glat+delta))
      iln = ind(x(k).eq.(glon+delta))
      if (.not.ismissing(ilt).and..not.ismissing(iln).and.\
          ilt.ge.0 .and. ilt.lt.nBin.and.iln.ge.0 .and.iln.lt.nBin) then
;;        print("other bin method: point " + y(k) + "/" + x(k) + " " + \
;;              "got binned into lat: " + (glat(ilt)-delta) + " to " + \
;;                                        (glat(ilt)+delta) + ", lon: " + \
;;                                        (glon(iln)-delta) + " to " + \
;;                                        (glon(iln)+delta))
        gcnt(ilt,iln) = gcnt(ilt,iln)+1
        gbin(ilt,iln) = gbin(ilt,iln)+z(k)
      end if
    end if
  end do

  gcnt!0 = "y"            ;---For plotting purposes
  gcnt!1 = "x"
  gcnt&y = glat
  gcnt&x = glon

  return([/gbin,gcnt/])
end


;----------------------------------------------------------------------
; Create contours of given count array
;----------------------------------------------------------------------
function plot_gcnt(wks, gcnt, title, binMin, binMax)
local res
begin
  res                      = True

  res@gsnDraw              = False
  res@gsnFrame             = False

  res@trXMinF              = -500
  res@trXMaxF              =  500
  res@trYMinF              = -500
  res@trYMaxF              =  500

  res@cnLevelSelectionMode = "ExplicitLevels"
  res@cnLevels             = ispan( 0, 6, 1 )

  res@cnFillOn             = True
  res@cnFillMode           = "CellFill"

  res@cnLinesOn            = False
  res@cnLineLabelsOn       = False

  res@cnInfoLabelOn        = False
  res@lbLabelBarOn         = False

  res@gsnLeftString        = ""+binMin
  res@gsnRightString       = ""+binMax
  res@tiMainString         = title

  plot = gsn_csm_contour(wks, gcnt, res)

  return(plot)
end

;----------------------------------------------------------------------
; Panel three plots
;----------------------------------------------------------------------
procedure plot_panel(wks, plotb, plotn, ploto)
local pres
begin
  pres                  = True
  pres@gsnMaximize      = True
  pres@gsnPanelLabelBar = True

  gsn_panel(wks,(/plotb, plotn, ploto/),(/1,3/),pres)

end

;----------------------------------------------------------------------
; Main code
;----------------------------------------------------------------------
begin

  wks = gsn_open_wks( "x11","bin")

  nData = 1000
  xData = random_normal( 0, 1000, nData )
  yData = random_normal( 0, 1000, nData )
  zData = new( nData, float )
  zData = 1. / nData

  printMinMax(xData,0)
  printMinMax(yData,0)
  printMinMax(zData,0)

;
; These are the various binning ranges to try. The
; rectilinear lat/lon arrays will be created based
; on these values.
;
; Note that bin_sum only seems to work for the
; -2900/3100 range, which is close to the original
; range of xData/yData.
;
  binMin = (/-2900,-500,-500,-600,-400,-400,-300/)
  binMax = (/ 3100, 500, 600, 500, 400, 300, 300/)
  binWdt = 50

  do nb=0,dimsizes(binMin)-1
    blist := bin_sum_builtin(xData, yData, zData, binMin(nb), binMax(nb), binWdt )
    gbin_blt := blist[01]
    gcnt_blt := blist[1]

    blist := bin_sum_ncl(xData, yData, zData, binMin(nb), binMax(nb), binWdt )
    gbin_ncl := blist[01]
    gcnt_ncl := blist[1]

    blist := bin_sum_other(xData, yData, zData, binMin(nb), binMax(nb), binWdt )
    gbin_other := blist[01]
    gcnt_other := blist[1]

    print("============================================================================")
    print("Apply the fix? " + APPLY_THE_FIX)
    print("Is gbin_blt equal to gbin_ncl?   " + all(gbin_blt.eq.gbin_ncl))
    print("Is gcnt_blt equal to gcnt_ncl?   " + all(gcnt_blt.eq.gcnt_ncl))
    print("Is gbin_ncl equal to gbin_other? " + all(gbin_ncl.eq.gbin_other))
    print("Is gcnt_ncl equal to gcnt_other? " + all(gcnt_ncl.eq.gcnt_other))
    print("Total # of values binned in gcnt_blt   = " + sum(gcnt_blt))
    print("Total # of values binned in gcnt_ncl   = " + sum(gcnt_ncl))
    print("Total # of values binned in gcnt_other = " + sum(gcnt_other))

    plot_blt   = plot_gcnt( wks, gcnt_blt,   "Built-in bin_sum", binMin(nb), binMax(nb))
    plot_ncl   = plot_gcnt( wks, gcnt_ncl,   "NCL version of  bin_sum", binMin(nb), binMax(nb))
    plot_other = plot_gcnt( wks, gcnt_other, "Different binning method", binMin(nb), binMax(nb))

    plot_panel(wks, plot_blt, plot_ncl, plot_other)
  end do  
end