PRO xcorr_fits, dfl, mode, gplot, xsave


;; calculates centroids of cross-correlation functions
;; from 2d array 

;; dfl is size of window to fit peak (in pixels / 2)
;;  mode = 0, fit Gaussian, 
;;       = 1, fit parabola
;;       = 2, compute moments
;; gplot = 0, no plots, 
;;       = 1, plot ccfs and fits 
;; xsave = 0, no save
;;       = 1, save last ccf and fit
 


COMMON fitsblock
COMMON xcorblock, ccfx, ccfy, ccf2y, ccf2n, nccf, afits


;; Announce programme and options

print, 'XCORR_FITS, dfl, model, gplot, xsave'
print
print, 'Analyse a series of CCFs'
print 
print, 'dfl   =',dfl, '; half-width (pixels) of region used for fit ' 
IF (mode eq 0) THEN $
    print, 'mode  =', mode, '; fit Gaussian'
IF (mode eq 1) THEN $
    print, 'mode  =', mode, '; fit parabola'
IF (mode eq 2) THEN $
    print, 'mode  =', mode, '; compute moments 0 - 3: 1 = centroid'
IF (gplot eq 0) THEN $
    print, 'gplot =', gplot, '; no graphics'
IF (gplot eq 1) THEN $
    print, 'gplot =', gplot, '; graphics on'
IF (xsave eq 0) THEN $
    print, 'xsave =', xsave, '; no saves'
IF (xsave eq 1) THEN $
    print, 'xsave =', xsave, '; save last xcorr and fit'
print


; Define window for fitting ccf
nfl=nccf/2-dfl
nfu=nccf/2+dfl
wfl=nccf/2-4*dfl
wfu=nccf/2+(4*dfl)

; Loop over spectra
FOR it = 0,nspec-1 DO BEGIN
   afit = [0.0,0.0,1.0,0.0,0.0,1.0E-8]
   afits(*,it) = afit

;  copy 1d ccf to work array and check valid
   ccfy = ccf2y(*,it) 
   IF ( stddev(ccfy) gt 1.e-6 ) THEN BEGIN

;    select central fraction of ccf
       v = ccfx(nfl:nfu)
       p = ccfy(nfl:nfu)


;    fit gaussian to peak
     IF (mode eq 0) THEN BEGIN
       gsb = p(0)                         ; get estimates for Gaussian
       gst = p(dfl)
       a0 = [gst-gsb,0.0,dfl*30,gsb,0.0,0.0]
       ccfit = gaussfit(v,p,afit,estimates=a0,nterms=6)

       IF (afit(1) GT dfl*30) THEN afit = [0.0,0.0,1.0,0.0,0.0,1.0E-6]
       IF (afit(0) LT 0.01)   THEN afit = [0.0,0.0,1.0,0.0,0.0,1.0E-6]
       IF TOTAL(FINITE(afit)) NE N_ELEMENTS(afit) THEN $
          afit = [0.0,0.0,1.0,0.0,0.0,1.0E-8]
       IF (ABS(afit(5)) LT 1.E-8)   THEN afit(5) = 1.0E-8
       afits(*,it) = afit
     ENDIF

;    fit parablola to peak
     IF (mode eq 1) THEN BEGIN
       prbola, v,p,a,as,t,ts,is,pbfit
       afits(1,it) = t(0)
       afits(0,it) = t(1)
     ENDIF

;    compute moments of line profile
     IF (mode eq 2) THEN BEGIN
       mom = fltarr(4)
       lp_moment, v,p,0.,mom
       afits(0:3,it) = mom
     ENDIF

;    plotting the individual ccfs - for debugging only
;     IF (gplot gt 0) THEN BEGIN
;        window,1
;        plot, v,p
;        oplot, v,p,psym=1
;        IF (mode eq 0) THEN oplot, v,ccfit, color=45
;        IF (mode eq 1) THEN oplot, v,pbfit, color=90
;     ENDIF


   ENDIF

ENDFOR

; IF (gplot gt 0) THEN window,0
; IF (gplot gt 0) THEN shade_surf, ccf2y(nfl:nfu,*)

IF (xsave eq 1) THEN BEGIN
   openw,lxfit,'xfit.data', /get_lun
   FOR ix=0,n_elements(v)-1 DO BEGIN
     IF (mode eq 0) THEN printf, lxfit, v(ix),p(ix),ccfit(ix), format='(3f10.3)'
     IF (mode eq 1) THEN printf, lxfit, v(ix),p(ix),pbfit(ix), format='(3f10.3)'
   ENDFOR
   free_lun, lxfit
ENDIF

END