PRO scf,file,dv

;
; Example: 
;
; scf,'hcl2_c18o.fits',0.05
;
;Compute and plot SCF from a given spectral map file with velocity
;channel width = dv, in km/s. 

!x.thick=2.0
!y.thick=2.0
!p.thick=2.0
!p.charthick=2.0
!p.charsize=1.4

;
;Read data file:
;---------------
spec=readfits(file)
;
;It is assumed velocity is first coordinate, otherwise:
;
;spec=transpose(spec,[2,0,1]) ;if velocity is third coordinate.
;spec(where(spec ne spec))=0.0 ;if need to turn NaN into zeros.
s=size(spec)
print,s
nv=s(1)
N1=s(2)
N2=s(3)
v=findgen(nv)*dv-(dv*nv/2.)                 ;velocity channels
noise=sqrt(aver(spec(where(spec lt 0.))^2)) ;1-sigma noise
signal=sqrt(total(spec^2)/(N1*N2*nv))       ;rms signal
print,'noise:',noise
print,'rms(S):',signal
;
;Average spectrum:
;-----------------
temp=fltarr(nv)
for j=0,nv-1 do begin     
     scr=aver(spec(j,*,*))
     temp(j)=temp(j)+scr
endfor
window,0,xsize=400,ysize=300
plot,v,temp,xtitle='v [km/s]',ytitle='T [K]', $ 
     xr=[min(v)-0.2,max(v)+0.2],xs=1,ys=1, $
     yr=[min(temp)-0.1,max(temp)+0.1]     ;plot the average spectrum
;
;Maps of integrated and peak temperatures:
;-----------------------------------------
map_p=fltarr(N1,N2)
map_i=fltarr(N1,N2)
width=fltarr(N1,N2)
for i=0,N1-1 do begin
   for j=0,N2-1 do begin
      map_p(i,j)=max(spec(*,i,j))       ;peak temperature
      map_i(i,j)=dv*total(spec(*,i,j))  ;integrated temperature
      width(i,j)=map_i(i,j)/map_p(i,j)  ;line width
  endfor
endfor
index=where(map_p gt 3.*noise and map_i gt 0.0) ;3-sigma detection for computing <linewidth>
print,'Average Line-Width:',aver(width(index))/2.355 ;/2.355 to compare with sigma(v)
;
;Statistics of mean spectrum:
;----------------------------
tt=total(temp)
mean_v=total(temp*v)/tt
sig_v=sqrt(total(temp*(v-mean_v)^2)/tt)
print,'<v>:',mean_v
print,'sigma(v)',sig_v
print,'sigma(v)/line-width:',sig_v/(aver(width(index))/2.355) ;comparison sigma_v of mean spectrum/<w>
;
;Parameters for SCF:
;-------------------
mean_vl=fltarr(N1,N2) ;local average velocity
Q=fltarr(N1,N2)       ;spectrum quality
;
widthv=2. ; velocity window for computing the SCF, given as a
          ; multiple of sigma_v of the mean spectrum. 
;
;Put intensity to zero in velocity channels outside the velocity window:
;
for i=0,N1-1 do begin
    for j=0,N2-1 do begin
        tot=total(spec(*,i,j))
        mean_vl(i,j)=total(spec(*,i,j)*v)/tot          
        scr=min(abs((v-(mean_vl(i,j)-widthv*sig_v))),v1)
        scr=min(abs((v-(mean_vl(i,j)+widthv*sig_v))),v2)
        spec(0:v1,i,j)=0.0
        spec(v2:nv-1,i,j)=0.0
        Q(i,j)=1.*sqrt(aver(spec(v1:v2,i,j)^2))/noise  ;Spectrum Quality
    endfor
endfor
;
print,'<Q>=',aver(Q) ;average Spectrum Quality 
;
scr=where(Q eq 0.0,c)
print,'Fraction of pixels with Q=0:',c/(1.0*N1*N2)
scr=where(Q gt 0.0 and Q lt 2,c)
print,'Fraction of pixels with 0<Q<2:',c/(1.0*N1*N2)
;
;IMAGE of integrated map:
;------------------------
window,2,xsize=N1*4+1,ysize=N2*2
tvscl,rebin(map_i,N1*2,N2*2),0,0
;
;Choose spatial lags for the SCF:
;--------------------------------
lag=3.0*exp(findgen(35)/5.) 
scr=where(lag lt max([N1,N2])-1.,c)
lag=lag(0:c)
lag=[1,2,lag] ;impose the first two spatial lags equal to 1 and 2.
lags=c+1+2    ;number of lags
;
;The SCF could be computed with differences of spectra along x or y
;or both. Here we compare the spectrum at one position with all the 
;spectra along a circle around it with radius equal to "lag". 
;This is much more expensive, but gives better results.
;
;Generate grid to compute shifts:
;--------------------------------
nx=fix(N1)*2+1   ;use 2D complete circle of shifts 
ny=fix(N2)*2+1 
x=rebin(reform(indgen(nx)-nx/2,nx,1),nx,ny)
y=rebin(reform(indgen(ny)-ny/2,1,ny),nx,ny)
rad=sqrt(x^2+y^2) 
;
;Define variables:
;-----------------
test=fltarr(nx,ny)       ;for testing shifts 
avlag=fltarr(lags)
scf=fltarr(lags,N1,N2)   ;SCF images (one image per lag)
scf_av=fltarr(lags)      ;average of the SCF over the whole map 
                         ;(one value per lag)

