; $Id: rmd_convol.pro,v 1.4 2005/02/08 17:57:15 dimeo Exp $
;###############################################################################
;
; NAME:
;  RMD_CONVOL
;
; PURPOSE:
;  Performs numerical convolution of two arrays.
;
; CATEGORY:
;  DAVE, Data Analysis, PAN, curve 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.
;
;###############################################################################
;+
; NAME:
;       RMD_CONVOL
;
; 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.
;
;   This routine is a modification on the earlier one named RMD_CONVOLUTION.
;   The difference here is that both a sequence of Dirac delta functions
;   and an array with which the resolution function is to be convoluted is
;   specified.  See examples below for syntax.
;
;
; 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_CONVOL(X,RESPTR,YFUN, $
;               RESLIMIT = RESLIMIT, $
;               PSEUDO_DELTA = PSEUDO_DELTA, $
;               NONORM = NONORM)
;
; INPUT PARAMETERS:
;
;       X -    A numerical vector.
;   RESPTR - pointer to 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],[0.5,-5.0],[0.5,5.0]]
;   The example above allows the user to specify three Dirac delta functions
;   located at 0.0, -5.0, and 5.0 with areas of 1.0, 0.5, and 0.5, respectively.
;
;  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)
;
; EXAMPLE   (example uses RMD_GAUSSIAN)
;
;  IDL>   xlo = -10.0 & xhi = 10.0 & npts = 200 & dx = (xhi-xlo)/(npts-1.0)
;  IDL>   x = xlo+dx*findgen(nx)
;  IDL>   y = rmd_gaussian(x,area = 2.0,center = 1.0,width = 2.0)
;  IDL>   yres = rmd_gaussian(x,area = 1.0,center = 0.0,width = 1.0
;  IDL>   ycon = rmd_convol(x,yres,y,pseudo_delta = [[0.5,-5.0],[0.5,-5.0]])
;  IDL>   plot,x,ycon,thick = 2.0
;
; MODIFICATION HISTORY:
;
;       Written by Rob Dimeo, August 11, 2003.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function rmd_convol,x,resPtr,yfunn,$
            reslimit = reslimit,$
            pseudo_delta = pseudo_delta
            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 = *resPtr
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]

;; Locate the peak and fit a gaussian to it
;npeaks = 1
;xpeaks = get_peak_pos(x,yres,npeaks,indices = indices)
;a1 = xpeaks & a2 = 10.0*(x[1]-x[0]) & a0 = yres[indices[0]]
;init_guess = [a0,a1,a2]
;fit = gaussfit(x,yres,parms,estimates = init_guess,nterms = 3)
;scale = sqrt(2.0*!dpi)*parms[2]*parms[0]

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
  delArea = pseudo_delta[0,*]
  delCen = pseudo_delta[1,*]
  n = n_elements(delArea)
  yout = fltarr(n_elements(x))
  for i = 0,n-1 do yout = yout+delArea[i]*INTERPOL(yres,x+delCen[i],x)
  if n_params() eq 2 then return,yout
endif

if n_elements(interp) eq 0 then interp = 0

  dx = x[1] - x[0]
  if n_elements(reslimit) eq 0 then reslimit = 0.95*min(abs([min(x),max(x)]))
  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.
  ycon = 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 ycon = shift(ycon,circShift)
   if n_elements(yout) ne 0 and n_elements(ycon) ne 0 then return,yout+ycon
return,ycon
end