;+ ***********************************************************************
; NAME:
;	GAUSS_RH
;
; PURPOSE:
;     Ajustement d'un profil entre les minimas de part et d'autre d'un
;       point d'indice imax par une gaussienne ou la somme d'une
;	gaussienne et d'une constante, une droite ou une parabole.
;
; CATEGORY:
;	NRH1 Positions
;
; CALLING SEQUENCE:
;	GAUSS_RH, xtab, ytab, imax, imin0, imin1, coeff, sigmaa, gsfctn, $
;	SIMPLE=simple, CONST=const, DROITE=droite, PARABOLE=parabole
;
; INPUTS:
;	xtab	Abscissa array (input)
;	ytab	Ordinate array (input)
;	imax	Index of the point around which local minima are searched for
;	
; KEYWORD PARAMETERS:
;		SIMPLE	Fit Gaussian (default value)
;		CONST	Fit Gaussian + constant
;		DROITE	Fit Gaussian + straight line
;		PARABOLE	Fit Gaussian + parabola
;
; OUTPUTS:
;	imax	Index of the highest ordinate value
;	imin0	Index of the minimum below point imax
;	imin1	Index of the minimum above point imax
;	coeff	Coefficients of the fit (6 elem. array)
;	sigmaa	Vector of standard deviations for the elements of coeff
;	gsfctn	Array of computed Gaussian values
;
; COMMON BLOCKS:
;	Non
;
; MODIFICATION HISTORY: (bonmartin@obspm.fr)
;	18/11/98 Adapte du logiciel XHELIO (KLK)
;       12/04/99 Redefinition de la procedure de recherche: chercher les
;                minima locaux autour du point no. imax dont l'exces par
;                rapport a la moyenne glissante est maximum (routines
;                SIGRMN.PRO, RUNRMN_RH.PRO);
;                ajuster une fonction en s'appuyant sur un intervalle
;                symetrique par rapport a ce point
;-*******************************************************************

pro gauss_rh, xtab, ytab, imax, imin0, imin1, coeff, sigmaa, gsfctn, $
	SIMPLE=simple, CONST=const, DROITE=droite, PARABOLE=parabole

; Arguments:

npts = N_ELEMENTS(ytab)

; 1) Index of the local minimum before point IMAX
;    Search the abscissa of that amplitude which is greater than its
;    predecessor (smaller abscissa) and smaller than the value at index imax,
;    supposed to be that point where the excess above the running mean
;    is maximum (cf. routine RUNMN_RH.PRO)

i = imax & i1 = imax-1
IF(i1 LT 0) THEN BEGIN
	imin0 = 0
ENDIF ELSE BEGIN
	while(ytab(i1) le ytab(i) OR ytab(i1) GE ytab(imax) ) do begin
		if(i le 0) then i = npts -1 + i else i = i - 1
		case 1 of
			i gt 0: i1 = i - 1
			i eq 0:	i1 = npts-1
		endcase
	endwhile
	imin0 = i
ENDELSE

; 2) Index of the local minimum after point IMAX

i = imax & i1 = imax+1
IF(i1 GE npts) THEN BEGIN
	imin1 = npts-1
ENDIF ELSE BEGIN
	while(ytab(i1) le ytab(i) OR ytab(i) GE ytab(imax)) do begin
		if(i ge npts-1) then i = i - npts + 1 else i = i + 1
		case 1 of
			i lt npts-1: i1 = i + 1
			i eq npts-1:	i1 = 0
		endcase
	endwhile
	imin1 = i
ENDELSE

; Symetriser l'intervalle de canaux a ajuster par rapport au point imax
aux = MIN([imax-imin0, imin1-imax])
imin0 = imax-aux
imin1 = imax+aux

case 1 of
	imin1 ge imin0: nptsfit = imin1-imin0+1
	imin1 lt imin0: nptsfit = imin1-imin0+1+npts
endcase

gsfctn = fltarr(npts) & gsfctn(*) = 0.

if(nptsfit lt 6) then begin
	coeff =fltarr(3) & coeff(2) = 0
endif else begin

; Construction of the subarray to be fit by a Gaussian + parabola 
; (IDL fct. GAUSSFIT)

	tabx = fltarr(nptsfit) & taby = fltarr(nptsfit)
	if(imin1 gt imin0) then begin
		tabx = xtab(imin0:imin1)
		taby = ytab(imin0:imin1)
	endif else begin
		tabx(0:npts-imin0-1) = xtab(imin0:npts-1) - npts
		tabx(npts-imin0:nptsfit-1) = xtab(0:imin1)
		taby(0:npts-imin0-1) = ytab(imin0:npts-1) 
		taby(npts-imin0:nptsfit-1) = ytab(0:imin1)
	endelse

	gauss = gaussfitklk(tabx, taby, coeff, sigmaa, $
			SIMPLE=simple, CONST=const, $
			DROITE=droite, PARABOLE=parabole)

; Computation of the fitted Gaussian function (gsfctn) in a field of npts 
; points centred around the point closest to the peak of the Gaussian:
	tabx = FINDGEN(npts) + (imax - npts/2)
	pas = xtab(1) - xtab(0) & canaux = xtab(0) + tabx*pas
	IF(coeff(2) LE 0.) THEN BEGIN
		gsfctn(*) = 0.		; zero gsfctn if half width 0 or <0
	ENDIF ELSE BEGIN
		FOR i = 0, npts-1 DO BEGIN
			CASE 1 OF
				tabx(i) LT 0: ipos = tabx(i) + npts
				tabx(i) GE 0 and tabx(i) lt npts: ipos=tabx(i)
				tabx(i) GE npts: ipos = tabx(i) - npts
			ENDCASE
			z = ( canaux(i) - coeff(1) ) / coeff(2)
			aux = ALOG(coeff(0)) - .5*(z*z)
			IF(aux GT -87) THEN aux = exp(aux) ELSE aux = 0.
			gsfctn(ipos) = aux
		ENDFOR
	ENDELSE


endelse


end