; $Id: gaussfitklk.pro,v 1.1 1993/04/02 19:43:31 idl Exp $

PRO	GAUSSPARAB_FUNCT,X,A,F,PDER
;+
; NAME:
;	GAUSSPARAB_FUNCT
;
; PURPOSE:
;	EVALUATE THE SUM OF A GAUSSIAN AND A 2ND ORDER POLYNOMIAL
;	AND OPTIONALLY RETURN THE VALUE OF IT'S PARTIAL DERIVATIVES.
;	NORMALLY, THIS FUNCTION IS USED BY CURVEFIT TO FIT THE
;	SUM OF A LINE AND A VARYING BACKGROUND TO ACTUAL DATA.
;
; CATEGORY:
;	NRH1 POSITIONS
; CALLING SEQUENCE:
;	FUNCT,X,A,F,PDER
; INPUTS:
;	X = VALUES OF INDEPENDENT VARIABLE.
;	A = PARAMETERS OF EQUATION DESCRIBED BELOW.
; OUTPUTS:
;	F = VALUE OF FUNCTION AT EACH X(I).
;
; OPTIONAL OUTPUT PARAMETERS:
;	PDER = (N_ELEMENTS(X),6) ARRAY CONTAINING THE
;		PARTIAL DERIVATIVES.  P(I,J) = DERIVATIVE
;		AT ITH POINT W/RESPECT TO JTH PARAMETER.
; COMMON BLOCKS:
;	NONE.
; SIDE EFFECTS:
;	NONE.
; RESTRICTIONS:
;	NONE.
; PROCEDURE:
;	F = A(0)*EXP(-Z^2/2) + A(3) + A(4)*X + A(5)*X^2
;	Z = (X-A(1))/A(2)
; MODIFICATION HISTORY: (bonmartin@obspm.fr)
;    16/11/98	adapte des routines de RSI par KLK (Ludwig Klein)
;
;-
;	WRITTEN, DMS, RSI, SEPT, 1982.
;	Modified, DMS, Oct 1990.  Avoids divide by 0 if A(2) is 0.
;	Added to Gauss_fit, when the variable function name to
;		Curve_fit was implemented.  DMS, Nov, 1990.
;
;	ON_ERROR,2                        ;Return to caller if an error occurs
	if a(2) ne 0.0 then Z = (X-A(1))/A(2) $	;GET Z
	else z= 10.
;	EZ = EXP(-Z^2/2.)*(ABS(Z) LE 7.) ;GAUSSIAN PART IGNORE SMALL TERMS

	EZ = FLTARR(N_ELEMENTS(X)) & EZ = EXP(-Z^2/2. > (-24.5)) 
	nuls = WHERE( EZ LT 2.2897e-11)	; coupure a EZ=exp(-49/2)
	IF ((SIZE(nuls))(0) GT 0) THEN EZ( nuls ) = 0.


	F = A(0)*EZ + A(3) + A(4)*X + A(5)*X^2 ;FUNCTIONS.
	IF N_PARAMS(0) LE 3 THEN RETURN ;NEED PARTIAL?
;
	PDER = FLTARR(N_ELEMENTS(X),6) ;YES, MAKE ARRAY.
	PDER(0,0) = EZ		;COMPUTE PARTIALS
	if a(2) ne 0. then PDER(0,1) = A(0) * EZ * Z/A(2)
	PDER(0,2) = PDER(*,1) * Z
	PDER(*,3) = 1.
	PDER(0,4) = X
	PDER(0,5) = X^2
	RETURN
END

PRO	GAUSSLIN_FUNCT,X,A,F,PDER
;+
; NAME:
;	GAUSSLIN_FUNCT
;
; PURPOSE:
;	EVALUATE THE SUM OF A GAUSSIAN AND A STRAIGHT LINE 
;	AND OPTIONALLY RETURN THE VALUE OF ITS PARTIAL DERIVATIVES.
;	NORMALLY, THIS FUNCTION IS USED BY CURVEFIT TO FIT THE
;	SUM OF A LINE AND A VARYING BACKGROUND TO ACTUAL DATA.
;
; CATEGORY:
;	NRH1 POSITIONS
; CALLING SEQUENCE:
;	FUNCT,X,A,F,PDER
; INPUTS:
;	X = VALUES OF INDEPENDENT VARIABLE.
;	A = PARAMETERS OF EQUATION DESCRIBED BELOW.
; OUTPUTS:
;	F = VALUE OF FUNCTION AT EACH X(I).
;
; OPTIONAL OUTPUT PARAMETERS:
;	PDER = (N_ELEMENTS(X),5) ARRAY CONTAINING THE
;		PARTIAL DERIVATIVES.  P(I,J) = DERIVATIVE
;		AT ITH POINT W/RESPECT TO JTH PARAMETER.
; COMMON BLOCKS:
;	NONE.
; SIDE EFFECTS:
;	NONE.
; RESTRICTIONS:
;	NONE.
; PROCEDURE:
;	F = A(0)*EXP(-Z^2/2) + A(3) + A(4)*X + A(5)*X^2
;	Z = (X-A(1))/A(2)
; MODIFICATION HISTORY: (bonmartin@obspm.fr)
;    16/11/98	adapte des routines de RSI par KLK (Ludwig Klein)
;
;-
;	WRITTEN, DMS, RSI, SEPT, 1982.
;	Modified, DMS, Oct 1990.  Avoids divide by 0 if A(2) is 0.
;	Added to Gauss_fit, when the variable function name to
;		Curve_fit was implemented.  DMS, Nov, 1990.
;
;	ON_ERROR,2                        ;Return to caller if an error occurs
	if a(2) ne 0.0 then Z = (X-A(1))/A(2) $	;GET Z
	else z= 10.
