; $Id: dm_load_inx.pro,v 1.3 2008/03/19 20:46:50 ymqiu Exp $
;#######################################################################
;
; NAME:
;  dm_load_inx
;
; PURPOSE:
;  load INX files
;
; 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.
;
;#######################################################################

pro dm_load_inx,unit,zdat=zdat,zerr=zerr,xdat=xdat,ydat=ydat,eief=eief,comment=comment,$
    tchannel=tchannel,temperature=temperature,headeronly=headeronly

    tmp0   = ''
    tmp1   = ''
    tmp1_1 = ''
    tmp1_2 = ''
    nth    = 0L
    while not eof(unit) do begin
        readf,unit,tmp0
        if nth eq 0 then begin
           temp = dm_to_number(strsplit(tmp0,' ',/extract),/long)
           tchannel = temp[[0,n_elements(temp)-1]]
           readf,unit,tmp1
           ncomment = 1
           readf,unit,tmp1_1
           temp = dm_to_number(strsplit(tmp1_1,' ',/extract))
           ;figure out home many comment lines
           while total(finite(temp)) ne n_elements(temp) do begin
              ncomment = ncomment+1
              tmp1 = strtrim(tmp1,2)+' '+strtrim(tmp1_1,2)
              readf,unit,tmp1_1
              temp = dm_to_number(strsplit(tmp1_1,' ',/extract))
           endwhile
           if n_elements(temp) ne 9 then begin
              readf,unit,tmp1_2
              tmp1_1 = tmp1_1+' '+tmp1_2
           endif
           tmp2 = dm_to_number(strsplit(tmp1_1,' ',/extract))
           comment = strtrim(tmp1,2)
           eief = tmp2[1]
           temperature = tmp2[3]

           if keyword_set(headeronly) then begin
              free_lun,unit
              return
           endif
           tmp3 = dblarr(3,tchannel[1])
           readf,unit,tmp3
           xdat = tmp2[0]
           ydat = reform(tmp3[0,*])
           zdat = reform(tmp3[1,*])
           zerr = reform(tmp3[2,*])
        endif else begin
           for i=0,ncomment-1 do $
               readf,unit,tmp1
           readf,unit,tmp2
           readf,unit,tmp3
           xdat = [xdat,tmp2[0]]
           zdat = [[zdat],[reform(tmp3[1,*])]]
           zerr = [[zerr],[reform(tmp3[2,*])]]
        endelse
        nth = nth+1L
    endwhile
    zdat = transpose(zdat)
    zerr = transpose(zerr)
    ydat = -ydat            ;change the sign of energy transfer
    free_lun,unit
end

