pro sc7

; 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

;------------------------------
; Next, if desired, the program smooths the spectra using the IDL
; routine SMOOTH
;------------------------------
if smw le 0 then goto,nosmooth

print,'Smoothing Data...'
print,'Width (in FWHMs) =',smw
data1=fltarr(nx0,ny0,nv0)

for a=0,nx0-1 do begin ; this loop smooths the spectra as appropriate 
for b=0,ny0-1 do begin	; 
vdisppt=fitparams(a,b,7)
wdth=(vdisppt*smw)
wdth=wdth*(wdth lt nv0)*(wdth ge 2)+2*(wdth ge nv0)+2*(wdth lt 2)
data1(a,b,*)=smooth(transpose(data0(a,b,*)),wdth,/edge_truncate)
endfor
endfor

data0=data1
print,'Done Smoothing, Begin Iterrations...'

nosmooth:

;------------------------------
; 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.
;------------------------------
if fitparams(x,y,5) ne 0 then begin
print,'Spectrum at',fix(x),fix(y),' rejected (Check fitparams data) '
outdata(x,y,*)=fltarr(16)+99999.0
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

incube=data0(xmin:xmax,ymin:ymax,*)	

spec0=transpose(data0(x,y,*)) ; get the originial spectrum
velpts=n_elements(spec0)	;number of velocity abcissae

;------------------------------
; The routine now fills a 4-D array with the central portions of
; each spectrum, shifting the lag as appropriate.  The subcubes
; with different shift values have different indexes in the 
; fourth dimension.
;------------------------------
stillmatrix=fltarr(delx,dely,2*q*vdisppt+1,2*q*vdisppt+1)
shiftmatrix=fltarr(delx,dely,2*q*vdisppt+1,2*q*vdisppt+1)

for i=0,2*q*vdisppt do begin
smw=q*vdisppt*(smw gt q*vdisppt)+smw*(smw le q*vdisppt)
shiftmatrix(*,*,*,i)=shiftcube(incube,q,vpeakpt,vdisppt,i-q*vdisppt,res,v0,$
	dv0);,inparams)
stillmatrix(*,*,*,i)=stillcube(spec0,q,vpeakpt,vdisppt,res,delx,dely)
endfor


;------------------------------
;Part of the Seljak Correction, see below
;------------------------------
;N1=replicate(fitparams(x,y,4)*fitparams(x,y,4),delx,dely,2*q*vdisppt+1)
;N2=fltarr(delx,dely,2*q*vdisppt+1)
;for i=0,2*q*vdisppt do N2(*,*,i)=fitparams(xmin:xmax,ymin:ymax,4)^2

;------------------------------
; Here the calculation of the scaling factor is performed, using
; the minimization methods discussed in the writeup.
;------------------------------

stau=total(shiftmatrix*stillmatrix,3)/$
(total(shiftmatrix*shiftmatrix,3));-(2*q*vdisppt+1)*N2)
;calculates s such that del^2 will be minimized

;------------------------------
; This portion calculates the correction factor due to noise proposed
; by U. Seljak.  Not implemented
;------------------------------
;uros=fltarr(delx,dely,2*q*vdisppt+1)
;uros=uros+(2*q*vdisppt+1)*(stau*stau*N1+N2)
;urosns=uros+(2*q*vdisppt+1)*(N1+N2)
uros=0
urosns=0

;------------------------------
; The calculation of the SCF proper begins now.  The function here
; is actually the deviation squared (See writeup).
;------------------------------
scf=(shiftmatrix-stillmatrix)*(shiftmatrix-stillmatrix)

staumult=fltarr(delx,dely,2*q*vdisppt+1,2*q*vdisppt+1) ; initializes s matrix
for i=0,2*q*vdisppt do staumult(*,*,i,*)=stau

scfmin=total((staumult*shiftmatrix-stillmatrix)*$
(staumult*shiftmatrix-stillmatrix),3) 
;calculates minimum
;value of del^2 for each value of tau with s chosen to minimize

norm=stau*stau*total(shiftmatrix*shiftmatrix,3)+$
total(stillmatrix*stillmatrix,3)
;calculates the maximum value of the deviation in the case of
;no absorption in the spectra.

scfmin=(scfmin-uros)/(norm-uros)
;application of the uros correction if implemented.


noscale=total(scf,3) ; takes s=1 minimum del^2
noscale=(noscale-urosns)/(total(shiftmatrix*shiftmatrix,3)$
+total(stillmatrix*stillmatrix,3)-urosns)

;------------------------------
; The routine prepares the output arrays and the mask
;------------------------------
sz=delx*dely
outd=fltarr(delx,dely,8)
normmask=mask

;------------------------------
; Now the SCF values among the different spectra in the subcube
; are performed.  
;------------------------------
for a=0,delx-1 do begin ; these loops examine each neighboring spectra
for b=0,dely-1 do begin	; and choose minimum values

if a+xmin eq x and b+ymin eq y then goto,endloop1

rej=fitparams(a+xmin,b+ymin,5)	; are the neighbors rejected?
del=(min(scfmin(a,b,*),sub)); gets min del and the subscript of the entry
tau=((sub-q*vdisppt)*dv0)	; converts sub to velocity
s=stau(a,b,sub)			; gets corresponding s value
del_l=(min(noscale(a,b,*),sub2)) ;repeat with s=1 tau=float
tau_l=((sub2-q*vdisppt)*dv0)
del_s=(scfmin(a,b,q*vdisppt))	; repeat for s=float tau=0
s_s=stau(a,b,q*vdisppt)
del_0=(noscale(a,b,q*vdisppt)) ; s=1 tau=0
normmask(a,b,*)=normmask(a,b,*)*(rej eq 0)
outd(a,b,*)=[del,tau,s,del_l,tau_l,del_s,s_s,del_0]*normmask(a,b)*(rej eq 0)


endloop1:

endfor
endfor

;------------------------------
; If there are no viable spectra for a base spectrum to be correlated 
; with that spectrum is then rejected
;------------------------------
if total(normmask) eq 0 then begin
print,x,y,'<------ Cannot correlate with nearby spectra, rejecting'
goto,endloop3
endif

;------------------------------
; The average is completed along with the normalization.
;------------------------------
out=total(total(outd,1),1)/(total(normmask)-1)
; -1 because the mask includes the center spectra where the data do not
out(0)=1-sqrt(out(0))
out(3)=1-sqrt(out(3))
out(5)=1-sqrt(out(5))
out(7)=1-sqrt(out(7))	;final step in normalization

outdata(x,y,*)=[out,Tmax0,vdisp0,vpeak0,Tint0,sig0,Tmax0/sig0,$
abs(out(1)),abs(out(4))]

endloop3:

endfor 
endfor
;------------------------------
; Finally, some timekeeping routines, FYI
;------------------------------
fin=systime(1)	;stop the clock
print,'Total time=',fin-st	;fun math!
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