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