;this procedure reads an INX file for dcs_mslice program
;Parameters:
;   open_path:  file path, string
;   open_file:  file name(s), string or string array
;Keywords:
; input:
;   avgsum:     0:average 1:sum, effective when diffuse keyword is sest, default is self.bin_avgsum
;   background: if set, the loaded files are regarded as background files and will be summed
;   diffetyp:   only needed when diffuse keyword is set 0:single energy  1:multiple energy
;   diffuse:    if set, diffuse scattering mode
;   parent:     message widget parent,  default is self.tlb
;   title:      message title
;   eint_yn:    only applicable for diffuse mode, if set, I=sum(I[i,j]*dE[i,j]) over the specified time channels
;   e_range:    energy range - diffuse keyword set, default is self.e_range
;   t_chan:     0:all 1:elastic 2:specify erange  - efffective when diffuse keyword set, default is self.e_bin[self.samp_typ]
; output:
;   error:      error flag 1 or 0
;   qty:        [ne,ndet] or [ne,ndet,nfiles]-diffuse keywords set [ne,ndet] if nfiles=1
;   dqty:       [ne,ndet] or [ne,ndet,nfiles]-diffuse keywords set [ne,ndet] if nfiles=1
;   ang:        [nfiles] - diffuse keyword set, not in inx file, need to be typed in by the user
;   eief:       fixed energy, Ei for direct geometry, Ef for inverse geometry
;   emean:      [ne]
;   ewid:       diffuse:ewid[tchans] ; powder,single:[ntchan], the corresponding energy width for the time channels
;   info:       {ne:ne, ndet:ndet, ei:ei, temperature:temperature}  used for checking file compatibility
;   weight:     weight or duration info , not in the file, need to be typed in by the user
;   comment:     header info
;   temperature:temerature
;   kfactor:    [ntchan]-(ki/kf)  or a scalor
pro dcs_mslice::dm_load_inx,open_path,open_file,ang=ang,avgsum=avgsum,background=background,$
    comment=comment,diffetyp=diffetyp,diffuse=diffuse,eief=eief,e_range=e_range,eint_yn=eint_yn,$
    emean=emean,error=error,ewid=ewid,info=info,kfactor=kfactor,parent=parent,qty=qty,dqty=dqty,$
    t_chan=t_chan,title=title

    if keyword_set(diffuse) then begin
       if n_elements(avgsum)  eq 0  then avgsum   = self.bin_avgsum
       if n_elements(t_chan)  eq 0  then t_chan   = self.e_bin[self.samp_typ]
       if n_elements(e_range) ne 2  then e_range  = self.e_range
       if n_elements(diffetyp) eq 0 then diffetyp = (self.samp_typ eq 3)
    endif
    if n_elements(parent) eq 0 then parent = self.tlb
    if n_elements(title)  eq 0 then title  = 'Loading data ...'
    error   = 0b  ;initialize no error
    n_files = n_elements(open_file)

    mesg = widget_base(title=title,/col,/base_align_center,group_leader=parent)

    ;catch and clear possible io errors
    catch, myerror
    if myerror ne 0 then begin
       error = 1b
       catch,/cancel
       if widget_info(mesg,/valid_id) then widget_control,mesg,/destroy
       ok = dialog_message(!error_state.msg,/error)
       if n_elements(unit) ne 0 then free_lun,unit
       return
    end

    ;ask for weight info
    myweight = 1.0
    if arg_present(weight) or (n_files ne 1) then begin
       myweight = !values.f_nan
       while(total(finite(myweight)) ne n_files) do $
            myweight = dm_dialog_listinput(open_file,title='Please input weight info:',default=1,$
                 label='weight',dialog_parent=parent,reset=1,info=['The weight is for summing files.',$
                 "Press OK if you don't know the information."])
    endif

    ;ask for rotation angles
    if keyword_set(diffuse) and (not keyword_set(background)) then begin
       ang = !values.f_nan
       while(total(finite(ang)) ne n_files) do $
          ang = dm_dialog_listinput(open_file,title='Please input rotation angles:',default=0,$
                label='angle',dialog_parent=parent,reset=0)
    endif

    for i=0L,n_files-1L do begin
        file = open_path+self.pathsep+open_file[i]
        if i eq 0 then begin
           label = widget_label(mesg,value='Loading '+file,/dynamic_resize,/align_center)
           dm_center_kid,mesg,parent
           widget_control,mesg,/realize,map=1
        endif else $
           widget_control,label,set_value='Loading '+file,tlb_set_title=strtrim(title,2)+'  '+$
                strtrim(string(i+1),2)+'/'+strtrim(string(n_files),2)

        openr,unit,file,/get_lun,error=fileerror
        if fileerror ne 0 then begin
           error = 1b
           widget_control,mesg,/destroy
           ok = dialog_message("Can't open "+file,/error,dialog_parent=parent)
           return
        endif

        ;read the inx file
        dm_load_inx,unit,zdat=zdat,zerr=dzdat,xdat=xdat,ydat=ydat,comment=comment,eief=eief,temperature=temperature    ;xdat:tth, ydat:en
        ndet = n_elements(xdat)
        nen  = n_elements(ydat)
        zdim = size(zdat,/dimension)

        if n_elements(info) eq 0 then begin
           info={ndet:ndet,nen:nen,eief:eief,temperature:temperature}
        endif else begin
           if (info.ndet ne ndet) or (info.nen ne nen) or (abs(info.eief-eief) gt 1e-4) then begin
              ok = dialog_message('Selected data files are of different types.',/error,dialog_parent=parent)
              error = 1b                       ;error occurs
              widget_control,mesg,/destroy     ;clean up
              return
           endif
        endelse

        if i eq 0 then begin
           emean = temporary(ydat)
           tth   = temporary(xdat)
           dtth  = fltarr(zdim[0])
           psi   = fltarr(zdim[0])
           dpsi  = fltarr(zdim[0])
           ewid  = emean[1:nen-1]-emean[0:nen-2]
           ewid  = [ewid[0], (ewid[0:nen-3]+ewid[1:nen-2])/2,ewid[nen-2]]   ;reconstruct ewid

           if keyword_set(diffuse) then begin
              weight  = myweight
              all_tch = 1
              kfactor = 1.0
              ;search the proper time channels to be summed or averaged
              case t_chan of
                    0:          ;all time channels
                    1:  begin   ;elastic peak
                        tmp_y    = total(zdat,1)
                        tmp_ny   = n_elements(tmp_y)
                        tmp_x    = findgen(tmp_ny)
                        tmp_ymax = max(tmp_y,imax)
                        ;shift the peak approximately to the center channel
                        tmp_y    = shift(tmp_y,tmp_ny/2.-imax)
                        estimates= [tmp_ymax,tmp_ny/2.,20.0,0.0]
                        result   = gaussfit(tmp_x,tmp_y,params,estimates=estimates,nterms=n_elements(estimates),chisq=chisq)
                        if keyword_set(debug) then begin
                           print,params
                           window,/free
                           plot,tmp_x,tmp_y,psym=6,xr=tmp_ny/2.+[-params[2],params[2]]*3.,title='chisq='+strtrim(string(chisq),2)
                           oplot,tmp_x,result
                           wait,2 & wdelete,!D.WINDOW
                        endif
                        tmp_x = -1 & tmp_y = -1
                        params[2] = abs(params[2])
                        if (abs(params[1]-tmp_ny/2.) gt tmp_ny/5.) or (params[2] le 1) then begin
                           ok=dialog_message(["Can't locate the elastic peak from the first file.","All time channels are used instead."],$
                                dialog_parent=parent)
                           widget_control,label,/show,/input_focus
                        endif else begin
                           all_tch = 0     ;fitting is good
                           tchans  = params[1]-(tmp_ny/2.-imax)-fix(2.5*params[2])+indgen(fix(5.0*params[2])+1)
                           index   = where(tchans lt 0, count)
                           if count ne 0 then tchans[index] = tchans[index]+tmp_ny
                           index   = where(tchans ge tmp_ny,count)
                           if count ne 0 then tchans[index] = tchans[index]-tmp_ny
                           emin    = min(emean[tchans],max=emax)
                           e_range = [emin,emax]
                        endelse
                        end
                    2:  begin   ;specified energy range
                        if n_elements(e_range) eq 2 then $
                           if total(finite(e_range)) eq 2 then begin
                              tchans = where((emean ge e_range[0]) and (emean le e_range[1]),count)
                              if count ne 0 then begin
                                 all_tch = 0
                                 if arg_present(kfactor) then begin
                                    if self.instrgeom eq 0 then begin   ;direct geometry
                                       kfactor = sqrt(eief/(eief-mean(e_range))) ;ki/kf_avg
                                    endif else begin   ;inverse geometry
                                       kfactor = sqrt((eief+mean(e_range))/eief) ;ki_avg/kf
                                    endelse
                                 endif
                                 break
                              endif
                           endif
                        ok = dialog_message(['The energy range specified is invalid.','All time channels are used instead.'],$
                             dialog_parent=parent)
                        widget_control,label,/show,/input_focus
                        end
                    else:
              endcase
              if all_tch eq 1 then begin
                 tchans = indgen(zdim[1])
                 emin   = min(emean,max=emax)
                 e_range = [emin,emax]
              endif
              ewid = ewid[tchans]
              nitem = (avgsum eq 0)?(1.0>(n_elements(ewid))):1.0
              if keyword_set(eint_yn) then factor=1.0 else factor = 1.0/nitem
              if (diffetyp eq 1) and arg_present(kfactor) then begin
                 if self.instrgeom eq 0 then begin   ;direct geometry
                    kfactor = (sqrt(eief/(eief-emean)))[tchans]
                 endif else begin   ;inverse geometry
                    kfactor = (sqrt((eief+emean)/eief))[tchans]
                 endelse
              endif
           endif else begin
              qty    = 0.0
              dqty   = 0.0
              weight = 0.0
           endelse
        endif

        ;add or append data
        if keyword_set(diffuse) then begin
           zdat  = zdat[*,tchans]
           dzdat = dzdat[*,tchans]
           if diffetyp eq 0 then begin
              if keyword_set(eint_yn) then begin
                 for j=0L,n_elements(zdat[*,0])-1 do begin
                     zdat[j,*] = zdat[j,*]*ewid
                     zdat[j,*] = zdat[j,*]*ewid
                 endfor
              endif
              if i eq 0 then begin
                 qty  = total(zdat,2)*factor
                 dqty = sqrt(total(dzdat^2,2))*factor
              endif else begin
                 qty  = [[qty],[total(zdat,2)*factor]]
                 dqty = [[dqty],[sqrt(total(dzdat^2,2))*factor]]
              endelse
           endif else begin
              if i eq 0 then begin
                 qty  = transpose(zdat)
                 dqty = transpose(dzdat)
              endif else begin
                 qty  = [[[qty]],[[transpose(zdat)]]]
                 dqty = [[[dqty]],[[transpose(dzdat)]]]
              endelse
           endelse
        endif else begin
           qty    = temporary(qty)+temporary(zdat*myweight[i])
           dqty   = temporary(dqty)+temporary(dzdat*myweight[i])
           weight = temporary(weight)+myweight[i]
        endelse
    endfor

    if not keyword_set(diffuse) then begin
       if arg_present(kfactor) then begin
          if self.instrgeom eq 0 then begin   ;direct geometry
             kfactor = sqrt(eief/(eief-emean))
          endif else begin   ;inverse geometry
             kfactor = sqrt((eief+emean)/eief)
          endelse
       endif
       if weight eq 0 then weight = 1.0
       qty  = qty/weight
       dqty = dqty/weight
       qty  = transpose(qty)   ;->[nen,ndet]
       dqty = transpose(dqty)  ;->[nen,ndet]
    endif else begin
       if diffetyp eq 1 then emean=emean[tchans] $
       else begin
          if t_chan eq 2 then emean=mean(emean[tchans]) $
          else emean = 0.0
       endelse
       if keyword_set(background) then begin    ;sum the files
          ind    = (size(qty))[0]
          totw   = total(weight)
          if totw le 0 then totw = 1.0
          weight = weight/totw
          for i=0,n_elements(weight)-1 do begin
             qty[*,i]  = qty[*,i]*weight[i]
             dqty[*,i] = (dqty[*,i]*weight[i])^2
          endfor
          qty    = total(qty,ind)
          dqty   = sqrt(total(dqty,ind))
          weight = totw
       endif
    endelse

    ptr_free,self.detPosPtr
    self.detPosPtr = ptr_new({two_theta:temporary(tth),dtwo_theta:temporary(dtth),$
                     psi:temporary(psi),dpsi:temporary(dpsi)})

    widget_control,mesg,/destroy
end

;the following procedure is for loading detector dependent background file for dcs_mslice
pro dcs_mslice::dm_load_inx_detbackground,open_path,open_file,parent=parent,error=error,qty=qty,$
    dqty=dqty,comment=comment
    self->dm_load_inx,open_path,open_file,qty=qty,dqty=dqty,t_chan=2,e_range=self.bgeran,avgsum=0,$
        title='loading background file...',parent=parent,error=error,/diffuse,/background,diffetyp=0,$
        comment=comment
end

;the following procedure is for loading time channel background file for dcs_mslice
pro dcs_mslice::dm_load_inx_tchanbackground,open_path,open_file,parent=parent,error=error,qty=qty,$
    dqty=dqty,comment=comment,emean=emean
    self->dm_load_inx,open_path,open_file,qty=qty,dqty=dqty,t_chan=2,e_range=self.dcor_eqran[0:1],$
        avgsum=0,title='loading background file...',parent=parent,error=error,/diffuse,/background,$
        diffetyp=1,comment=comment
end
