;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; NAME: 
;      varismooth.pro
; PURPOSE:
;      To smooth spectra with variable resolution
; INPUTS:
;      x   - array of x values (wavelengths)
;      y   - array of y values (fluxes to be smoothed)
;      r   - array of resolutions (R = lambda/FWHM)
; OUTPUT:
;      ysmooth - array of smoothed y values
; MODIFICATION HISTORY:
;     Written 09 Nov. 2001 by WDV
;
; CALLING SEQUENCE:
;     varismooth, x, y, r, ysmooth
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

pro varismooth, x, y, r, ysmooth, SILENT=silent

siglim  = 6.0
const   = 2.0*sqrt(2.0*alog(2.0))
sigma   = x/(r*const)

nx      = n_elements(x)
ysmooth = fltarr(nx)+!values.f_nan
for i=0L,long(nx)-1L do begin

; --- Generate Gaussian array

    if finite(sigma[i]) eq 0 then goto, cont
    ggaus = fltarr(nx)
    wgaus = x - x[i]
    zg    = where(abs(wgaus/sigma[i]) le siglim)
    ggaus = (exp(-0.5*((wgaus[zg]/sigma[i])^2)))/(sigma[i]*sqrt(2.0*!pi))
    scale = total(ggaus)
    ggaus = temporary(ggaus)/scale
    if (n_elements(zg) gt nx/2) then begin
       print, 'Warning: kernel too large; only part of the input array will be correctly convolved'
       ysmooth[i] = !values.f_nan
    endif else begin

; --- Pad input array

       addon = fltarr(n_elements(zg)/2)
       yfull = [y,addon]

; --- Convolve spectrum with Gaussian

       if not(keyword_set(SILENT)) then print, 'Convolving...', i, ' Resolution...',r[i],' Sigma...',sigma[i]
       yconv      = convol(yfull,ggaus,total(ggaus))
       ysmooth[i] = yconv[i]
   endelse

   cont:

endfor

end