; $Id: dm_choose_slicedata.pro,v 1.25 2008/03/19 20:46:50 ymqiu Exp $
;#######################################################################
;
; NAME:
;  dm_choose_slicedata
;
; PURPOSE:
;  choose slice data
;
; CATEGORY:
;  dcs_mslice
;
; AUTHOR:
;  Yiming Qiu
;  NIST Center for Neutron Research
;  100 Bureau Drive, Gaithersburg, MD 20899-6102
;  United States
;  yiming.qiu@nist.gov
;  February, 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:
;   None
; keywords:
;   (xyz)dat: data to be returned
;   zerr:     error bar of zdat
;   (xy)tit:  label of the data
;   zran:     for z data only, contour plot color bar range
;   zmin:     for z data only
pro dcs_mslice::dm_choose_slicedata,xdat=xdat,ydat=ydat,zdat=zdat,zerr=zerr,zran=zran,xtit=xtit,ytit=ytit,zmin=zmin,zmax=zmax,debug=debug
    if keyword_set(debug) then begin
       print,'&&choose slice data...'
       current=systime(/sec)
    endif
    zdat   = (*self.projStrPtr).qty
    uniq_E = 0 ;uniq Energy value
    if arg_present(zerr) then begin
       zerr    = (*self.projStrPtr).err
       ok_zerr = 1b
    endif else ok_zerr = 0b

    ;ewid
    if self.eint_yn and (self.samp_typ ne 2) then $
       self->dm_choose_trajdata,-2,data=ewid
    ok_ewid = (n_elements(ewid) ne 0)
    viewlbs = [self.view_u1lab[(self.samp_typ ge 1)],self.view_u2lab[(self.samp_typ ge 1)],self.view_u3lab]
    case self.samp_typ of
       0:begin      ;powder
         self->dm_choose_trajdata,self.view_u1typ,data=xdat,label=xtit,debug=debug
         self->dm_choose_trajdata,self.view_u2typ,data=ydat,label=ytit,debug=debug
         if xtit eq 'E (meV)' or xtit eq 'Energy' then uniq_E=xdat[*,0]
         if ytit eq 'E (meV)' or ytit eq 'Energy' then uniq_E=ydat[*,0]
         end
       1:begin      ;single crystal
         self->dm_choose_trajdata,3+((self.disp_perp+1) mod 3),data=xdat,label=xtit,debug=debug
         self->dm_choose_trajdata,3+((self.disp_perp+2) mod 3),data=ydat,label=ytit,debug=debug
         if xtit eq 'E (meV)' or xtit eq 'Energy' then uniq_E=xdat[*,0]
         if ytit eq 'E (meV)' or ytit eq 'Energy' then uniq_E=ydat[*,0]
         ;thickness
         if total(finite(self.disp_tran)) ne 0 then begin
            self->dm_choose_trajdata,3+self.disp_perp,index=index1,range=self.disp_tran,empty=empty,debug=debug
            if empty then begin
               emptyitem = viewlbs[self.disp_perp]
               index1 = -1 & tmp = temporary(index1)
            endif
         endif
         end
       2:begin      ;diffuse scattering single E
         self->dm_choose_trajdata,3,data=xdat,label=xtit,debug=debug
         self->dm_choose_trajdata,4,data=ydat,label=ytit,debug=debug
         end
       3:begin      ;diffuse scattering multiple E
         self->dm_choose_trajdata,3+((self.disp_perp+1) mod 3),data=xdat,label=xtit,debug=debug
         self->dm_choose_trajdata,3+((self.disp_perp+2) mod 3),data=ydat,label=ytit,debug=debug
         if xtit eq 'E (meV)' or xtit eq 'Energy' then uniq_E=xdat[*,0,0]
         if ytit eq 'E (meV)' or ytit eq 'Energy' then uniq_E=ydat[*,0,0]
         ;thickness
         if total(finite(self.disp_tran)) ne 0 then begin
            self->dm_choose_trajdata,3+self.disp_perp,index=index1,range=self.disp_tran,empty=empty,debug=debug
            if empty then begin
               emptyitem = viewlbs[self.disp_perp]
               index1 = -1 & tmp = temporary(index1)
            endif
         endif
         end
       else:
    endcase

    ;check extra hkl range
    n_hklran = total(finite(self.disp_hklran))
    if (self.samp_typ gt 0) and (self.hklran_yn eq 1) and (n_hklran ne 0) then begin
       self->dm_choose_trajdata,self.disp_hkl,index=index2,range=self.disp_hklran,empty=empty,debug=debug
       if empty then begin
          if n_elements(emptyitem) eq 0 then emptyitem = (['H','K','L'])[self.disp_hkl] $
          else emptyitem = emptyitem+', '+(['H','K','L'])[self.disp_hkl]
          index2 = -1 & tmp = temporary(index2)
       endif
    endif

    n_index1 = n_elements(index1)
    n_index2 = n_elements(index2)
    if (n_index1 ne 0) or (n_index2 ne 0) then begin
       if n_index1 eq 0 then index = temporary(index2) $
       else if n_index2 eq 0 then index=temporary(index1) $
       else index = dm_common(index1,index2,/no_copy)
       if index[0] eq -1 then begin
          emptyitem = viewlbs[self.disp_perp]+' and '+(['H','K','L'])[self.disp_hkl]
       endif else if n_elements(zdat) ne n_elements(index) then begin
          xdat = xdat[index]
          ydat = ydat[index]
          zdat = zdat[index]
          if ok_zerr then zerr = zerr[index]
          if ok_ewid then ewid = ewid[index]
       endif
       index = -1
    endif

    if n_elements(emptyitem) ne 0 then $
       ok = dialog_message(['The specified '+emptyitem+' range covers no data.', $
               'The entire possible range is used instead.'],dialog_parent=self.tlb)

    ;check folding flag
    if (self.disp_fold eq 1) and ((self.samp_typ eq 1) or (self.samp_typ eq 3)) then begin
       ;determine on which side of the axis the data has wider range
       ;always fold to greater range side unless specified in the hrange
       index = where(xdat gt 0,count)
       if count ne 0 then begin
          x_max = max(xdat[index],min=x_min)
          p_ran = x_max-x_min
          index = -1
       endif else p_ran=0.0
       index = where(xdat lt 0,count)
       if count ne 0 then begin
          x_max = max(xdat[index],min=x_min)
          n_ran = x_max-x_min
          index = -1
       endif else n_ran=0.0
       if p_ran ge n_ran then p_n=1 else p_n=0
       if total(finite(self.disp_hran)) eq 2 then $
          if self.disp_hran[1] gt self.disp_hran[0] then $
             if self.disp_hran[1] gt (-self.disp_hran[0]) then p_n=1 else p_n=0
       if p_n then $
          index = where(xdat lt 0,count) $
       else $
          index = where(xdat gt 0,count)
       if count ne 0 then xdat[index]=(-xdat[index])
       index = -1
    endif

    ;horizontal range and vertical range
    n_xran = total(finite(self.disp_hran))
    n_yran = total(finite(self.disp_vran))
    if n_xran eq 2 then begin
       if self.disp_hran[1] gt self.disp_hran[0] then begin
          index=where((xdat le self.disp_hran[1]) and (xdat ge self.disp_hran[0]),count)
          if count gt 0 then begin
             xdat = xdat[index]
             ydat = ydat[index]
             zdat = zdat[index]
             if ok_ewid then ewid = ewid[index]
             if ok_zerr then zerr = zerr[index]
             index = -1
          endif
       endif
    endif
    if n_yran eq 2 then begin
       if self.disp_vran[1] gt self.disp_vran[0] then begin
          index  = where((ydat le self.disp_vran[1]) and (ydat ge self.disp_vran[0]),count)
          if count gt 0 then begin
             xdat = xdat[index]
             ydat = ydat[index]
             zdat = zdat[index]
             if ok_ewid then ewid = ewid[index]
             if ok_zerr then zerr = zerr[index]
             index = -1
          endif
       endif
    endif
    if n_xran lt 2 then begin
       if finite(self.disp_hran[0]) then begin
          index = where(xdat ge self.disp_hran[0],count)
          if count gt 0 then begin
             xdat = xdat[index]
             ydat = ydat[index]
             zdat = zdat[index]
             if ok_ewid then ewid = ewid[index]
             if ok_zerr then zerr = zerr[index]
             index = -1
          endif
       endif
       if finite(self.disp_hran[1]) then begin
          index = where(xdat le self.disp_hran[1],count)
          if count gt 0 then begin
             xdat = xdat[index]
             ydat = ydat[index]
             zdat = zdat[index]
             if ok_ewid then ewid = ewid[index]
             if ok_zerr then zerr = zerr[index]
             index = -1
          endif
       endif
    endif
    if n_yran lt 2 then begin
       if finite(self.disp_vran[0]) then begin
          index = where(ydat ge self.disp_vran[0],count)
          if count gt 0 then begin
             xdat = xdat[index]
             ydat = ydat[index]
             zdat = zdat[index]
             if ok_ewid then ewid = ewid[index]
             if ok_zerr then zerr = zerr[index]
             index = -1
          endif
       endif
       if finite(self.disp_vran[1]) then begin
          index = where(ydat le self.disp_vran[1],count)
          if count gt 0 then begin
             xdat = xdat[index]
             ydat = ydat[index]
             zdat = zdat[index]
             if ok_ewid then ewid = ewid[index]
             if ok_zerr then zerr = zerr[index]
             index = -1
          endif
       endif
    endif

    ;rebinning according to the steps
    ;rebin even if hstep and vstep are not specified for faster contour plotting
    nogrid = 0
    xmax   = max(xdat, min=xmin)
    ymax   = max(ydat, min=ymin)
    hstep  = (xmax-xmin)/1000.
    vstep  = (ymax-ymin)/1000.
    if (xtit eq 'E (meV)') or (xtit eq 'Energy') then begin
       nogrid = 1
       index  = where((uniq_E ge xmin) and (uniq_E le xmax),count)
       if count ne 0 then uniq_E=uniq_E[index]
       index  = -1
    endif
    if (ytit eq 'E (meV)') or (ytit eq 'Energy') then begin
       nogrid = 2
       index  = where((uniq_E ge ymin) and (uniq_E le ymax),count)
       if count ne 0 then uniq_E=uniq_E[index]
       index  = -1
    endif
    if finite(self.disp_hstep) then begin
       if (self.disp_hstep ne 0) then begin
          hstep  = abs(self.disp_hstep)
          if (xtit eq 'E (meV)') or (xtit eq 'Energy') then nogrid=0
       endif
    endif
    if finite(self.disp_vstep) then begin
       if (self.disp_vstep ne 0) then begin
          vstep  = abs(self.disp_vstep)
          if (ytit eq 'E (meV)') or (ytit eq 'Energy') then nogrid=0
       endif
    endif
    if keyword_set(debug) then help,xdat,nogrid,uniq_E
    if nogrid eq 1 then extra={uniq_xval:uniq_E} $
    else if nogrid eq 2 then extra={uniq_yval:uniq_E}

    if ok_ewid then begin
       dE = 0
       if (xtit eq 'E (meV)') or (xtit eq 'Energy') then dE=hstep
       if (ytit eq 'E (meV)') or (ytit eq 'Energy') then dE=vstep
       if ((self.samp_typ eq 1) and $
           (([self.view_u1[3],self.view_u2[3],self.view_u3[3]])[self.disp_perp] eq 1)) or $
          ((self.samp_typ eq 3) and (self.disp_perp eq 2)) then begin
             tmp_ran = self.disp_tran
             emax = max((*self.projStrPtr).en,min=emin)
             if finite(tmp_ran[0]) eq 0 then tmp_ran[0]=-1e30
             if finite(tmp_ran[1]) eq 0 then tmp_ran[1]=1e30
             if tmp_ran[0] lt emin then tmp_ran[0]=emin
             if tmp_ran[1] gt emax then tmp_ran[1]=emax
             if tmp_ran[1] gt tmp_ran[0] then dE=tmp_ran[1]-tmp_ran[0] $
             else dE=emax-emin
       endif
       if dE le 0 then dE=1.0
    endif

    dm_stepgrid_bin,hstep,vstep,xdat,ydat,zdat,zerr=zerr,avgsum=self.bin_avgsum, conststep=self.bin_conststep,ewid=ewid,debug=debug,_extra=extra

    if ok_ewid then zdat=zdat*dE

    ;smooth regular grid data
    index  = where(1b-finite(zdat),count)   ;bins with no data
    refill = (count ne 0)  ;flag for refilling
    if finite(self.disp_smooth) and (not ok_zerr) then begin    ;smooth data only on plotting
       tmp = size(zdat,/dimension)
       if (n_elements(tmp) gt 1) and (min(tmp) gt 3) and (self.disp_smooth gt 0) then begin
          if float(!version.release) ge 5.6 then begin
             refill = 0
             for i=0, self.disp_smooth-1 do $
                 zdat = convol(zdat,self.disp_sm_kernel,total(self.disp_sm_kernel),/edge_truncate,/nan)
          endif else begin
             ;put 0 in bins with no data
             if refill then zdat[index]=0
             for i=0, self.disp_smooth-1 do $
                 zdat = convol(zdat,self.disp_sm_kernel,total(self.disp_sm_kernel),/edge_truncate)
          endelse
       endif
    endif
    ;intensity range
    oldzmax  = max(zdat,min=oldzmin,/nan)
    if refill and (float(!version.release) lt 5.6) then begin   ;refill the empty bins
       zdat[index] = oldzmin-1e5
       index = -1
    endif
    zmax = oldzmax & zmin = oldzmin
    if finite(self.disp_iran[0]) then $
       if self.disp_iran[0] lt oldzmax then zmin=self.disp_iran[0]
    if finite(self.disp_iran[1]) then $
       if self.disp_iran[1] gt oldzmin then zmax=self.disp_iran[1]
    zran = [zmin,zmax]
    if keyword_set(debug) then print,'&&choose slice data finished in ',systime(/sec)-current,' sec.'
end