;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; test_bin_sum.ncl
; Carl Schreck (cjschrec@ncsu.edu)
; February 2015
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Description: Trying to replicate some odd behavior with bin_sum
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/diagnostics_cam.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl"
load "$CJS_NCL_LIB/lib.cjs_graphics.ncl"
load "$CJS_NCL_LIB/print_clock.ncl"


procedure manual_histogram( \
  io_wks[1] : graphic, \
  i_xData[*] : numeric, \
  i_yData[*] : numeric, \
  i_binMin[*] : numeric, \
  i_binMax[*] : numeric \
)
begin ; manual_histogram

  print_clock( "Manual " + i_binMin + " " + i_binMax )

  binWidth = 50
  nBin = 1 + round( ( i_binMax - i_binMin ) / binWidth, 3 )
  bins = fspan( i_binMin, i_binMax, nBin )

  count = new( (/ nBin, nBin /), integer )
  count = 0

  count!0 = "y"
  count!1 = "x"
  count&x = bins
  count&y = bins

  do xCounter = 0, nBin-1
    do yCounter = 0, nBin-1
      isInXbin = ( i_xData.ge.( bins(xCounter) - binWidth ) ) \
          .and.( i_xData.lt.( bins(xCounter) + binWidth ) )
      isInYbin = ( i_yData.ge.( bins(yCounter) - binWidth ) ) \
          .and.( i_yData.lt.( bins(yCounter) + binWidth ) )
      count(xCounter,yCounter) = num( isInXbin.and.isInYbin )
    end do
  end do

  res = True
  res@trXMinF = -500
  res@trXMaxF =  500
  res@trYMinF = -500
  res@trYMaxF =  500
  res@cnFillOn        = True
  res@cnFillMode      = "CellFill"
  res@cnLinesOn = False
  res@cnLevelSelectionMode = "ExplicitLevels"
  res@cnLevels = ispan( 0, 6, 1 )
  res@gsnLeftString = i_binMin
  res@gsnRightString = i_binMax
  res@gsnCenterString = "Manual"

  plot = gsn_csm_contour( io_wks, count, res )
end ; manual_histogram


procedure plot_histogram( \
  io_wks[1] : graphic, \
  i_xData[*] : numeric, \
  i_yData[*] : numeric, \
  i_binMin[*] : numeric, \
  i_binMax[*] : numeric \
)
begin ; plot_histogram

  print_clock( "bin_sum " + i_binMin + " " + i_binMax )

  binWidth = 50
  nBin = 1 + round( ( i_binMax - i_binMin ) / binWidth, 3 )
  bins = fspan( i_binMin, i_binMax, nBin )

  count = new( (/ nBin, nBin /), integer )
  count = 0

  count!0 = "y"
  count!1 = "x"
  count&x = bins
  count&y = bins

  percent = new( (/ nBin, nBin /), float )
  percent = 0
  copy_VarMeta( count, percent )

  z = new( dimsizes(i_xData), float )
  z = 1. / dimsizes(i_xData)

  bin_sum( percent, count, bins, bins, i_xData, i_yData, z )

  res = True
  res@trXMinF = -500
  res@trXMaxF =  500
  res@trYMinF = -500
  res@trYMaxF =  500
  res@cnFillOn        = True
  res@cnFillMode      = "CellFill"
  res@cnLinesOn = False
  res@cnLevelSelectionMode = "ExplicitLevels"
  res@cnLevels = ispan( 0, 6, 1 )
  res@gsnLeftString = i_binMin
  res@gsnRightString = i_binMax
  res@gsnCenterString = "bin_sum"

  plot = gsn_csm_contour( io_wks, count, res )
end ; plot_histogram


begin ; main
  print_clock( "Here we go!" )

  ; These are some parameters that could be useful to have up top

  plotType = "pdf"
  plotName = "figures/test_bin_sum"
  plotDpi  = 200

  wks = gsn_open_wks( plotType, plotName )

  xData = random_normal( 0, 1000, 1000 )
  yData = random_normal( 0, 1000, 1000 )

  binMin = -600  
  binMax =  600  
  manual_histogram( wks, xData, yData, binMin, binMax )

  binMin = -300  
  binMax =  300  
  manual_histogram( wks, xData, yData, binMin, binMax )

  binMin = -500  
  binMax =  500  
  plot_histogram( wks, xData, yData, binMin, binMax )

  binMin = -500  
  binMax =  600  
  plot_histogram( wks, xData, yData, binMin, binMax )

  binMin = -600  
  binMax =  500  
  plot_histogram( wks, xData, yData, binMin, binMax )

  binMin = -400  
  binMax =  400  
  plot_histogram( wks, xData, yData, binMin, binMax )

  binMin = -400  
  binMax =  300  
  plot_histogram( wks, xData, yData, binMin, binMax )

  binMin = -300  
  binMax =  300  
  plot_histogram( wks, xData, yData, binMin, binMax )

  print_clock( "Thank you, come again." )
end; main