;-------------------------------------
; The GAUSSIANF.PRO file contains three programs.
;	LFIT - Prepares and fits a Gaussian to the spectra
;	LOOPFIT - Raster scans through the cube and returns the 
;	fit parameters. Also handles rejection and noise calculations. 
;	GAUSSI - The target function for the fitting routine.
;-------------------------------------
function lfit,v0,dv0,velp,inspec,sv

;a leaner, meaner version of lfit.pro
g0=max(inspec,maxsubs)
g1=velp(maxsubs)
del=g0/exp(1.)	;1/e value
;this next part of the routine is from gaussfit.pro
i=0
while ((maxsubs+i+1) lt sv) and $    
        ((maxsubs-i) gt 0) and $
        (abs(inspec(maxsubs+i)) gt abs(del)) and $
        (abs(inspec(maxsubs-i)) gt abs(del)) do i=i+1
;g2=abs(velp(maxsubs)-velp(maxsubs+i))
g2=i*dv0
g2=g2+1-(g2 GT 0)
D=[g0,g1,g2,0]			;assume constant displacement ~=0
D=double(D)
Wt=fltarr(sv)+1.	;take uniform curve weighting
outp=stgaussfit(velp,inspec,Wt,D,function_name='gaussi')
return,D
end
;------------------------------------
function loopfit,incube,v0,dv0,q,s2n

sz=size(incube)		;initial dimensions
sx=sz(1)		
sy=sz(2)
sv=sz(3)
outcube=fltarr(sx,sy,9)+1.
fctn=fltarr(sv)
velp=findgen(sv)*dv0+v0
indexes=where(finite(incube) ne 1)
if total(indexes) gt -1 then incube(indexes)=99999.0


for i=0,sx-1 do begin
for j=0,sy-1 do begin
inspec=(transpose(incube(i,j,*)))
rej=0.
if max(inspec) eq min(inspec) then begin
outcube(i,j,*)=[replicate(99999.0,5),1,replicate(99999.0,3)]
print,i,j,1
goto,nodat
endif
params=lfit(v0,dv0,velp,inspec,sv)
params(2)=abs(params(2))
vdisp=sqrt(alog(2))*params(2)
vpeakpt=round(((params(1)-v0)/dv0))
vdisppt=(ceil(abs(vdisp/dv0)))

siggy=abs(noise(inspec,vpeakpt,vdisppt,4,params(3)))

;if s2n*siggy gt params(0) then rej=1
if round(vpeakpt+2*q*vdisppt) gt (sv-1) or $
round(vpeakpt-2*q*vdisppt) lt 0. then rej=2

rej=rej+(s2n*siggy gt params(0))*4+(siggy eq 99999.0)*8
outcube(i,j,*)=[params,siggy,rej,vpeakpt,vdisppt,params(0)/siggy]
print,i,j,fix(rej)
;gaussi,velp,params,fctn
;plot,vel,fctn,xsty=5,ysty=5

;oplot,vel,inspec
nodat:
endfor
endfor
print,round(total(outcube(*,*,5) ne 0)), ' Spectra Rejected of'
print,sx*sy,' Total Spectra'
return,outcube
end
;-------------------------------------
pro gaussi,X,A,F,pder
z=(x-a(1))/a(2)
bx=exp(-z^2)
F=A(3)+A(0)*bx
sub1=2*a(0)/a(2)*z*bx
sub2=2*a(0)/a(2)*z^2*bx
if N_PARAMS() GE 4 THEN $
pder=[[bx],[sub1],[sub2],[replicate(1.0,N_ELEMENTS(X))]]
end
;--------------------------------------