; $Id: dm_calc_shielding.pro,v 1.8 2008/03/19 20:46:50 ymqiu Exp $
;#######################################################################
;
; NAME:
;  dm_calc_shielding
;
; PURPOSE:
;  this program calculates the absorption correction for equatorial scattering only
;  annular, cylindrical, and spherical geometry are implemented
;
; CATEGORY:
;  dcs_mslice
;
; 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.
;
;#######################################################################

; calculate the traveling path in an annular or cylindrical or speherical geometry
; equatorial scattering only
; input:
;   type:   0-cylinder
;           1-annulus
;           2-sphere
;   tth:    scalar or array of scattered angle tth=0 for north in DCS   (+:cw)
;   r:      radius of the sampling point
;   phi:    angle of the sampling point in cylindrical or spherical coordinate system,
;           phi=0 for north in DCS (+:ccw)  range:[0,360)
;   psi:    angle of the sampling point in spherical coordinate range:[0,180)
;   rin:    inner radius of the annulus
;   rout:   outer radius of the annulus or the radius of the cylinder or the radius of the sphere
; output:
;   l1:     incoming path, same number of elements as that of tth
;   l2:     outgoing path, same number of elements as that of tth
;   cb:     1: incoming beam is in the 1.5 central region
;           0: outside that region
;   orc:    if set, calculate the effect of oscillating radial collimator, results are saved in passrc
;   passrc: flag indicating whether the neutron will pass the radial collimator,1-pass 0-stopped
;           same number of elements as that of tth
;           assuming the sample and radial collimator are concentric and the center is at the beam center
pro dm_shielding_calcl1l2,tth,type,r,phi,psi,rin=rin,rout=rout,l1=l1,l2=l2,cb=cb,orc=orc,passrc=passrc
    ntth   = n_elements(tth)
    l2     = dblarr(ntth)
    l1     = dblarr(ntth)
    passrc = intarr(ntth)+1
    if ntth eq 1 then begin
        l2 = l2[0] & l1 = l1[0] & passrc = passrc[0]
    endif
    cb = 0

    if r gt rout then return ;outside the sample region

    tmp = phi*!dpi/180d
    if type eq 2 then begin ;sphere
        tmp1 = psi*!dpi/180d
        rp = r*sin(tmp1)
        x  = rp*cos(tmp)    ;sampling point in cartesian coordinates
        y  = rp*sin(tmp)
        z  = r*cos(tmp1)
        rp = sqrt(rout^2-z^2)
    endif else begin        ;cylinder & annulus
        x  = r*cos(tmp)     ;sampling point in cartesian coordinates
        y  = r*sin(tmp)
        rp = rout
    endelse
    x0 = -sqrt(rp^2-y^2)     ;entry intercept point in cartesian coordinates

    if abs(y) le 0.75 then cb=1 else cb=0
    l1[*] = x-x0

    if type eq 1 then $
       if (x gt 0) and (abs(y) lt rin) then $  ;intercept the inner circle
          l1[*]  =   l1[*]-2*sqrt(rin^2-y^2)

    ;rotate the cartesian coordinate system by tth degree cw
    tmp    = tth*!dpi/180d
    sintth = sin(tmp)
    costth = cos(tmp)
    ;new coordinates for the sampling point
    xp = x*costth-y*sintth
    yp = x*sintth+y*costth
    ;exit intercept point in this new coordinate system
    x1 = sqrt(rp^2-yp^2)
    l2 = x1-xp
    if type eq 1 then begin
       ind = where((xp lt 0) and (abs(yp) lt rin),count)    ;intercept the inner circle
       if count gt 0 then $
          l2[ind] = l2[ind]-2*sqrt(rin^2-yp[ind]^2)
    endif
    if keyword_set(orc) then begin
       delta     = 0.0349066       ;=2*!dtor; DCS radial collimator 2-degree separation
       rrc       = 20.0            ;inside radius of the radial collimator
       tmp       = asin(abs(yp)/rrc)
       passrc[*] = (tmp le delta)
    endif
end

