; $Id: ref_delta.pro,v 1.6 2006/07/26 14:38:16 azuah Exp $
;###############################################################################
;
; NAME:
;  REF_DELTA
;
; PURPOSE:
;  Implementation of a Dirac delta function for use in fitting
;  data in RAINS.
;
; CATEGORY:
;  DAVE, Data Analysis, RAINS, surface fitting
;
; AUTHOR:
;   Robert M. Dimeo, Ph.D.
;   NIST Center for Neutron Research
;   100 Bureau Drive
;   Gaithersburg, MD 20899
;   Phone: (301) 975-8135
;   E-mail: robert.dimeo@nist.gov
;   http://www.ncnr.nist.gov/staff/dimeo
;
; 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 of if the code in this file is
;  included in another product.
;
;###############################################################################
; This function returns the resolution function, normalized to unit area,
; evaluated over the range of x but whose center is located at CENTER.
function ref_delta,xvals,center,resolution,_Extra = extra
if n_params() ne 3 then return,0
nx = extra.nx & ny = extra.ny
x = xvals[0:nx-1] & y = xvals[nx:nx+ny-1]

;RTA - change implementation

; !!! Old Implementation
;dx = x[1]-x[0]
;wh = where(center eq 0.0,count_zero,complement = ok,ncomplement = nok)
;nshift = 0*center
;if nok gt 0 then nshift[ok] = fix(center[ok]/dx)+abs(center[ok])/center[ok]
;;if center ne 0.0 then begin
;;  nshift = fix(center/dx)+abs(center)/center
;;endif else begin
;;  nshift = 0
;;endelse
;yout = resolution
;; Normalize the resolution function
;for i = 0,ny-1 do begin
;   area = int_tabulated(x,resolution[*,i])
;   yout[*,i] = shift(resolution[*,i]/area,nshift[0,i])
;endfor

;!!! New Implementation
; interpol() does not extrapolate well, so only interpol resolution in
; region of overlap between xRes+center and x. The convolution should
; be zero otherwise.
xRes = resolution.x
zRes = resolution.z

con = fltarr(nx,ny)             ; init to 0.0
nr = n_elements(xRes)
xr0 = xRes[0]
xrn = xRes[nr-1]
for i=0,ny-1 do begin
   ind1 =  where(x ge xr0+center[i] and x le xrn+center[i],cnt);,complement=ind2,ncomplement=ncnt)
   if (cnt gt 0) then con[ind1,i] = interpol(zRes[*,i],xRes+center[i],x[ind1])
   ;if (ncnt gt 0) then con[ind2,i] = 0.0
endfor

return,con
end
