;$Id: find_peak.pro,v 1.2 2005/02/08 18:30:58 dimeo Exp $
;###############################################################################
;
; NAME:
;  FIND_PEAK
;
; PURPOSE:
;  Simple peak location estimator using MPFIT routines.
;
; 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 find_peak,x,z,dz,npeaks
; Find the peak of the first curve
; Find the maximum data point
i = 0
z_in = z & x_in = x
xmax = max(x,min = xmin)
fwhm = 0.1*(xmax-xmin)
while i lt npeaks do begin
	z_max = max(z_in,max_index,min = z_min)
	; Fit a gaussian to the curve
	center = x_in[max_index] & sigma = fwhm/2.354
	parms = [z_max,center,sigma,z_min]
	p = 0
	yfit = mpfitpeak(x_in, z_in, p, /gaussian,nterms = 4,status = status)
	if i eq 0 then x_array = p[1] else x_array = [x_array,p[1]]
	;result = gaussfit(x,z_in,parms,nterms = 4)
	; Mask out the region encompassed by the 2*FWHM
	lo = p[1]-1.5*2.354*p[2] & hi = p[1]+1.0*2.354*p[2]
	mask = where((x_in le hi) and (x_in ge lo) ,count)
;	window,/free
;	plot,x_in,z_in,title = string(i)
;	plots,[lo,lo],!y.crange,linestyle = 2
;	plots,[hi,hi],!y.crange,linestyle = 2
	if count gt 0 then z_in[mask] = z_min
	i = i + 1
endwhile
;window,/free
;plot,x_in,z_in,psym = 0
return,x_array
end