for i=0,lags-2 do begin
    ;
    ;Compute shifts:
    ;---------------
    ind=where(rad ge lag(i) and rad lt min([lag(i+1),lag(i)+1.]),c)
    
    if (c gt 0) then begin
       avlag(i)=aver(rad(ind))
       ind_x=x(ind)  ;coordinates in the shift grid
       ind_y=y(ind)
       ;
       ;IMAGE for testing shift circle:
       ;-------------------------------
       ;test(ind)=i+2               
       ;imsz,nx,ny
       ;tvscl,rebin(test,nx,ny)
       ;
       ;Loop to shift over all values of "ind" (over the circle of radius = lag):
       ;
       buffer=fltarr(N1,N2)
       count=fltarr(N1,N2)

       for j=0,c-1 do begin
           ;
           scr=shift(spec,0,ind_x(j),ind_y(j)) ;shifted spectral map
           Q_shift=shift(Q,ind_x(j),ind_y(j))  ;shifted Q
           Q_scr=sqrt((Q^2+Q_shift^2)/2.)  ;take average of the two Q^2 values 
    
           ;
           ;Kill periodicity of 'shift' (for non periodic 
           ;boundary, such as observational data.):
           ;---------------------------------------------
           if (ind_x(j) gt 0) then scr(*,0:ind_x(j)-1,*)=0.0   
           if (ind_y(j) gt 0) then scr(*,*,0:ind_y(j)-1)=0.0   
           if (ind_x(j) lt 0) then scr(*,N1+ind_x(j):N1-1,*)=0.0   
           if (ind_y(j) lt 0) then scr(*,*,N2+ind_y(j):N2-1)=0.0 
           ;
           ;Compute SCF (individual shift) for each pixel on the map:
           ;---------------------------------------------------------
           num=(spec-scr)^2
           num=total(num,1)
           den=total(spec^2,1)+total(scr^2,1)
           scr=1.-sqrt(num/den)
           scr=scr/(1.-1./Q)    ;correction for spectrum quality Q
           ;
           ;exclude all pixels with Q<2 or negative SCF:
           index=where(Q gt 2.0 and Q_shift gt 2. and scr gt 0.0,c) 
           ;
           ;the above condition (scr gt 0) is necessary to avoid 
           ;counting the periodicity of the shift killed earlier.
           ;
           if (c gt 0) then begin
              count(index)=count(index)+1.  ;statistical weight of individual shift 
              buffer(index)=buffer(index)+scr(index) ;sum SCF of this individual shift
           endif
           ;

       endfor   ;j loop (shifts)
       ;
       index=where(count gt 0)
       buffer(index)=buffer(index)/(1.0*count(index))  ;from sum to average 
       scf(i,*,*)=buffer
       scf_av(i)=aver(buffer(index)) ;average over the map for pixels with Q>2.0 
       
       print,'DX, <SCF>:',avlag(i),scf_av(i)
       ;
       ;IMAGE of SCF:
       ;-------------
       index=where(Q lt 2.,c)
       if (c gt 0) then buffer(index)=0.0
       tvscl,rebin(buffer,N1*2,N2*2),N1*2,0
       ;
    endif
endfor          ;i loop (lags)

index=where(scf_av gt 0.0)
print,'DX:',avlag(index)
print,'<SCF>:',scf_av(index)

;Plot SCF versus lag:
;--------------------
scf_av=scf_av(index)
avlag=avlag(index)
window,3,xsize=500,ysize=300
plot_oo,avlag,scf_av,psym=6,ys=1,yr=[min(scf_av)-0.2,max(scf_av)+0.2], $
        xs=1,xtitle='r',ytitle='<S!D0!N(r)>',xr=[min(avlag)-0.2,max(avlag)+0.2]


STOP
END