Viewing contents of file '../idllib/iuedac/iuelib/pro/g1fit.pro'
;************************************************************************
;+
;*NAME:
;
;   	G1FIT
;  
;*CLASS:
;  
;  	Menu Procedure
;	Spectral Fitting
;  
;*CATEGORY:
;
;*PURPOSE:
;
;    	Procedure to fit gaussians and a linear baseline to data points (from 
;	GAUSSFIT.PRO)
;  
;*CALLING SEQUENCE:
;
;	G1FIT,X,Y,NCEN,CENTER,A,SIGMA
;  
;PARAMETERS:
;
;	X	(REQ) (I) (1) (I L F D)
;		Vector of independent variable elements.
;
;       Y	(REQ) (I) (1) (I L F D)
;		Vector of dependent variable elements.
;
;       NCEN	(REQ) (I) (0) (I)
;		Number of orders in signal region.
;
;       CENTER	(REQ) (I) (0) (I)
;		Indice of central order (e.g. center order number - 1).
;
;     	A	(REQ) (O) (1) (F)
;		Vector of function parameters.
;
;       SIGMA	(REQ) (O) (1) (F)
;		Sigmas from background fit.
;
;*SYSTEM VARIABLES USED:
;  
;  	NOPRINT	-   !NOPRINT
;
;*INTERACTIVE INPUT:
;  
;*SUBROUTINES CALLED:
;  
;	PARCHECK
;	G1BCKGND
;  	WFIT
;  
;*FILES USED:
;  
;*SIDE EFFECTS:
;  
;*RESTRICTIONS:  
;
;	modified for unix/sun idl version 1.1
;  
;*NOTES:
;
;       tested with IDL Version 2.1.0 (sunos sparc)    07 Aug 91
;       tested with IDL Version 2.1.0 (ultrix mispel)  07 Aug 91
;       tested with IDL Version 2.1.0 (vax vms)        07 Aug 91
;  
;*PROCEDURE:
;
;    	Procedure to fit gaussians and a linear baseline to data points (from 
;	GAUSSFIT.PRO)	
;  
;*I_HELP nn:
;  
;*MODIFICATION HISTORY:
;  
;   nov-21-89 jtb@gsfc the vax/vms version was modified for unix/sun idl
;   2 July 91 LLT cleaned up; tested on VAX; updated prolog
;   15 Jul 91 PJL added PARCHECK; tested on SUN and VAX; updated prolog
;   06 Aug 91 GRA removed unused parameter "0." fro, the WFIT calling
;                 statements.
;  
;-
;************************************************************************
 pro g1fit,x,y,ncen,center,a,sigma
;
 npar = n_params()
 if npar eq 0 then begin
    print,' G1FIT, X,Y,NCEN,CENTER,A,SIGMA'
    retall 
 endif  ; npar
 parcheck,npar,6,'G1FIT'
;
 b=fltarr(3) & deltaa=b
 a=fltarr(6)
 range = ncen/2
 lind = center - range   ; indice of right edge of signal region
 rind = center + range   ; indice of left edge of signal region
 mr = range/2            ; center region for avs calculation
 fwid = fix(ncen/9)      ; first estimate of width (1 or 2)
;
; perform linear lsf to bkgnd, extract feature and subtract bkgnd
;
 g1bckgnd,x,y,lind,rind,a,sigma
 xg = x(lind:rind)
 yb = a(3) + a(4)*xg + a(5)*xg*xg                ; calc bg vector
 yg = y(lind:rind) - yb                          ; subtract bckgnd   
;
; if average signal over source region lt sigma/2,
; store zeroes in a(0), a(1), & a(2) & return
;
 avs = total(yg(range-mr:range+mr)) / (2.0*mr+1)   ; center 2*mr+1 lines
 if avs lt sigma/2. then begin
    a([0,1,2]) = 0.0
    if !noprint eq 0 then begin
       print,'Average signal =',avs,'  sigma/2 =',sigma/2.
       print,'Average source signal too low! Return zeroes'          
    endif  ; !noprint
    return 
 endif  ; avs
; 
; estimate gaussian parameters
;
 b([0,1,2]) = [x(center),fwid,y(center)]  ; 1st est. for center, width & height
 deltaa([0,1]) = 1.0/(sqrt(ncen+1.)*3.)
 deltaa(2) = sqrt(sigma)
; 
; use wfit to fit gaussian to feature (use 1/sigma for wt)
;
 ifit = intarr(3) + 1
 wfit,xg,yg,yg*0+1/sigma,ifit,deltaa,b,ygf,sig      
;
; check for points > fac*sigma from best gaussian fit
;
 fac = 1
 absnet = abs(yg-ygf)
 repeat begin
    fac = fac + 1
    dif = absnet - fac*sigma           ; dif is excess over fac*sigma
    ind = where(dif le 0.0,count)      ; count will equal number of pts. kept
 endrep until count gt range          ; continue til pts saved > range 
 if (!noprint eq 0) then $
    print,format="(1x,i4,a,i2,a,i2,a)", $
    count,' of ',ncen,' points kept; ',fac,' sigma deviations excluded'
;
; re-fit gaussian for better parameters
;
 if count ne ncen then begin
    yg = yg(ind)
    xg = xg(ind)
    if (!noprint eq 0) then print,'refitting gaussian'
    wfit,xg,yg,yg*0+1/sigma,ifit,deltaa,b,ygf,sig      
 endif  ; count ne ncen
;
; add gaussian parameters to a
;
 b(1)=abs(b(1))
 a(0) = b 
 return
 end  ; g1fit