;*************************************************
; histo_6.ncl
;
; Concepts illustrated:
;   - Comparing two sets of histograms
;   - Generating dummy data using "rand"
;   - Specifying the bar fill colors in a histogram
;
;************************************************
;
; This file is loaded by default in NCL V6.2.0 and newer
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   
;************************************************
 begin
  int_radar = ((/-9,-6,-3,0,3,6,9,12,15,18,21,24,27,30,33,36,39,42,45,48,51,54,57,60,63/))
  print("dimsizes(int_radar) = " + dimsizes(int_radar))
;************************************************
; Generate some random data.
;************************************************
; here it is
   num_o = (/0,0,13,160,2339,4632,7374,11044,14245,17172,24537,38696,51255,43006\
                    ,33918,32011,36194,46269,43924,26732,12078,4392,746,62,0/)
  num_m =        (/0,13873,13853,16611,14577,14316,14329,15375,17018,19546,24192,36264,\
         60084,66934,38879,83443,21433,6238,1702,499,56,9,14,0,0/)
; max_o = 3*max(num_o)
; max_m = 3*max(num_m)

;   wtype = "png"
  wtype = "pdf"
; for png --> gives "error" for pdf plot
  wtype@wkWidth  = 7500    ; produces a better looking PNG
  wtype@wkHeight = 7500
  dims_num_o = dimsizes(num_o)
  dims_num_m = dimsizes(num_m)
  obs_x = new(dims_num_o-1,float)
  mod_x = new(dims_num_m-1,float)
  ii_int = 0
; this does a loop over all array values
  obs_x = num_o(1:)
  mod_x = num_m(1:)
; note: the histogram plots data between intervals, so you must pad
; the data to the right by an increment of i_int
  dims_obs_x = dimsizes(obs_x)
  dims_mod_x = dimsizes(mod_x)
  print("dims_obs_x = " + dimsizes(obs_x))
  print("dims_mod_x = " + dimsizes(mod_x))
   
  sum_obs = floattointeger(sum(obs_x))
  sum_mod = floattointeger(sum(mod_x))
; check for 0 values
  do i_int = 0,dims_obs_x-1
  if (obs_x(i_int).eq.0)then
   sum_obs = sum_obs + 1
  end if
  if (mod_x(i_int).eq.0)then
   sum_mod = sum_mod + 1
  end if
  end do
  print("sum_obs = " + sum_obs)
  print("sum_mod = " + sum_mod)
   
  x_new = new(sum_obs,float)
  y_new = new(sum_mod,float)
  i_num = 0
  i_end = 0
  i_beg = 0
  do i_int = 0,dims_obs_x-1
  i_end = i_beg + floattointeger(obs_x(i_int))
  do x_int = i_beg,i_end-1
  x_new(i_num) = i_int
  i_num = i_num + 1
  end do
  if (i_beg.eq.i_end)then
   x_new(i_num) = i_int
   i_num = i_num + 1
  end if
  i_beg = i_end
  end do
  i_num = 0
  i_end = 0
  i_beg = 0
  do i_int = 0,dims_mod_x-1
  i_end = i_beg + floattointeger(mod_x(i_int))
  do x_int = i_beg,i_end-1
  y_new(i_num) = i_int
  i_num = i_num + 1
  end do
  if (i_beg.eq.i_end)then
   y_new(i_num) = i_int
   i_num = i_num + 1
  end if
  i_beg = i_end
  end do
  
  x = x_new
  y = y_new


  wks = gsn_open_wks(wtype,"histo_mod")       ; 9 bars of values 3,2,5,3,4,1,0,1,4

  res = True
; colors = (/(/"maroon4","firebrick","yellow","paleturquoise1","darksalmon","darkkhaki"/),\
;            (/ "orange","forestgreen","tan3","navyblue","hotpink","plum3"/)/)
; colors         = (/ "grey40","grey45","grey50","grey55","grey60","grey65","grey70","grey75","grey80","grey85"/)
  colors         = read_colormap_file("NCV_rainbow2")


  res@gsnHistogramBarColors = colors
  res@tiMainString = "Obs vs Full-SBM"
  res@tiXAxisString         = "dBZ"
  res@tiYAxisString         = "Number of Observations"
  res@trYMinF = 0.
  res@trYMaxF = 1.1*max((/obs_x,mod_x/))
;---Specify the bin intervals to use.
  nx  = dimsizes(x)
  ny  = dimsizes(y)
  lxy = max( (/ nx,ny /) )

  xy = new( (/2,lxy/), typeof(x), getVarFillValue(x))

  xy(0,0:nx-1) = (/ x /)     ;   (/... /) no meta data
  xy(1,0:ny-1) = (/ y /)  
  
; x := xy(0,0:nx-1)

;Let x and y be 1D but of different lengths:
  res@gsnDraw = False
  res@gsnFrame = False
  res@tmXBMode = "Explicit"
; res@tmXBLabels = int_radar(::3)
  res@tmXBLabels = int_radar
  res@tmXBLabelStride                = 1     ; 2 would be every other x-label
  res@gsnHistogramBinIntervals  = ispan(0,dimsizes(obs_x),1)
  res@tmXBLabelFontHeightF     = 0.01   ; make labels smaller
  plot = gsn_histogram(wks,x,res)
; res@tmXBValues = plot@MidBarLocs ; this makes the labels come out in the middle of the interval.
  print("plot@MidBarLocs = " + plot@MidBarLocs)
  print("res@tmXBLabels = " + res@tmXBLabels )
; print("tmXBValues = " + res@tmXBValues)
  res@gsnHistogramCompare        = True
  delete(res@gsnDraw)
  delete(res@gsnFrame)
  plot = gsn_histogram(wks,xy,res)
; plot = gsn_histogram(wks,x,res)
  end