; $Id: dm_step_bin.pro,v 1.10 2008/03/19 20:46:50 ymqiu Exp $
;#######################################################################
;
; NAME:
;  dm_step_bin
;
; PURPOSE:
;  bin data according to given step
;
; 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.
;
;#######################################################################

; parameter:
;   in_xstep: step of binning
;   xdat:     array of xdat
; keyword:
;   ydat:     array of ydat, using the same binning as ydat
;   yerr:     array of yerr, using the same binning as xdat
;   uniq_val: if set, then use the given uniq values instead of the xstep to determine the bins
;   avgsum:   0:average 1:sum
;   conststep: 0: use the average as the resulted binning values(default) 1: constant binning intervals
;   ewid:     energy width, used when e-integrated intensity option is choosed
;             input ydat=sum(I*ewid)
;             output ydat=sum(I*ewid)/sum(ewid)
pro dm_step_bin,in_xstep,xdat,ydat=ydat,yerr=yerr,uniq_val=uniq_val,debug=debug,avgsum=avgsum,conststep=conststep,ewid=ewid
    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(debug) then current=systime(/sec)
    xstep = abs(in_xstep)
    if xstep eq 0 then xstep=1.
    nxdat    = n_elements(xdat)
    ok_ydat  = (n_elements(ydat) eq nxdat)
    ok_yerr  = (n_elements(yerr) eq nxdat)
    ok_eint  = (n_elements(ewid) eq nxdat)
    tmp_xmax = max(xdat,min=tmp_xmin)
    xmax  = temporary(tmp_xmax)+0.5*xstep
    xmin  = temporary(tmp_xmin)-0.5*xstep
    xsize = double(xmax-xmin)/xstep
    if xsize gt fix(xsize,type=3) then xsize=(2.>(xsize+0.5))
    if xsize eq 2 then conststep=0
    xsize = fix(xsize,type=3)
    if n_elements(uniq_val) ne 0 then begin
       index = where((uniq_val ge xmin) and (uniq_val lt xmax),count)
       if count ne 0 then uniq_val=uniq_val[index]
       index = -1
       xlookup = uniq_val[uniq(uniq_val,sort(uniq_val))]
       xlookup = [xlookup,xmax]
    endif else $
       xlookup = xmin+findgen(xsize+1)*xstep
    xsize = n_elements(xlookup)
    tmp_x = fltarr(xsize)
    if ok_ydat then tmp_y = fltarr(xsize)
    if ok_yerr then tmp_e = fltarr(xsize)
    if ok_eint then tmp_w = fltarr(xsize)
    num_x = lonarr(xsize)

    for i=0L,nxdat-1L do begin
        ix = value_locate(xlookup , xdat[i])
        tmp_x[ix] = tmp_x[ix]+xdat[i]
        if ok_ydat then tmp_y[ix] = tmp_y[ix]+ydat[i]
        if ok_yerr then tmp_e[ix] = tmp_e[ix]+yerr[i]*yerr[i]
        if ok_eint then tmp_w[ix] = tmp_w[ix]+ewid[i]
        num_x[ix] = num_x[ix]+1
    endfor
    index = where(num_x gt 0,count,complement=zindex,ncomplement=zcount)
    if count ne 0 then begin
       if keyword_set(conststep) and n_elements(uniq_val) eq 0 then begin
          if zcount ne 0 then begin
             num_x[zindex] = 1.0
             if ok_eint then tmp_w[zindex] = 1.0
          endif
          xdat = xmin+(findgen(xsize-1)+0.5)*xstep
          num_x = num_x[0:xsize-2]
          if ok_ydat then tmp_y = tmp_y[0:xsize-2]
          if ok_yerr then tmp_e = tmp_e[0:xsize-2]
          if ok_eint then tmp_w = tmp_w[0:xsize-2]
       endif else begin
          tmp_x = tmp_x[index]
          num_x = num_x[index]
          if ok_ydat then tmp_y = tmp_y[index]
          if ok_yerr then tmp_e = tmp_e[index]
          if ok_eint then tmp_w = tmp_w[index]
          xdat  = tmp_x/num_x
       endelse
       index = -1 & zindex=-1
       if avgsum eq 1 then num_x = 1
       if ok_eint then num_x = tmp_w
       if ok_ydat then ydat  = tmp_y/num_x
       if ok_yerr then yerr  = sqrt(tmp_e)/num_x
    endif
    if keyword_set(debug) then print,'step_bin finished in ',systime(/sec)-current,' secs.'
end