; $Id: dm_proj_spec_samp.pro,v 1.10 2008/03/19 20:46:50 ymqiu Exp $
;#######################################################################
;
; NAME:
;  dm_proj_spec_samp
;
; PURPOSE:
;  transform data points in the form (detector_angles, energy)->(Qx,Qy,Qz)
;  in the spectrometer frame where [1,0,0] is along ki, [0,1,0] in the horizontal scattering plane
;  and [0,0,1] pointing upward
;  then into the sample reference frame by a rotation of the coordinate system
;
; 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:
;   eief:     Ei or Ef according to instrument geometry
;   en:       energy transfer: [ne] for powder and single crystal type, scalor or [ne] for diffuse scattering type
;   det_tt:   two_theata:[ndet]
;   det_dtt:  dtwo_theta:[ndet]
;   det_psi:  psi:[ndet]
;   det_dpsi: dpsi:[ndet]
; keywords:
;   psi:      scalor for single crystal type, [nfile] array for diffuse scattering type; positive if ccw rotation of u from ki
;   type:     0:powder 1:single crystal 2:diffuse scattering
;   diffetyp: combined with type=2   0:single energy 1:multiple energy
;   is_spe:   for matlab style spe file
;   instrgeom:0-direct geometry, 1-inverse geometry
function dm_proj_spec_samp,eief,en,det_tt,det_dtt,det_psi,det_dpsi,psi=psi,type=type,debug=debug,$
  is_spe=is_spe,diffetyp=diffetyp,instrgeom=instrgeom
    if keyword_set(debug) then begin
       print,'##proj spec samp...'
       current=systime(/sec)
    endif
    E_ksq = 2.0721246  ;E_ksq=!dcs_hsq2mn/4./!dpi/!dpi;=2.0721246
    nen = n_elements(en)
    if keyword_set(instrgeom) then begin
       ki = sqrt((eief+en)/E_ksq)
       kf = sqrt(eief/E_ksq)*(fltarr(nen)+1.0)
    endif else begin
       ki = sqrt(eief/E_ksq)*(fltarr(nen)+1.0)
       kf = sqrt((eief-en)/E_ksq)
    endelse
    if n_elements(ki) eq 1 then ki=ki[0]
    if n_elements(kf) eq 1 then kf=kf[0]
    ndet= n_elements(det_tt) & nfile=n_elements(psi)
    if n_elements(type) eq 0 then type=0      ;default-powder type
    case type of
       0:begin  ;powder type -> Q[nen,ndet]
         Q = sqrt((ki*ki+kf*kf)#(fltarr(ndet)+1.0)-2.*(ki*kf)#cos(det_tt*!dtor))
         end
       1:begin  ;single crystal type -> Q[3,nen,ndet]
         Q        = fltarr(3,nen,ndet,/nozero)
         dettt    = det_tt*!dtor
         detpsi   = det_psi*!dtor
         cos2th   = cos(dettt)
         sin2th   = sin(temporary(dettt))
         sinpsi   = sin(detpsi)
         cospsi   = cos(temporary(detpsi))
         Q[0,*,*] = temporary(ki)#(fltarr(ndet)+1.0)-kf#temporary(cos2th)
         if keyword_set(is_spe) then begin  ;spe style definition of psi
            Q[1,*,*] = -kf#(sin2th*temporary(cospsi))
            Q[2,*,*] = -temporary(kf)#(temporary(sin2th)*temporary(sinpsi))
         endif else begin    ;dcs style defininiton of psi (7.6285 for upper bank, -7.6285 for lower bank) and 2th(+/-)
            Q[1,*,*] = kf#(sin2th*sqrt(1-(sinpsi/sin2th)^2))
            Q[2,*,*] = -temporary(kf)#temporary(sinpsi)
         endelse
         sinpsi = sin(psi[0]*!dtor)
         cospsi = cos(psi[0]*!dtor)
         R      = [[cospsi,-sinpsi,0.0],[sinpsi,cospsi,0.0],[0.0,0.0,1.0]]
         ;old and slower code
         ;tmp     = Q
         ;Q[0,*,*]= tmp[0,*,*]*R[0,0]+tmp[1,*,*]*R[0,1]+tmp[2,*,*]*R[0,2]
         ;Q[1,*,*]= tmp[0,*,*]*R[1,0]+tmp[1,*,*]*R[1,1]+tmp[2,*,*]*R[1,2]
         ;Q[2,*,*]= tmp[0,*,*]*R[2,0]+tmp[1,*,*]*R[2,1]+tmp[2,*,*]*R[2,2]
         Q  = reform(Q,3,nen*ndet)
         Q  = temporary(Q)##temporary(R)
         Q  = reform(Q,3,nen,ndet)
         end
       2:begin  ;diffuse scattering type -> Q[3,ndet, nfile] or Q[3,nen,ndet,nfile]
         dettt  = det_tt*!dtor
         detpsi = det_psi*!dtor
         cos2th = cos(dettt)
         sin2th = sin(temporary(dettt))
         sinpsi = sin(detpsi)
         cospsi = cos(temporary(detpsi))
         if diffetyp eq 0 then begin    ;single energy
            Q0      = fltarr(3,ndet,/nozero)
            Q0[0,*] = temporary(ki)-kf*temporary(cos2th)
            if keyword_set(is_spe) then begin  ;spe style definition of psi
               Q0[1,*] = -kf*(sin2th*temporary(cospsi))
               Q0[2,*] = -temporary(kf)*(temporary(sin2th)*temporary(sinpsi))
            endif else begin                   ;dcs style definition of psi
               Q0[1,*] = kf*(sin2th*sqrt(1-(sinpsi/sin2th)^2))
               Q0[2,*] = -temporary(kf)*temporary(sinpsi)
            endelse
         endif else begin               ;multiple energy
            Q0        = fltarr(3,nen,ndet,/nozero)
            Q0[0,*,*] = temporary(ki)#(fltarr(ndet)+1.0)-kf#temporary(cos2th)
            if keyword_set(is_spe) then begin  ;spe style definition of psi
               Q0[1,*,*] = -kf#(sin2th*temporary(cospsi))
               Q0[2,*,*] = -temporary(kf)#(temporary(sin2th)*temporary(sinpsi))
            endif else begin                   ;dcs style definition of psi
               Q0[1,*,*] = kf#(sin2th*sqrt(1-(sinpsi/sin2th)^2))
               Q0[2,*,*] = -temporary(kf)#temporary(sinpsi)
            endelse
         endelse
         R        = fltarr(3,3,nfile)
         sinpsi   = sin(psi*!dtor)
         cospsi   = cos(psi*!dtor)
         R[0,0,*] = cospsi
         R[1,0,*] = -sinpsi
         R[0,1,*] = temporary(sinpsi)
         R[1,1,*] = temporary(cospsi)
         R[2,2,*] = 1.0
         if diffetyp eq 0 then begin    ;single energy
            Q = fltarr(3,ndet,nfile,/nozero)
            for i=0,nfile-1 do $
                Q[*,*,i] = Q0##R[*,*,i]
         endif else begin               ;multiple energy
            Q  = fltarr(3,nen,ndet,nfile,/nozero)
            Q0 = reform(Q0,3,nen*ndet)
            for i=0,nfile-1 do begin
                tmp = Q0##R[*,*,i]
                Q[*,*,*,i] = reform(tmp,3,nen,ndet)
            endfor
         endelse
         end
       else:
    endcase
    if keyword_set(debug) then begin
       print,'##proj spec samp finished in ',systime(/sec)-current,' sec.'
    endif
    return, Q
end