; $Id: dm_det_spurion.pro,v 1.6 2008/03/19 20:46:50 ymqiu Exp $
;#######################################################################
;
; NAME:
;  dm_det_spurion
;
; PURPOSE:
;  plot spurion caused by the scattering from detectors
;
; CATEGORY:
;  dcs_mslice
;
; AUTHOR:
;  Yiming Qiu
;  NIST Center for Neutron Research
;  100 Bureau Drive, Gaithersburg, MD 20899-6102
;  United States
;  yiming.qiu@nist.gov
;  March, 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 dcs_mslice::dm_det_spurion,debug=debug
    case self.ftype of
        0:  begin
            psi_0_only = 1b                           ;calculate only psi=0 detectors
            d_sd  = 4010.                             ;sample to detector distance, in cm
            ;detector is made of stainless steel
            a_det = 3.594                             ;lattice parameter for fcc stainless steel
            q_det = sqrt([3,4,8,11,12])*2.*!pi/a_det  ;[1,1,1],[2,0,0],[2,2,0],[3,1,1],[2,2,2]
            i_det = [0.005,0.003,0.003,0.003,0.003]   ;relative intensity
            end
        1:  begin
            psi_0_only = 0b                           ;calculate all detectors
            d_sd  = 4010.                             ;sample to detector distance
            ;assume stainless steel detector
            a_det = 3.594                             ;lattice parameter for fcc stainless steel
            q_det = sqrt([3,4,8,11,12])*2.*!pi/a_det  ;[1,1,1],[2,0,0],[2,2,0],[3,1,1],[2,2,2]
            i_det = [0.005,0.003,0.003,0.003,0.003]   ;relative intensity
            end
        else:
    endcase

    ;in case the necessary info is absent
    if n_elements(psi_0_only) eq 0 then psi_0_only = 0b
    if n_elements(d_sd)  eq 0 then d_sd  = 4010.    ;sample to detector distance in cm
    if n_elements(a_det) eq 0 then a_det = 3.594    ;lattice parameter in A
    if n_elements(q_det) eq 0 then q_det = sqrt([3,4,8,11,12])*2.*!pi/a_det
    if n_elements(i_det) eq 0 then i_det = fltarr(n_elements(q_det))+0.003

    dPtr    = self.dataStrPtr
    pPtr    = self.projStrPtr
    tPtr    = self.detPosPtr
    kPtr    = self.kfactor
    tth_all = (*tPtr).two_theta
    case self.ftype of
        0:  fwhm = dcs_jcc_resn((*dPtr).info.ch_wl,(*dPtr).info.ch_wl,(*dPtr).info.ch_res,(*dPtr).info.ch_ms)/1000.    ;in meV
        else:begin
            ;fit the scattering data
            tmp_y    = total((*dPtr).qty,2)
            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)
            tmp_x = -1 & tmp_y = -1
            params[2] = abs(params[2])
            if (abs(params[1]-tmp_ny/2) lt tmp_ny/5.) and (params[2] gt 1) then begin
                tmp = long(params[2]*1.17741)
                tchans = [imax-tmp,imax+tmp]
                fwhm = (*dPtr).e[tchans[1]]-(*dPtr).e[tchans[0]]
            endif else $
                fwhm = 0.5
            end
    endcase
    index   = where(((*dPtr).e gt -1.2*fwhm) and ((*dPtr).e lt 1.2*fwhm), count1)  ;elastic channels
    if count1 eq 0 then return
    qty     = (*dPtr).qty[index,*]
    en      = (*dPtr).e[index]
    tth     = ((*tPtr).two_theta)
    if psi_0_only then begin        ;use only psi = 0 data, central bacnk for DCS
        index = where((*tPtr).psi eq 0,count2)
        if count2 ne 0 then begin
            qty = qty[*,index]
            tth = ((*tPtr).two_theta)[index]
        endif else return
    endif else $
        count2 = n_elements(tth)


    tth_all = tth[sort(tth)]
    max_en  = max((*dPtr).e)
    max_tt1 = 0.0 & min_tt1 = 0.0 & max_tt2 = 0.0 & min_tt2 = 0.0
    ipos    = where(tth_all gt 0, cpos)
    if cpos ne 0 then max_tt1 = max(tth_all[ipos],min=min_tt1)
    ineg    = where(tth_all lt 0, cneg)
    if cneg ne 0 then max_tt2 = max(tth_all[ineg],min=min_tt2)


    ext_qty = 0.0
    ext_en  = 0.0
    ext_tth = 0.0

    current = systime(/sec)

    ei      = (*dPtr).eief
    conv    = 81.8042/4./!pi^2    ;2.072
    hpi     = !pi/2.
    for i=0L,count1-1 do begin
        ki  = sqrt((ei-en[i])/conv)
        for j=0L,n_elements(q_det)-1 do begin
            if q_det[j]/2. lt ki then begin
               dtth = 2.*asin(q_det[j]/2./ki)
               if dtth gt hpi then begin
                  new_en = ei-(ki*d_sd/(d_sd-2*d_sd*cos(dtth)))^2*conv
                  if new_en lt max_en then begin
                     d_tth = 2*dtth/!dtor-180.
                     for k=0L,count2-1 do begin
                         new_tth = tth[k]+d_tth
                         if (new_tth lt max_tt1 and new_tth gt min_tt1) or (new_tth lt max_tt2 and new_tth gt min_tt2) then begin
                            ext_tth = [ext_tth,new_tth]
                            ext_en  = [ext_en,new_en]
                            ext_qty = [ext_qty,qty[i,k]*i_det[j]]
                         endif
                         new_tth = tth[k]-d_tth
                         if (new_tth lt max_tt1 and new_tth gt min_tt1) or (new_tth lt max_tt2 and new_tth gt min_tt2) then begin
                            ext_tth = [ext_tth,new_tth]
                            ext_en  = [ext_en,new_en]
                            ext_qty = [ext_qty,qty[i,k]*i_det[j]]
                         endif
                     endfor
                  endif
               endif
            endif
        endfor
    endfor

    if keyword_set(debug) then print,'Spurion calculation is done is ',systime(/sec)-current,' secs.'

    n_ext   = n_elements(ext_qty)
    if n_ext eq 1 then begin
        ok  = dialog_message("There's no spurion caused by the detectors.",dialog_parent=self.tlb,/info)
        return
    endif

    en      = en#(fltarr(count2)+1)
    tth     = (fltarr(count1)+1)#tth
    qty     = [qty[*],ext_qty[1:n_ext-1]]
    en      = [en[*], ext_en[1:n_ext-1]]
    tth     = [tth[*],ext_tth[1:n_ext-1]]
    max_en  = max(en,min=min_en)
    itmp    = where((*dPtr).e le min_en,count)
    min_en  = (*dPtr).e[itmp[count-1]]
    itmp    = where((*dPtr).e ge max_en)
    max_en  = (*dPtr).e[itmp[0]]
    uniq_en = (*dPtr).e[where((*dPtr).e ge min_en and (*dPtr).e le max_en)]
    ewid    = (*dPtr).ewid
    time    = (*dPtr).time
    dm_stepgrid_bin,0.0,0.0,tth,en,qty,uniq_xval=tth_all,uniq_yval=uniq_en,/leave_blank
    qty     = transpose(qty)

    ;replace dataStrPtr, projStrPtr, and detPosPtr by the new spurion data
    titl_lab        = self.titl_lab
    self.titl_lab   = 'Spurion'
    self.dataStrPtr = ptr_new({qty:qty,err:0*qty,e:en,weight:0,eief:ei,info:(*dPtr).info,ewid:ewid,time:time})
    self.detPosPtr  = ptr_new({two_theta:tth,dtwo_theta:0*tth,psi:0*tth,dpsi:0*tth})
    self.projStrPtr = ptr_new()
    self.kfactor    = ptr_new(sqrt(ei/(ei-en)))
    self->dm_calc_projection    ;recalculate projection
    self->dm_plot_slice         ;plot the spurion

    ;return to initial state
    self.titl_lab   = titl_lab
    ptr_free,self.dataStrPtr,self.projStrPtr,self.detPosPtr,self.kfactor
    self.dataStrPtr = dPtr
    self.projStrPtr = pPtr
    self.detPosPtr  = tPtr
    self.kfactor    = kPtr
end