pro scf0,q1

; SC7 is the main processing routine of the SCF package.  It raster
; scans through the data cube and excises square boxes containing 
; the resolution elements.  The SCF values are calculated and inserted
; in the OUTDATA variable. This is the seventh version of the SCF
; program.  There remain some coding artificats of previous versions
; in the program.
; Note: this version of the SCF requires a distance mask to work!

@varld.cb	;access variable blocks
r=(res-1)/2	;convert resolution radius (odd) to subcube radius
count=0		;keeps track of successful interations
outdata=replicate(99999.0,nx0,ny0,16)	;give us our data matrix
st=systime(1)	;time the routine for performance purposes
;------------------------------
; First, the program replaces all non-finite entries with 99999.0
;------------------------------
finitecbe=fltarr(nx0,ny0,nv0)+99999.0
finsubs=where(finite(data0))
finitecbe(finsubs)=data0(finsubs)
data0=finitecbe
signal=fltarr(nx0,ny0)+99999.0
nse=fltarr(nx0,ny0)+99999.0
q=q1

v=findgen(nv0)*dv0+v0
avg_sp=total(total(data0,1),1)/(nx0*ny0)
avg_centroid=total(avg_sp*v)/total(avg_sp)
avg_sigma=sqrt(abs(total(avg_sp*(v-avg_centroid)^2))/total(avg_sp))
sigma=fitparams(*,*,2)
fitparams(*,*,2)=replicate(avg_sigma,nx0,ny0)
errmat=fitparams(*,*,5)
;------------------------------
; These loops set up the raster scan through the cube.
;------------------------------
for x=r,nx0-r-1 do begin	; begin looping through basis spectra
	print,'Now Calculating row',fix(x+1),' of',fix(nx0) 
	; keeps the user from getting bored
for y=r,ny0-r-1 do begin

;disp0=fitparams(x,y,2)		;get the line parameters from the fit
;vdisp0=sqrt(alog(2))*disp0	;parameters data
;vpeak0=fitparams(x,y,1)
vpeakpt=fitparams(x,y,6)
vdisppt=fitparams(x,y,7)
;sig0=abs(fitparams(x,y,4))
;Tmax0=fitparams(x,y,0)
;Tint0=2*sqrt(!pi)*vdisp0*gaussint(q)*Tmax0

;------------------------------
;The rejection of data is handled by loopfit.pro.  If a datum is not
;good, loopfit.pro puts a non-zero flag in the last entry of the fit
;parameters this next step recognizes this and puts in a bad data
;column in the output file.
;------------------------------
error_key=fix(errmat(x,y))
if error_key ne 0 then begin
	goto,endloop3
endif

count=count+1 ;aha! Successful spectra, let's get to work!

;------------------------------
; A subcube of the original cube containing the resolution element must
; be extracted for calculations. This section also cacluates the corner 
; of the box in pixel coordinates
;------------------------------
ymin=(y-r)*( (y-r) ge 0)
ymax=(y+r)*((y+r) le ny0-1)+(ny0-1)*((y+r) gt ny0-1)
xmin=(x-r)*( (x-r) ge 0)
xmax=(x+r)*((x+r) le nx0-1)+(nx0-1)*((x+r) gt nx0-1)
delx=xmax-xmin+1	;calcs how big the box is
dely=ymax-ymin+1
v=findgen(nv0)
; Paolo's a genius!!:
hold=min(abs(v-(vpeakpt-q*vdisppt)),v1)
hold=min(abs(v-(vpeakpt+q*vdisppt)),v2)
incube=data0(xmin:xmax,ymin:ymax,v1:v2)	
if (v1 eq v2) then goto,endloop3
spec0=transpose(data0(x,y,*)) ; get the originial spectrum
spec0=spec0(v1:v2)
subs=where(spec0 lt 0)
if (subs(0) eq -1) then goto,endloop3
;	print,'Spectrum at '+string(fix(x))+' and '+string(fix(y))+' rejected'
;	print,v1,v2
	
;endif

velpts=v2-v1+1	;number of velocity abcissae

nse(x,y)=sqrt(total(spec0(where(spec0 lt 0))^2)/total(spec0 lt 0))
signal(x,y)=sqrt(total(spec0^2)/velpts)

base_matrix=fltarr(delx,dely,velpts)
for i=0,delx-1 do begin
	for j=0,dely-1 do begin
		base_matrix(i,j,*)=spec0
	endfor
endfor
scf=1.-sqrt(total((incube-base_matrix)^2,3)/total((incube^2+base_matrix^2),3))

;------------------------------
; The routine prepares the output arrays and the mask
;------------------------------
sz=delx*dely
outd=fltarr(delx,dely,8)
normmask=mask
normmask(x-xmin,y-ymin)=0.
normmask=normmask*(fitparams(xmin:xmax,ymin:ymax,5) eq 0)
if total(normmask) eq 0 then begin
	print,x,y,'<------ Cannot correlate with nearby spectra, rejecting'
	goto,endloop3
endif
output=total(normmask*scf)/total(normmask)
outdata(x,y,7)=output
endloop3:

endfor 
endfor

outdata(*,*,8:12)=fitparams(*,*,0:4)
s2nm=signal/nse
if total(where(errmat ne 0.)) ne -1 then s2nm(where(errmat ne 0.))=99999.0
outdata(*,*,13)=s2nm
outdata(*,*,11.)=99999.0
v=findgen(nv0)*dv0+v0
T_int=total(data0,3)
weight=fltarr(nx0,ny0,nv0)
for i=0,nx0-1 do begin
        for j=0,ny0-1 do begin
                weight(i,j,*)=v
        endfor
endfor

centroid=total(weight*data0,3)/T_int
centroid_mat=fltarr(nx0,ny0,nv0)
for i=0,nv0-1 do begin        
        centroid_mat(*,*,i)=centroid
endfor
outdata(*,*,9)=sigma
fitparams(*,*,2)=sigma
outdata(*,*,10)=total(data0*(weight-centroid_mat)^2,3)/T_int

;------------------------------
; Finally, some timekeeping routines, FYI
;------------------------------
fin=systime(1)	;stop the clock
total_time=fix(fin-st)
hours=total_time/3600
minutes_and_seconds=total_time mod 3600
minutes=minutes_and_seconds/60
seconds=minutes_and_seconds mod 60
print, 'Total Analysis time:'
if (hours eq 1) then print,'1 Hour' else print,hours,' Hours'
if (minutes eq 1) then print, '1 Minute' else print, minutes,' Minutes'
if (seconds eq 1) then print, '1 Second' else print, seconds,' Seconds'
print,'Avg. Time / spectra=', (fin-st)/count ; statistics!
print,'Total Spectra Rejected=',nx0*ny0-count ;tidy up
save,file='scf.data',outdata,fitparams	;keep the results
end