; $Id: hfbs_convol.pro,v 1.5 2005/02/09 16:58:01 dimeo Exp $
;###############################################################################
;
; NAME:
;  HFBS_CONVOL
;
; PURPOSE:
;  Early version of the convolution routine used by the newer data analysis
;  modules.
;
; CATEGORY:
;  DAVE, HFBS, analysis
;
; 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.
;
;###############################################################################
; hfbs_convol.pro
; Written by R.M.Dimeo (11/11/01)
;
; NAME:
;       HFBS_CONVOL
;
; PURPOSE:
;       This function calculates the proper convolution of
;       two (one-dimensional) functions that have been defined
;       over the same domain. YRES is automatically normalized
;		to unit area within the routine.
;
;		This routine corrects the fact that the built-in IDL CONVOL.PRO
;		function actually calculates the autocorrelation function of the two
;		functions.
;
;		The use of a PSEUDO_DELTA function requires that YRES be present but
;		YFUN need not be present.
;
;
; CATEGORY:
;       Signal analysis
;
; CALLING SEQUENCE:
;       Result = HFBS_CONVOL(X,YRES,YFUN, $
;							RESLIMIT = RESLIMIT, $
;							PSEUDO_DELTA = PSEUDO_DELTA,$
;							XDIFF = XDIFF)
;
; INPUTS:
;       X:			Domain of YRES and YFUN that must include zero within
;					 it's bounds.  Note that the domain does not have to
;					 be symmetric nor does it have to be equally-spaced.
;
;       YRES:	Resolution function or kernel
;
;		YFUN:	Function to be convoluted with the kernel
;
; KEYWORDS:
;
;		NONORM:	  If keyword is set then the resolution kernel is not normalized
;
;
;		RESLIMIT:	If keyword set, then the resolution function is truncated
;					so that it is only defined over the range from -RESLIMIT to
;					RESLIMIT.
;
;		PSEUDO_DELTA:
;					If keyword set, then the resolution function is convolved analytically
;					with a narrow gaussian.  A boxcar approximation to the resolution function
;					is used as the continuous representation of the actual sampled function.
;					pseudo_delta is a 3-element vector composed of the delta-function area,center
;					and width.  Width is a silly parameter for a delta function but it is used
;					because of the parametrization of the delta function as a narrow Gaussian.
; EXAMPLES:
;
;
; PROCEDURE:
;
; REFERENCE:
;       CONTINUOUS AND DISCRETE SIGNALS AND SYSTEMS
;       S.S.Soliman and M.D.Srinath
;       ISBN 0-13-171257-8
;
; MODIFICATION HISTORY:
;
;	11/30/01 (RMD): Added reslimit and pseudo_delta keywords
;	12/01/01 (RMD): Added array multiplication to speed up calculation
;								  involving a pseudo_delta function
;	12/03/01 (RMD):	Improved speed with which the pseudo-delta function
;					convolution is computed.
;	01/03/02 (RMD):	Added nonorm keyword for resolution function
;	01/28/02 (RMD): Removed the interpolation option and the double option
;					in order to gain speed.
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function hfbs_convol,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
nsort = sort(x)
scale = int_tabulated(x[nsort],yres[nsort])
if n_elements(nonorm) eq 0 then yres = yres/scale

xsort = sort(x) & x = x[xsort] & yres = yres[xsort]
if n_params() eq 3 then yfun = yfun[xsort]

et = 1
ctr = 1
npts = 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
  endif else begin
     wherevalid = where((x ge -reslimit) and (x le reslimit),mpts)
     xres = fltarr(mpts)
     xres = x[wherevalid]
     yres = yres[wherevalid]
  endelse
  delArea = pseudo_delta[0]
  delCen = pseudo_delta[1]
  delWid = pseudo_delta[2]
  delSig = delWid/2.354
  dx = x[1]-x[0]
  if n_elements(xdiff) eq 0 then begin
	; vectorized method for calculating xdiff
	nx1 = nx
	nx2 = mpts
	xdiff = fltarr(nx1,nx2)
	u1 = 1+bytarr(nx1)
	u2 = 1+bytarr(nx2)
	xdiff = transpose(u1#x[wherevalid]-x#u2)
;	loop method for calculating xdiff
;    xdiff = fltarr(mpts,npts)
;    for i = 0,mpts-1 do begin
;     	xdiff[i,0:npts-1] = (shift(xres,mpts))[i]-x[*]
;    endfor
  endif
  arg1 = dblarr(mpts,npts) & arg2 = dblarr(mpts,npts)
  arg1 = ((xdiff)+0.5*dx+delCen)/sqrt(2.0*delSig^2)
  arg2 = ((xdiff)-0.5*dx+delCen)/sqrt(2.0*delSig^2)
  yout = 0.5*delArea*transpose(yres)#(errorf(arg1)-errorf(arg2))
  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 auto-correlation.
; 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