; $Id: dm_stepgrid_bin.pro,v 1.13 2008/03/19 20:46:50 ymqiu Exp $
;#######################################################################
;
; NAME:
;  dm_stepgrid_bin
;
; PURPOSE:
;  bin data according to give xstep and ystep
;
; CATEGORY:
;  general
;
; AUTHOR:
;  Yiming Qiu
;  NIST Center for Neutron Research
;  100 Bureau Drive, Gaithersburg, MD 20899-6102
;  United States
;  yiming.qiu@nist.gov
;  Feb, 2008
;
; LICENSE:
;  The software in this file is written by an employee of
;  National Institute of Standards and Technology
;  as part of the DAVE software project.
;
;  The DAVE software package is not subject to copyright protection
;  and is in the public domain. It should be considered as an
;  experimental neutron scattering data reduction, visualization, and
;  analysis system. As such, the authors assume no responsibility
;  whatsoever for its use, and make no guarantees, expressed or
;  implied, about its quality, reliability, or any other
;  characteristic. The use of certain trade names or commercial
;  products does not imply any endorsement of a particular product,
;  nor does it imply that the named product is necessarily the best
;  product for the stated purpose. We would appreciate acknowledgment
;  if the DAVE software is used or if the code in this file is
;  included in another product.
;
;#######################################################################