;	EZ = EXP(-Z^2/2.)*(ABS(Z) LE 7.) ;GAUSSIAN PART IGNORE SMALL TERMS

	EZ = FLTARR(N_ELEMENTS(X)) & EZ = EXP(-Z^2/2. > (-24.5)) 
	nuls = WHERE( EZ LT 2.2897e-11)	; coupure a EZ=exp(-49/2)
	IF ((SIZE(nuls))(0) GT 0) THEN EZ( nuls ) = 0.

	F = A(0)*EZ + A(3) + A(4)*X  ;FUNCTIONS.
	IF N_PARAMS(0) LE 3 THEN RETURN ;NEED PARTIAL?
;
	PDER = FLTARR(N_ELEMENTS(X),5) ;YES, MAKE ARRAY.
	PDER(0,0) = EZ		;COMPUTE PARTIALS
	if a(2) ne 0. then PDER(0,1) = A(0) * EZ * Z/A(2)
	PDER(0,2) = PDER(*,1) * Z
	PDER(*,3) = 1.
	PDER(0,4) = X
	RETURN
END


PRO	GAUSSCST_FUNCT,X,A,F,PDER
;+
; NAME:
;	GAUSSCST_FUNCT
;
; PURPOSE:
;	EVALUATE THE SUM OF A GAUSSIAN FUNCTION AND A CONSTANT 
;	AND OPTIONALLY RETURN THE VALUE OF ITS PARTIAL DERIVATIVES.
;	NORMALLY, THIS FUNCTION IS USED BY CURVEFIT TO FIT THE
;	SUM OF A LINE AND A VARYING BACKGROUND TO ACTUAL DATA.
;
; CATEGORY:
;	NRH1 POSITIONS
; CALLING SEQUENCE:
;	FUNCT,X,A,F,PDER
; INPUTS:
;	X = VALUES OF INDEPENDENT VARIABLE.
;	A = PARAMETERS OF EQUATION DESCRIBED BELOW.
; OUTPUTS:
;	F = VALUE OF FUNCTION AT EACH X(I).
;
; OPTIONAL OUTPUT PARAMETERS:
;	PDER = (N_ELEMENTS(X),5) ARRAY CONTAINING THE
;		PARTIAL DERIVATIVES.  P(I,J) = DERIVATIVE
;		AT ITH POINT W/RESPECT TO JTH PARAMETER.
; COMMON BLOCKS:
;	NONE.
; SIDE EFFECTS:
;	NONE.
; RESTRICTIONS:
;	NONE.
; PROCEDURE:
;	F = A(0)*EXP(-Z^2/2) + A(3)
;	Z = (X-A(1))/A(2)
; MODIFICATION HISTORY: (bonmartin@obspm.fr)
;    16/11/98	adapte des routines de RSI par KLK (Ludwig Klein)
;
;-
;	WRITTEN, DMS, RSI, SEPT, 1982.
;	Modified, DMS, Oct 1990.  Avoids divide by 0 if A(2) is 0.
;	Added to Gauss_fit, when the variable function name to
;		Curve_fit was implemented.  DMS, Nov, 1990.
;
;	ON_ERROR,2                        ;Return to caller if an error occurs
	if a(2) ne 0.0 then Z = (X-A(1))/A(2) $	;GET Z
	else z= 10.
;	EZ = EXP(-Z^2/2.)*(ABS(Z) LE 7.) ;GAUSSIAN PART IGNORE SMALL TERMS

	EZ = FLTARR(N_ELEMENTS(X)) & EZ = EXP(-Z^2/2. > (-24.5)) 
	nuls = WHERE( EZ LT 2.2897e-11)	; coupure a EZ=exp(-49/2)
	IF ((SIZE(nuls))(0) GT 0) THEN EZ( nuls ) = 0.

	F = A(0)*EZ + A(3)  ;FUNCTIONS.
	IF N_PARAMS(0) LE 3 THEN RETURN ;NEED PARTIAL?
;
	PDER = FLTARR(N_ELEMENTS(X),4) ;YES, MAKE ARRAY.
	PDER(0,0) = EZ		;COMPUTE PARTIALS
	if a(2) ne 0. then PDER(0,1) = A(0) * EZ * Z/A(2)
	PDER(0,2) = PDER(*,1) * Z
	PDER(*,3) = 1.
	RETURN
