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.
;------------------------------
error_key=fix(fitparams(x,y,5))
if error_key ne 0 then begin
print,'Spectrum at',fix(x),fix(y),' rejected.'
if (error_key/16 ne 0) then print,'Rejected because of spectrum could not be smoothed to the desired signal-to-noise.'
error_key=error_key mod 16
if (error_key/8 ne 0) then print,'Rejected because signal-to-noise could not be calculated'
error_key=error_key mod 8
if (error_key/4 ne 0) then print,'Rejected because signal-to-noise is too low.'
error_key=error_key mod 4
if (error_key/2 ne 0) then print,'Rejected because of insufficient baseline on one side of the peak.'
error_key=error_key mod 2
if (error_key/1 ne 0) then print,'Rejected because no spectrum is in array at this point.'


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


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

;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

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)
if total(where(outd ne outd)) ne -1 then outd(where(outd ne outd))=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
outdata(*,*,8:12)=fitparams(*,*,0:4)
outdata(*,*,13)=fitparams(*,*,8)
;------------------------------
; 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