; parameters:
;   in_xstep:   xstep
;   in_ystep:   ystep
;   xdat:       input xdat[ne,ndet]    output xdat[xsize]
;   ydat:       input ydat[ne,ndet]    output ydat[ysize]
;   zdat:       input zdat[ne,ndet]    output zdat[xsize,ysize]
; keyword:
;   uniq_xval:  if set, then use the given unique values instead of the steps to determine xbin
;   uniq_yval:  if set, then use the given unique values instead of the steps to determine ybin
;   avgsum:     0: average 1:sum
;   conststep:  0: use the average as the resulted binning values(default) 1: constant binning intervals
;   leave_blank: if set, leave empty column or rows as they are
;   ewid:       energy width, used when e-integrated intensity option is chose
;               input ydat=sum(I*ewid)
;               output ydat=sum(I*ewid)/sum(ewid)
;   zerr:       input error bar of zdat[ne,ndet], ouput zerr[xsize,ysize]
pro dm_stepgrid_bin,in_xstep,in_ystep,xdat,ydat,zdat,uniq_xval=uniq_xval,uniq_yval=uniq_yval,$
    debug=debug,avgsum=avgsum,conststep=conststep,ewid=ewid,leave_blank=leave_blank,zerr=zerr
    if n_elements(avgsum) eq 0 then avgsum=0   ;default using average in binning
    if n_elements(conststep) eq 0 then conststep = 0 ;default is using average as the binned value
    if keyword_set(conststep) then leave_balnk=1
    if keyword_set(debug) then current=systime(/sec)
    nxdat   = n_elements(xdat)
    ok_eint = n_elements(ewid) eq nxdat
    ok_zerr = n_elements(zerr) eq nxdat
    hstep   = abs(in_xstep)
    vstep   = abs(in_ystep)
    if (hstep eq 0) or (n_elements(uniq_xval) ne 0) then hstep=0.1
    if (vstep eq 0) or (n_elements(uniq_yval) ne 0) then vstep=0.1
    xmax    = max(xdat,min=xmin)
    ymax    = max(ydat,min=ymin)
    xmax    = xmax+hstep/2.
    xmin    = xmin-hstep/2.
    ymax    = ymax+vstep/2.
    ymin    = ymin-vstep/2.
    xsize   = (xmax-xmin)/hstep
    ysize   = (ymax-ymin)/vstep
    if xsize gt fix(xsize,type=3) then xsize=(2.>(xsize+0.5))
    if ysize gt fix(ysize,type=3) then ysize=(2.>(ysize+0.5))
    xsize   = fix(xsize,type=3) & ysize=fix(ysize,type=3)
    xlookup = xmin+findgen(xsize+1)*hstep
    ylookup = ymin+findgen(ysize+1)*vstep
    if n_elements(uniq_xval) ne 0  then begin
       xlookup = uniq_xval[uniq(uniq_xval,sort(uniq_xval))]
       xlookup = [xlookup,xmax]
    endif
    if n_elements(uniq_yval) ne 0 then begin
       ylookup = uniq_yval[uniq(uniq_yval,sort(uniq_yval))]
       ylookup = [ylookup,ymax]
    endif
    xsize = n_elements(xlookup)
    ysize = n_elements(ylookup)
    tmp_x = fltarr(xsize)
    tmp_y = fltarr(ysize)
    tmp_z = fltarr(xsize,ysize)
    num_z = lonarr(xsize,ysize)
    num_x = lonarr(xsize) & num_y=lonarr(ysize)
    if ok_eint then tmp_w = fltarr(xsize,ysize)
    if ok_zerr then tmp_e = fltarr(xsize,ysize)
    for i=0L,nxdat-1L do begin
        ix = value_locate(xlookup , xdat[i])
        iy = value_locate(ylookup , ydat[i])
        tmp_x[ix] = tmp_x[ix]+xdat[i] & num_x[ix] = num_x[ix]+1
        tmp_y[iy] = tmp_y[iy]+ydat[i] & num_y[iy] = num_y[iy]+1
        tmp_z[ix,iy] = tmp_z[ix,iy]+zdat[i] & num_z[ix,iy] = num_z[ix,iy]+1
        if ok_eint then tmp_w[ix,iy] = tmp_w[ix,iy]+ewid[i]
        if ok_zerr then tmp_e[ix,iy] = tmp_e[ix,iy]+zerr[i]*zerr[i]
    endfor
    indx = where(num_x ne 0,countx,complement=indx1,ncomplement=countx1)
    indy = where(num_y ne 0,county,complement=indy1,ncomplement=county1)
    if keyword_set(leave_blank) then begin
       if countx1 ne 0 then begin
          num_x[indx1]=1.0 & tmp_x[indx1]=xlookup[indx1]
       endif
       if county1 ne 0 then begin
          num_y[indy1]=1.0 & tmp_y[indy1]=ylookup[indy1]
       endif
       indx=-1 & indx1=-1 & indy=-1 & indy1=-1
    endif else begin
       if countx ne 0 then begin
          tmp_x = tmp_x[indx]   & num_x = num_x[indx]
          tmp_z = tmp_z[indx,*] & num_z = num_z[indx,*]
          if ok_eint then tmp_w = tmp_w[indx,*]
          if ok_zerr then tmp_e = tmp_e[indx,*]
          indx = -1
       endif
       if county ne 0 then begin
          tmp_y = tmp_y[indy]   & num_y=num_y[indy]
          tmp_z = tmp_z[*,indy] & num_z=num_z[*,indy]
          if ok_eint then tmp_w = tmp_w[*,indy]
          if ok_zerr then tmp_e = tmp_e[*,indy]
          indy = -1
       endif
    endelse
    if keyword_set(conststep) and n_elements(uniq_xval) eq 0 then begin
       xdat  = xmin+(findgen(xsize-1)+0.5)*hstep
       tmp_z = tmp_z[0:xsize-2,*]
       num_z = num_z[0:xsize-2,*]
       if ok_zerr then tmp_e = tmp_e[0:xsize-2,*]
       if ok_eint then tmp_w = tmp_w[0:xsize-2,*]
    endif else $
       xdat = temporary(tmp_x)/temporary(num_x)
    if keyword_set(conststep) and n_elements(uniq_yval) eq 0 then begin
       ydat = ymin+(findgen(ysize-1)+0.5)*vstep
       tmp_z = tmp_z[*,0:ysize-2]
       num_z = num_z[*,0:ysize-2]
       if ok_zerr then tmp_e = tmp_e[*,0:ysize-2]
       if ok_eint then tmp_w = tmp_w[*,0:ysize-2]
    endif else $
       ydat = temporary(tmp_y)/temporary(num_y)
    zdat = temporary(tmp_z)
    if ok_zerr then zerr = sqrt(temporary(tmp_e))
    index = where(num_z ne 0,countz)
    if (countz ne 0) then begin
       if ok_eint then begin
          zdat[index] = zdat[index]/tmp_w[index]
          if ok_zerr then zerr[index] = zerr[index]/tmp_w[index]
          tmp_w = 0
       endif else if avgsum eq 0 then begin
          zdat[index] = zdat[index]/num_z[index]    ;average
          if ok_zerr then zerr[index] = zerr[index]/num_z[index]
       endif
    endif
    index = -1
    index = where(num_z eq 0,countz)
    if (countz ne 0) and (not keyword_set(leave_blank)) then begin
       zdat[index] = !values.f_nan
       if ok_zerr then zerr[index] = !values.f_nan
    endif
    index = -1
    if keyword_set(debug) then print,'stepgrid_bin finished in ',systime(/sec)-current,' secs.'
end