END



Function Gaussfitklk, x, y, a, sigmaa, SIMPLE=simple, CONST=const, $
	DROITE=droite, PARABOLE=parabole
;+
; NAME:
;	GAUSSFIT
;
; PURPOSE:
; 	Fit the equation y=f(x) where:
;
; 		F(x) = A0*EXP(-z^2/2) + A3 + A4*x + A5*x^2
; 			and
;		z=(x-A1)/A2
;
;	A0 = height of exp, A1 = center of exp, A2 = sigma (the width).
;	A3 = constant term, A4 = linear term, A5 = quadratic term.
; 	The parameters A0, A1, A2, A3 are estimated and then CURVEFIT is 
;	called.
;	
;	Modif KLK 1/2/1996: Ajuster, selon le mot cle, une Gaussienne, sa somme
;	avec une constante, une droite ou une parabole.
;
; CATEGORY:
;	NRH1 Positions
;
; CALLING SEQUENCE:
;	Result = GAUSSFIT(X, Y [, A])
;
; INPUTS:
;	X:	The independent variable.  X must be a vector.
;	Y:	The dependent variable.  Y must have the same number of points
;		as X.
;
; OUTPUTS:
;	The fitted function is returned.
;
; OPTIONAL OUTPUT PARAMETERS:
;	A:	The coefficients of the fit.  A is a six-element vector as 
;		described under PURPOSE.
;
; COMMON BLOCKS:
;	None.
;
; SIDE EFFECTS:
;	None.
;
; RESTRICTIONS:
;	The peak or minimum of the Gaussian must be the largest
;	point in the Y vector.
;
; PROCEDURE:
;	If the (MAX-AVG) of Y is larger than (AVG-MIN) then it is assumed
;	that the line is an emission line, otherwise it is assumed there
;	is an absorbtion line.  The estimated center is the MAX or MIN
;	element.  The height is (MAX-AVG) or (AVG-MIN) respectively.
;	The width is found by searching out from the extrema until
;	a point is found less than the 1/e value.
;
; MODIFICATION HISTORY: (bonmartin@obspm.fr)
;    16/11/98	adapte des routines de RSI par KLK (Ludwig Klein)
;
;-
;	DMS, RSI, Dec, 1983.
;
; on_error,2                      ;Return to caller if an error occurs

n = n_elements(y)		;# of points.

c1 = (y(n-1) - y(0)) / (x(n-1)-x(0)) & c0 = y(0) - c1*x(0) ; straight baseline.
yd = (y-c1*x) - c0 			;difference.

ymax=max(yd, i0) & xmax=x(i0)	;x,y and subscript of extrema
; a=fltarr(6)			;coefficient vector

dy=yd(i0)			;diff between extreme and mean

del = .367879*dy		;1/e value
i=0
while ((i0+i+1) lt n) and $	;guess at 1/e width.
	((i0-i) gt 0) and $
	(abs(yd(i0+i)) gt abs(del)) and $
	(abs(yd(i0-i)) gt abs(del)) do i=i+1


; Tentative de gestion des erreurs dans CURVEFIT

	IF(KEYWORD_SET(SIMPLE)) THEN BEGIN ; Fit Gaussian
		a = [yd(i0), x(i0), abs(x(i0)-x(i0+i))] ;estimates
		return,M_CURVEFIT(x,y,replicate(1.,n),a,sigmaa, $
			function_name = "GAUSSSIM_FUNCT") ;call curvefit
Print,"Code d'erreur = ", !ERROR

	ENDIF

	IF(KEYWORD_SET(CONST)) THEN BEGIN ; Fit Gaussian+constant value
		a = [yd(i0), x(i0), abs(x(i0)-x(i0+i)), c0 + c1*x(i0)] ; est.
		return,M_CURVEFIT(x,y,replicate(1.,n),a,sigmaa, $
			function_name = "GAUSSCST_FUNCT") ;call curvefit
	ENDIF

	IF(KEYWORD_SET(DROITE)) THEN BEGIN ; Gaussian+straight line
;	Print,'Fit Gaussian + straight line'
		a = [yd(i0), x(i0), abs(x(i0)-x(i0+i)), c0, c1] ;estimates
		return,M_CURVEFIT(x,y,replicate(1.,n),a,sigmaa, $
			function_name = "GAUSSLIN_FUNCT") ;call curvefit
	ENDIF

	IF(KEYWORD_SET(PARABOLE)) THEN BEGIN ; Gaussian+parabola
;		Print,'Fit Gaussian + parabola'
		a = [yd(i0), x(i0), abs(x(i0)-x(i0+i)), c0, c1, 0.] ;estimates
		return,M_CURVEFIT(x,y,replicate(1.,n),a,sigmaa, $
			function_name = "GAUSSPARAB_FUNCT") ;call curvefit
	ENDIF

Print,"Code d'erreur = ", !ERROR

end