; $Id: convolute.pro,v 1.2 2007/11/16 21:58:36 azuah Exp $
;###############################################################################
;
; NAME:
;  convolute
;
; PURPOSE:
;   Perform a numerical convolution between two 1D datasets.
;   Based on convol_as_approx written by Rob Dimeo.
;
; CATEGORY:
;   DAVE Data Analysis
;
; AUTHOR:
;   Richard Tumanjong Azuah
;   NIST Center for Neutron Research
;   100 Bureau Drive, Gaithersburg, MD 20899
;   United States
;   azuah@nist.gov; (301) 9755604
;   July, 2006
;
; 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.
;
;###############################################################################

;===============================================================================
; convolute
;
; PURPOSE:
;  Perform a numerical convolution between two 1D datasets.
;
; PARAMETERS:
;   xFun [in] - The independent variable of the first dataset or
;               intrinsic data. It is also assummed to be the
;               independent axis of the convoluted result.
;
;   yFun [in] - The dependent variable of the first dataset or
;               intrinsic data
;
;   xRes [in] - The independent variable of the second dataset or
;               resolution data
;
;   yRes [in] - The dependent variable of the second dataset or
;               resolution data
;
; RETURN VALUE:
;   - The convolution integral. The length of the result is the same
;     as the xFun parameter.
;
function convolute, xFun,yFun,xRes,yRes;,xr,r,xs,s
compile_opt idl2
;
; int[ yFun(xFun-xRes) * yRes(xRes) * delta_xRes ]
;
nxf = n_elements(xFun)
nxr = n_elements(xRes)

uXf = xFun*0.+1                   ;'unit' vector of length xFun
uXr = xRes*0.+1                   ;'unit' vector of length xRes

; yFun(xFun-xRes)
xf_xr = uXr#xFun - xRes#uXf
yf_xf_xr = interpol(yFun,xFun,xf_xr)

; interpol() does not extrapolate well so set all results for which
; xf_xr  fall outside xFun_min and xFun_max to zero.
xmin = min(xFun,max=xmax)
ind = where(xf_xr lt xmin or xf_xr gt xmax,cnt)
if (cnt gt 0) then yf_xf_xr[ind] = 0.0

; yRes(xRes)
yr_xr = yRes#uXf

; delta_xRes
dxr = xRes[1:nxr-1] - xRes[0:nxr-2]
dxr = [dxr,dxr[nxr-2]]
dxr = dxr#uXf

; integrate
conv = total(yf_xf_xr * yr_xr * dxr, 1)

return, conv 

end


;===============================================================================
; test_convolute
;
; PURPOSE:
;  Test program for convolute function
pro test_convolute
compile_opt idl2
!except = 0
nxs = 164
xlo = -10.0 & xhi = 10.0
dxs = (xhi-xlo)/(nxs-1.0)
xs = xlo+dxs*findgen(nxs)
xs = xs + 0.5*dxs*randomn(s,nxs)
xs = xs[sort(xs)]

nxr = 101
dxr = (xhi-xlo)/(nxr-1.0)
xr = xlo+dxr*findgen(nxr)
xr = xr+0.5*dxr*randomn(s,nxr)
xr = xr[sort(xr)]

delta = 5.0
h = 4.0*1.0/(3.0*delta)
rcen = 0.0
r = h*(rmd_heaviside(xr+0.5*delta-rcen)-0.5*rmd_heaviside(xr-rcen)-0.5*rmd_heaviside(xr-0.5*delta-rcen))
area = 1.0
center = 0.0
fwhm1 = 2.5
r = gauss(xr,area,center,fwhm1)

area = 1.0
center = 3.5
fwhm2 = 2.0
s = gauss(xs,area,center,fwhm2)
;s = xs*0.5 + 5.0

yhi = 1.2*max([max(r),max(s)])


con = convolute(xs,s,xr,r)
expected = gauss(xr,area,center,sqrt(fwhm1^2+fwhm2^2))

window,/free,xsize = 600,ysize = 400
wid = !d.window

device,decomposed = 1
white = fsc_color('white')
black = fsc_color('black')
red = fsc_color('red')
blue = fsc_color('blue')


plot,xs,s,yrange = [0.0,yhi],/ysty,linestyle=0,color=black,background=white $
     ,xrange = [xlo,xhi],/xsty, /nodata ,xtitle = 'x'
oplot,xs,s,psym = 0,linestyle = 2,color = black
oplot,xr,r,color = black
oplot,xr,expected,color = red, psym=2
oplot,xs,con,psym = 4,color = blue,thick = 2.0
;xyouts,-3.0,1.7,'R(x)',color = black,/data,charsize = 2.
;xyouts,0.0,0.3,'S(x)',color = black,/data,charsize = 2.0
;xyouts,-8.0,0.3,'I(x)',color = black,/data,charsize = 2.0

ch = ''
print,' '
print,' '
read,'To exit, use carriage return:',ch

wdelete, wid
end

; *********************************************************** ;
function gauss,x,area,center,fwhm
compile_opt idl2
sigma = fwhm/2.354
g = (area/sqrt(2.0*!dpi*sigma^2))*exp(-0.5*((x-center)/sigma)^2)
return,g
end