; usage:
;   ei=5 & en=findgen(1000)/1000.*10-5 & tth=findgen(900)/10.+10 & info=[1,5,6,1.5,1.25]
;   factor=dm_calc_shielding(info,ei,en,tth,npoints=npoints,/interp,/debug)
; input:
;   info:           an array of the sample geometry and absorption information
;                   [type,li,la,rout,(rin),...]
;                   type:   0-cylinder
;                           1-annulus
;                           2-sphere
;   eief:           fixed Ei or Ef based on instrument geometry
;   en:             [nen] energy transfer corresponds to the time channels
;   tth:            [ndet] two-theta angle of the detectors
;   group_leader:   group leader widget id
;   interp:         if set,calculate at most 100 en and 100 tth, the others are interpolated
;   npoints:        [nr,nphi,(npsi)] number of points used in MC-calculation default[10,36,(18)]
;   orc:            if set, consider the effect of DCS oscillating radial collimator
;   instrgeom:      0 for direct geometry, 1 for inverse geometry
; return:
;   factor:         [nen,ndet] 1/A  A is the transmission coefficient equ (6.3.3.1)
function dm_calc_shielding,info,eief,en,tth,group_leader=group_leader,interp=interp,npoints=npoints,$
    orc=orc,instrgeom=instrgeom,debug=debug
    if keyword_set(debug) then  current=systime(/sec)
    if n_elements(npoints) lt 2 then npoints=[10,36]
    if info[0] eq 2 and n_elements(npoints) lt 3 then npoints=[npoints,18]
    en1  = en  & intp_en  = 0
    tth1 = tth & intp_tth = 0

    if keyword_set(interp) then begin
        nint = 100
        if n_elements(en) gt nint then begin
            emax    = max(en,min=emin)
            en1     = emin+(findgen(nint)+0.5)/nint*(emax-emin)
            intp_en = 1
        endif
        if n_elements(tth) gt nint then begin
            tmax    = max(tth,min=tmin)
            tth1    = tmin+(findgen(nint)+0.5)/nint*(tmax-tmin)
            intp_tth= 1
        endif
    endif

    nen1     = n_elements(en1)
    ndet1    = n_elements(tth1)
    li       = info[1] & la=info[2] & rout=info[3]
    e0       = 25.3d        ;meV, thermal neutron energy 2200m/s
    phi      = findgen(npoints[1])*360./npoints[1]
    area     = !pi*rout^2
    factor   = dblarr(nen1,ndet1)
    if keyword_set(instrgeom) then begin
       efei  = eief+en1
       muin  = 1d/li+1d/(la*sqrt(efei/e0))   ;1/(incoming mean free path)
       muout = 1d/li+1d/(la*sqrt(eief/e0))   ;1/(outgoing mean free path)
    endif else begin
       efei  = eief-en1
       muin  = 1d/li+1d/(la*sqrt(eief/e0))   ;1/(incoming mean free path)
       muout = 1d/li+1d/(la*sqrt(efei/e0))   ;1/(outgoing mean free path)
    endelse
    tit =  ''
    case info[0] of
        0:  begin ;cylinder
            r   = rout*(findgen(npoints[0])+0.5)/npoints[0]
            tit = ' for cylinder'
            end
        1:  begin ;annulus
            rin = info[4]
            r   = rin+(rout-rin)*(findgen(npoints[0])+0.5)/(npoints[0])
            tit = ' for annulus'
            end
        2:  begin ;sphere
            r   = rout*(findgen(npoints[0])+0.5)/npoints[0]
            psi = findgen(npoints[2])*180./npoints[2]
            tit = ' for sphere'
            end
        else:
    endcase

    vol = 0.0

    ;message bar
    tit   = 'Absorption correction'+tit
    mesg  = 'Calculating absorption correction. Please wait...'
    base  = widget_base(title=tit,/col,/base_align_center,group_leader=group_leader,tlb_frame_attr=9,map=0)
    label = widget_label(base,value=mesg,/dynamic_resize,/align_center)
    dm_center_kid,base,group_leader
    widget_control,base,/realize,/map

    ;pseudo MC calculation
    if info[0] eq 2 then begin
       for m=0L,npoints[2]-1L do begin
           sinpsi = sin(psi[m]*!dtor)
           for i=0L,npoints[0]-1L do begin
               rsq = r[i]*r[i]
               for j=0L,npoints[1]-1L do begin
                   dm_shielding_calcl1l2,tth1,info[0],r[i],phi[j],psi[m],rout=rout,l1=l1,l2=l2,cb=cb,orc=orc,passrc=passrc
                   weight = cb*0.8+0.2  ;1 for central beam, 0.2 for outer beam
                   vol    = vol+weight
                   weight = weight*passrc
                   for k=0L,nen1-1L do begin
                       if keyword_set(instrgeom) then begin
                          factor[k,*] = factor[k,*]+exp(-muin[k]*l1-muout*l2)*weight*rsq*sinpsi
                       endif else begin
                          factor[k,*] = factor[k,*]+exp(-muin*l1-muout[k]*l2)*weight*rsq*sinpsi
                       endelse
                   endfor
               endfor
           endfor
       endfor
    endif else begin
       for i=0L,npoints[0]-1L do begin
           for j=0L,npoints[1]-1L do begin
               dm_shielding_calcl1l2,tth1,info[0],r[i],phi[j],rin=rin,rout=rout,l1=l1,l2=l2,cb=cb,orc=orc,passrc=passrc
               weight = cb*0.8+0.2  ;1 for central beam, 0.2 for outer beam
               vol    = vol+weight
               weight = weight*passrc
               for k=0L,nen1-1L do begin
                   if keyword_set(instrgeom) then begin
                      factor[k,*] = factor[k,*]+exp(-muin[k]*l1-muout*l2)*weight*r[i]
                   endif else begin
                      factor[k,*] = factor[k,*]+exp(-muin*l1-muout[k]*l2)*weight*r[i]
                   endelse
               endfor
           endfor
       endfor
    endelse

    case info[0] of
        0:  begin  ;cylinder
            factor = factor*2./rout
            end
        1:  begin  ;annulus
            factor = factor*2./(rout+rin)
            end
        2:  begin  ;sphere
            factor = factor*1.5*!pi/rout^2
            end
        else:
    endcase

    if keyword_set(interp) then begin
       nen  = n_elements(en)
       ndet = n_elements(tth)
       ;first interploate over all en
       if intp_en eq 1 then begin
          factor1 = fltarr(nen,ndet1)
          for i=0L,ndet1-1 do $
              factor1[*,i] = spline(en1,factor[*,i],en)
       endif else factor1=temporary(factor)
       ;then interpolate over all tth
       if intp_tth eq 1 then begin
          ind    = sort(tth)
          tth2   = tth[ind]
          factor = fltarr(nen,ndet)
          for i=0L,nen-1 do $
              factor[i,*] = spline(tth1,factor1[i,*],tth2)
          tmp    = factor
          factor[*,ind[*]] = tmp[*,*]
       endif else factor=temporary(factor1)
    endif
    widget_control,base,/destroy
    if keyword_set(debug) then  print,'calculation finished in ',systime(/sec)-current,' secs.'
    return,reform(vol/float(factor))
end