; $Id: pan_convolute.pro,v 1.1 2006/07/25 20:53:00 azuah Exp $
;###############################################################################
;
; NAME:
;  PAN_CONVOLUTE
;
; PURPOSE:
;  Performs numerical convolution of two arrays.
;  Replacement for rmd_convol which only worked for data with equally
;  spaced x bins and with both the convolution kernel and the intrinsic 
;  function assumed to share the same x. This procedure makes use of 
;  of the convolute procedure which elliminates these restrictions.
;
; CATEGORY:
;  DAVE, Data Analysis, PAN, curve fitting
;
; AUTHOR:
;   Richard Tumanjong Azuah
;   NIST Center for Neutron Research
;   100 Bureau Drive
;   Gaithersburg, MD 20899
;   Jul, 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.
;
;###############################################################################
;
; CALLING SEQUENCE:
;
;       Y = pan_convolute(X,RESPTR,YFUN, $
;               RESLIMIT = RESLIMIT, $
;               PSEUDO_DELTA = PSEUDO_DELTA)
;
; INPUT PARAMETERS:
;
;   X -    A numerical vector.
;   RESPTR - pointer to the convolution kernel structure containing
;            its ind and dep values.
;   Y -   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.
;
; OUTPUT KEYWORDS:
;
; REQUIRED PROGRAMS:
;
;  convolute.pro
;
; EXAMPLE   (example uses RMD_GAUSSIAN)
;
;  IDL>   xlo = -10.0 & xhi = 10.0 & npts = 200 & dx = (xhi-xlo)/(npts-1.0)
;  IDL>   xres = xlo+dx*findgen(nx)
;  IDL>   x = xres + 20
;  IDL>   y = rmd_gaussian(x,area = 2.0,center = 1.0,width = 2.0)
;  IDL>   yres = rmd_gaussian(xres,area = 1.0,center = 0.0,width = 1.0)
;  IDL>   resPtr = ptr_new({x:xres,data:y})
;  IDL>   ycon = pan_convolute(x,resPtr,y,pseudo_delta = [[0.5,-5.0],[0.5,-5.0]])
;  IDL>   plot,x,ycon,thick = 2.0
;  IDL>   ptr_free, resPtr
;
; MODIFICATION HISTORY:
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function pan_convolute, x,resPtr,y,reslimit = reslimit, pseudo_delta = pseudo_delta
yres = (*resPtr).data
xres = (*resPtr).x

nx = n_elements(x)
ny = n_elements(y)

; If a pseudo_delta function is requires, evaluate it
yout = x*0.0                    ; init to 0.0
delta = 0
if n_elements(pseudo_delta) ne 0 then begin
   delta = 1
   delArea = pseudo_delta[0,*]
   delCen = pseudo_delta[1,*]
   n = n_elements(delArea)
   yout = fltarr(nx)
   nr = n_elements(xRes)
   xr0 = xRes[0]
   xrn = xRes[nr-1]
   for i = 0,n-1 do begin
      ind1 =  where(x ge xr0+delCen[i] and x le xrn+delCen[i],cnt) ;,complement=ind2,ncomplement=ncnt)
      if (cnt gt 0) then $
         yout[ind1] = yout[ind1] + delArea[i]*interpol(yres,xres+delCen[i],x[ind1])
   endfor

   ;; if only the delta is required, return now
   if n_params() eq 2 then return,yout
endif

if ((ny gt 1)) then begin
   conv = convolute(x,y,xres,yres) ;xres,yres,x,y)
   return, (delta)? yout+conv : conv 
endif else $
   return, yout

end
