;+
; NAME:
;       RMD_CONVOLUTION
;
; PURPOSE:
;
;       Implementation of a one-dimensional convolution integral.  This is
;	essentially a wrapper for IDL's built-in routine, CONVOL.PRO.  However
;	the behavior of CONVOL in which the autocorrelation function is actually
;	calculated rather than the mathematical convolution has been corrected
;	here.
;
;	In addition, convolution of a kernel with a resolution function has been
;	implemented via the PSEUDO_DELTA keyword in which the user can specify
;	the integrated area and center of the delta function prior to
;	convolution.
;
;
; 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
;
; CATEGORY:
;
;       Mathematics
;
; CALLING SEQUENCE:
;
;       Y = RMD_CONVOLUTION(X,YRES,YFUN, $
;					RESLIMIT = RESLIMIT, $
;					PSEUDO_DELTA = PSEUDO_DELTA, $
;					NONORM = NONORM)
;
; INPUT PARAMETERS:
;
;       X -	A numerical vector.
;	YRES - 	the convolution kernel
;	YFUN -	function to convolve with the kernel
;
; RETURNS:
;
;   The convolution product of the function with the kernel.
;
; INPUT KEYWORDS:
;
;	RESLIMIT	-positive limit that defines the truncation of the
;			 resolution function, i.e. the resolution function
;			 will be defined over [-RESLIMIT,+RESLIMIT].  This
;			 is useful to set if the resolution function is zero
;			 for a range of X values.
;
;	PSEUDO_DELTA	-parameters of the pseudo-delta function
;
;		Example: PSEUDO_DELTA = [1.0,0.0]
;
;	NONORM			-if set then the resolution kernel is NOT normalized to one
;				 The default (i.e. if NONORM is not set) is to normalize
;				 the resolution kernel to unit integrated area.
;
; OUTPUT KEYWORDS:
;
;
; REQUIRED PROGRAMS:
;
;        CONVOL.PRO (in the IDL distribution)
;	 RMD_HEAVISIDE.PRO
;
; EXAMPLE	(example uses RMD_MAKEPOINTS, RMD_HEAVISIDE, and RMD_GAUSSIAN)
;
;	IDL>	xlo = -5.0 & xhi = 5.0 & npts = 200
;	IDL>	x = rmd_makepoints(xlo = xlo,xhi = xhi, npts = npts)
;	IDL>	y = rmd_heaviside(x+2.0)-rmd_heaviside(x-1.0)
;	IDL>	yres = rmd_gaussian(x,area = 1.0,center = 0.0,width = 1.0
;	IDL>	ycon = rmd_convolution(x,yres,y)
;	IDL>	plot,x,ycon,yrange = [-2.0,2.0],ystyle = 1
;
;	DISCLAIMER
;
;		This software is provided as is without any warranty whatsoever.
;		Permission to use, copy, modify, and distribute modified or
;		unmodified copies is granted, provided this disclaimer
;		is included unchanged.
;
; MODIFICATION HISTORY:
;
;       Written by Rob Dimeo, September 6, 2002.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function rmd_convolution,x,yress,yfunn,$
						reslimit = reslimit,$
						pseudo_delta = pseudo_delta,$
						xdiff = xdiff,$
						nonorm = nonorm

; The next two lines are necessary because we will be modifying
; the parameters yres and yfun.  We don't want the calling program
; to see the truncated version of the resolution function!
yres = yress
if n_params() eq 3 then yfun = yfunn
; Normalize the resolution function if not already done
xsort = sort(x) & x = x[xsort] & yres = yres[xsort]
scale = int_tabulated(x,yres)
if n_elements(nonorm) eq 0 then yres = yres/scale
if n_params() eq 3 then yfun = yfun[xsort]

et = 1
ctr = 1
npts = n_elements(x)
nx = n_elements(x)
; Test if we are using a pseudo_delta function for YFUN...
if n_elements(pseudo_delta) ne 0 then begin
  if n_elements(reslimit) eq 0 then begin
    mpts = npts
    xres = x
    reslimit = max(x)
  endif
  wherevalid = where((x ge -reslimit) and (x le reslimit),mpts)
  xres = x[wherevalid]
  yres = yres[wherevalid]
  delArea = pseudo_delta[0]
  delCen = pseudo_delta[1]
  dx = x[1]-x[0]
  ; vectorized method for calculating xdiff
  nx1 = nx
  nx2 = mpts
  u1 = 1+bytarr(nx1)
  u2 = 1+bytarr(nx2)
  xdiff = transpose(u1#x[wherevalid]-x#u2)
  yout = delArea*transpose(yres)#(rmd_heaviside(xdiff+0.5*dx+delCen)-$
         rmd_heaviside(xdiff-0.5*dx+delCen))
  return,yout
endif

; First make sure that the functions are sorted properly
;xsort = sort(x)
;x = x[xsort] & yres = yres[xsort] & yfun = yfun[xsort]
if n_elements(interp) eq 0 then interp = 0
;if npts mod 2 eq 0 then interp = 1

  dx = x[1] - x[0]
  if n_elements(reslimit) eq 0 then begin
    n = 3 & nlo = n & nhi = npts - (n+1)		; symmetric cut points for kernel
  endif else begin
    wherevalid = where((x ge -reslimit) and (x le reslimit),count)
    nlo = wherevalid[0] & nhi = wherevalid[count-1]
    if count eq 0 then begin
      n = 3 & nlo = n & nhi = npts-(n+1)
    endif
  endelse
;
; Now we use the IDL built-in routine CONVOL.PRO but with
; few additions to make it work properly.  First, the kernel is
; cut symmetrically at the edges so that the number of points in
; the kernel is less than in the function.  Next, we reverse the
; the kernel so we are doing a convolution, not an autocorrelation.
; Finally we select appropriate keywords and multiply by dx so that
; the convolution sum becomes a convolution integral.

  yout = convol(yfun,reverse(yres[nlo:nhi]),edge_truncate = et,center=ctr)*dx
; Wait!...now calculate the amount by which we need to shift the result
; in case the interval of x is asymmetric....
  midpoint = min(x) + (max(x)-min(x))/2.0
  circShift = round(midpoint/dx)
; ...and return the shifted result.
  if n_elements(reslimit) eq 0 then yout = shift(yout,circShift)

return,yout
end