Viewing contents of file '../idllib/contrib/harris/gaussian_fit.pro'
FUNCTION GAUSSIAN_FIT,X,Y,YFIT,YBAND,WEIGHT = WEIGHT
;+
; NAME:
;	GAUSSIAN_FIT
; PURPOSE:
;	Fits a gaussian using a least square polynomial fit (poly_fit) 
;	to the log of the data.
; CATEGORY:
;	?? - CURVE FITTING.
; CALLING SEQUENCE:
;	COEFF = GAUSSIAN_FIT(X,Y[,YFIT,YBAND,WEIGHT = WEIGHT] )
; INPUTS:
;	X = independent variable vector.
;	Y = dependent variable vector, should be same length as x.
;
; OUTPUTS:
;	Function result= Coefficient vector, length 2.
;	first coeff is the 1/e width,
;	second coeff is the x value of the peak
;	third coeff is the peak y value
;   OPTIONAL PARAMETERS:
;	YFIT = Vector of calculated Y's.  Has error + or - YBAND.
;	YBAND = Error estimate for each point = 1 sigma (untested)
;	WEIGHT = weighting of data for fitting
;
; COMMON BLOCKS:
;	none.
; SIDE EFFECTS:
;	none.
; MODIFICATION HISTORY:
;	Written by: Trevor Harris, Physics Dept., University of Adelaide,
;		July, 1990.
;
;-
;	ON_ERROR,2		;RETURN TO CALLER IF ERROR
	XX = X*1.		;BE SURE X IS FLOATING OR DOUBLE
	N = N_ELEMENTS(X) 	;SIZE
	IF N NE N_ELEMENTS(Y) THEN $
	  message,'X and Y must have same # of elements'
;
	case (n_params(0)) of
		0	: message,'X and Y must be given'
		1	: message,'X and Y must be given'
		else	: ;continue
	endcase
	
	coeff = fltarr(3)
	yvar  = fltarr(3)
	YFIT = fltarr(N)
	YBAND = fltarr(N)

	coeff = [0.0,0.0,0.0]
	goodpts = where(Y gt 0, count)

	if (count gt 0) then begin

		ylimit = min(Y(goodpts))/10.  
		yy = alog(Y > ylimit)

		if (not keyword_set(WEIGHT)) then begin
			ymax = max(Y,maxi)
			xwt = 1/(1+(X - X(maxi))^2*10.)
			;WEIGHT = Y > ylimit
			;WEIGHT = WEIGHT*xwt
			WEIGHT = xwt
		endif
	
		polycoeff = svdfit(X,yy,3,weight=WEIGHT,yfit=YFIT,var=yvar)

		coeff(0) = 1/sqrt(-polycoeff(2))
		coeff(1) = -0.5*polycoeff(1)/polycoeff(2)
		
		YFIT = exp(YFIT)
		tmp = exp( -(X-coeff(1))^2 / coeff(0)^2 )
		coeff(2) = median (YFIT/tmp)
		YBAND = (sqrt(yvar)#[X*X,X,1])*YFIT
	
	endif 

	return,coeff
END