; $Id: ref_dirac_delta.pro,v 1.3 2006/07/26 14:41:35 azuah Exp $
;###############################################################################
;
; NAME:
;  REF_DIRAC_DELTA
;
; PURPOSE:
;  Implementation of the Dirac delta function model for use in fitting
;  data in RAINS.
;
; CATEGORY:
;  DAVE, Data Analysis, RAINS, surface 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.
;
;###############################################################################
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ref_dirac_delta,xvals,parms,          $
                     parmnames = parmnames,         $
                     init_guess = init_guess,         $
                     qty = qty,               $
                     single_parmnames = single_parmnames,   $
                     multi_parmnames = multi_parmnames,   $
                     res_ptr = res_ptr,            $
                     all_y = all_y,             $
                     step = step,               $
                     xunits = xunits,            $
                     _Extra = extra
; all_y keyword is not relevant for this model
if n_elements(all_y) eq 0 then all_y = 1

nx = extra.nx & ny = extra.ny
x = xvals[0:nx-1] & y = xvals[nx:nx+ny-1]
ux = 1+bytarr(nx) & uy = 1+bytarr(ny)
parmnames = ['Intensity:'+strtrim(string(indgen(ny)),2), $
      'Center:'+strtrim(string(indgen(ny)),2)]

single_parmnames = ''
multi_parmnames = ['Intensity','Center']
if arg_present(step) then begin
   step = dblarr(2*extra.ny)
endif
if (n_params() eq 0) or (n_params() eq 1) then begin
   if arg_present(init_guess) then begin
     intensity = max(qty,dimension = 1,index)
     qty_size = size(qty)
     for i = 0,qty_size[2]-1 do begin
      qty_max = max(qty[*,i],index)
      if i eq 0 then index_array = index else $
         index_array = [index_array,index]
     endfor
     ; Guess the 1/2 width at the highest y-value
     ; Get the low value first
     amp = max(qty[*,qty_size[2]-1])
     init_guess = [intensity,x[index_array]]
   endif
   return,-1
endif

intensity = ux#(parms[0:ny-1])
;center = ux#(parms[ny:2*ny-1])
center = parms[ny:2*ny-1]       ; RTA
resolution = *res_ptr
delta_function = ref_delta(xvals,center,resolution,nx = nx,ny = ny)
yout = intensity*delta_function
bad_index = where(~finite(yout),count_bad)
if count_bad gt 0 then yout[bad_index] = 0d
return,yout
end
