;+
; FUNCTION: THM_QFIT,  data, phs, phsf=phsf, es=es, ec=ec, zero=zero,
;                  do_sigma=do_sigma, sigma=sigma,
;                  period=period, slide=slide, n_fitpts=n_fitpts, 
;                  out=out, max_err=max_err, bad_pts=bad_pts                 
;       
; PURPOSE: User unfriendly fit routine. 
;          NOTE: REQUIRES CONTINUOUS BUFFERS!
;          NOT FOR GENERAL USE. 
;
; INPUT: 
;       data -        REQUIRED. A DATA ARRAY -  NOT A STRUCTURE!
;                     RECOMMEND: dat_in does not have NANS. One nan
;                                may destroy entire period.
;       phs -         REQUIRED. An array of phases.
;                     IMPORTANT! All phases must be valid! No NANs!
;
; KEYWORDS: 
;       period   -    INPUT. Fit period. DEFAULT = 4pi
;       slide    -    NO LONGER AN OPTION!
;                     NOTE: Slide will be forced to be pi/2.
;       n_fitpts -    INPUT. Number of points per fit. DEFAULT = 64.
;       do_sigma -    OPTION. /do_sigma fills sigma.
;       out      -    OPTION. /out fills bad_pts.
;                     NOTE: FOR OUTLYER REJECTION, SEE BELOW!
;       max_err  -    INPUT. Maximum allowable error of bad_pts. DEFAULT=25 nT
;
; OUTPUT: 
;       phsf    -     OUTPUT. Phase of fit -> Time of fit.
;       es      -     OUTPUT. Sin phase of fit.
;       ec      -     OUTPUT. Cos phase of fit.
;       zero    -     OUTPUT. Zero level of fit.
;       sigma   -     OUTPUT. Deviation in nT. Only filled if /do_sigma
;       bad_pts -     OUTPUT. Deviation in nT. Only filled if /out
;
; CALLING: 
;       SEE ff_magdc for an example. MUST USE fa_fields_bufs, and 
;       ff_zero_crossing first.
;
; OUT-LYING POINTS REJECTION: 
;       The program must be iterated. For example, below shows one interation:
;
;       ff_qfit,data,phs,phsf=phsf,es=es,ec=ec,zero=zero,/out, bad_pts=bad_pts
;       data(bad_pts) = !values.f_nan
;       index = where(finite(data), n_index)
;       IF n_index GT 0 then BEGIN
;           data = data(index)
;           phs  = phs(index)
;           ff_qfit,data,phs,phsf=phsf,es=es,ec=ec,zero=zero
;       ENDIF ELSE print, "Big trouble! No valid points."
;
; HISTORY: 
;   INITIAL VERSION: REE 97-10-05
;   Space Sciences Lab, UCBerkeley
;   2011-05-20: Added the keyword FAIL. 
;         Jianbao Tao (JBT), CU/LASP.
; 
;
; 
;-
pro thm_qfit, data, phs, phsf=phsf, es=es, ec=ec, zero=zero, $
             do_sigma=do_sigma, sigma=sigma, $
             period=period, slide=slide, n_fitpts=n_fitpts, $
             out=out, max_err=max_err, bad_pts=bad_pts, $
             fail = fail

fail = 0

two_pi    = 2.d*!dpi
pi_over_2 = !dpi/2.d
ec        = !values.d_nan  ; ERROR MARKER

; CHECK KEYWORDS
if not keyword_set(period) then period = 2.d*two_pi $
ELSE BEGIN 
    period = double(long(period/two_pi + 0.0001)) * two_pi
    if period LT two_pi then period = two_pi
ENDELSE
n_per =  long(period/two_pi + 0.0001)
if not keyword_set(n_fitpts) then n_fitpts = 64 
n_fitpts = long( n_fitpts/(n_per*4) ) * n_per*4
if (n_fitpts LT n_per*4) then n_fitpts = n_per*4
n_slide = n_fitpts/(n_per*4)


; START BY CHECKING INPUTS.
; CHECK PHS.
npts = n_elements(phs)
IF npts LT 4 then BEGIN
    dprint, "STOPPED!"
    dprint, "Need to give phase array."
    fail = 1
    return
ENDIF

IF size(/type,phs) NE 5 then BEGIN
    dprint, 'STOPPED! '
    dprint, 'Phs must be double array.'
    fail = 1
    return
ENDIF

IF phs(0) LT 0 THEN BEGIN
    n_skip = long( abs(phs(0)) / two_pi ) + 1l
    phs  = phs + n_skip*two_pi
    phs0 = n_skip*two_pi
ENDIF ELSE phs0 = 0

; CHECK DAT.
dat_type = size(/type,data)
IF dat_type gt 5 or dat_type lt 1 then BEGIN
    dprint, 'STOPPED! '
    dprint, 'Dat byte, int, long, float, or double array.'
    fail = 1
    return
ENDIF
IF n_elements(data) NE npts then BEGIN
    dprint, "STOPPED!"
    dprint, "Data and phase must have the same number of points."
    fail = 1
    return
ENDIF

; CREATE TEMPORARY FIT ARRAY.
max_nfit = long( (phs(npts-1) - phs(0))/pi_over_2 + 9 )
qfit = dblarr(max_nfit)*!values.d_nan
base = dblarr(max_nfit)*!values.d_nan
phsf = dblarr(max_nfit)

; SET UP BASIS FUNCTION MATRIX A
;fix_phs = dindgen(n_fitpts)*n_per*two_pi/n_fitpts $
;    - n_per*two_pi/2.d + n_per*!dpi/n_fitpts
d_phs      = period/(n_fitpts*2)                  ; d_phs needed to detrend!
fix_phs_sm = dindgen(n_fitpts)*period/n_fitpts + d_phs
fix_cos    = cos(fix_phs_sm)


; MAKE UP FIXED PHASE ARRAY
phs_zero    = phs(0) - ( phs(0) mod pi_over_2 ) + pi_over_2 + d_phs
phs_end     = phs(npts-1) - ( phs(npts-1) mod pi_over_2 ) + d_phs
n_phases    = long( (phs_end - phs_zero) * n_fitpts / period + 0.0001)

IF n_phases le 0 then BEGIN
    dprint, "STOPPED!"
    dprint, "Negative phase values"
    fail = 1
    return
ENDIF

fix_phs     = dindgen(n_phases)*period/n_fitpts + phs_zero

; DETERMINE WHICH PHASE QUADRANT WE ARE IN
phs_quad    = ( phs(0) / pi_over_2 ) + 1.001
n_quad      = (long(phs_quad) mod 4)

; FORCE THE DATA TO MATCH THE FIXED PHASE ARRAY.
fix_dat  = interpol(data,phs,fix_phs, /spline)
fix_strt = 0l
fix_end  = long(n_fitpts-1)
nfits       = 0l

; LOOP THROUGH THE BUFFER
WHILE fix_end LT n_phases DO BEGIN

    ; DO THE QUICK FIT!
    qfit(nfits) = total( fix_dat(fix_strt:fix_end) * fix_cos ) * 2.d / n_fitpts 
    base(nfits) = total( fix_dat(fix_strt:fix_end) ) / n_fitpts
    phsf(nfits) = fix_phs(fix_strt) + period/2.0
    nfits       = nfits + 1
    fix_strt    = fix_strt + n_slide
    fix_end     = fix_end  + n_slide

ENDWHILE

; STOP HERE IF THERE ARE NOT ENOUGH FITS
if nfits LT 5 then begin
   dprint, 'nfits = ', nfits
   dprint, 'Not enough fits. Exiting...'
    fail = 1
   return
endif

; EXTRACT THE DATA - THIS IS A TRUE BRAIN TWISTER - TALK TO REE TO MODIFY
kerz  = [0.125d, 0.25d, 0.25d, 0.25d, 0.125d]
kers  = [0.5d, 0.5d]

; EXTRACT COS FITS
di     = n_quad mod 2				; di = 1 is odd case.
n_ec   = (nfits+1-di)/2
index  = lindgen(n_ec)             		; INDEX IS UNIFORM
sgn    = (index mod 2) * 2 - 1 			; sgn is actually -sgn
if (n_quad EQ 0 OR n_quad EQ 3) then sgn = -sgn ; Flip if first COS is 0
ec     = double( qfit(index*2 + di) ) * sgn	; COS is even qfits
n_ec   = (nfits+1-di)/2

; EXTRACT OFFSETS - Keeps edges correct.
base   = [base(2), base(3), base(0:nfits-1), base(nfits-4), base(nfits-3)] 
base   = convol(base,kerz, /edge_t) 
zero   = base(index*2 + di + 2)			; See edge correction.

; EXTRACT PHASES
phsf   = phsf(index*2 + di) - d_phs  ; STILL NEED TO SUBTRACT phs0!

; EXTRACT SIN FITS
n_es   = (nfits+di)/2
index  = lindgen(n_es) 		               	; INDEX IS UNIFORM
sgn    = (index mod 2) * 2 - 1 			; sgn is actually -sgn
if (n_quad EQ 0 OR n_quad EQ 1) then sgn = -sgn ; Flip if first SIN is 1
es     = double( qfit(index*2 + 1 - di) ) * sgn  ; SIN is the odd qfits

; BELOW IS A BLOODY MESS TO AVERAGE SIN TO PROPER ARRAY SIZE
if (n_es LT n_ec) then es = [es,es(n_es-1)]
if (di EQ 1) AND (n_es EQ n_ec) then es = [es,es(n_es-1)]
es     = convol(es,kers, /edge_t)		; AVERAGE PAIRS
if (n_es GT n_ec) then es = es(1:n_es-1)
if (di EQ 1) AND (n_es EQ n_ec) then es = es(1:n_es)

; WE ARE DONE EXCEPT FOR OUT-LYING ERRORS AND SIGMA
if NOT keyword_set(max_error) then max_error=25.d ; nT 

; CHECK FOR do_sigma or out
IF keyword_set(do_sigma) or keyword_set(out) then BEGIN

    ; SET UP ARRAYS FOR SIGMA AND OUT
    z_big   = interpol(zero,phsf,fix_phs)
    c_big   = interpol(ec,phsf,fix_phs)
    s_big   = interpol(es,phsf,fix_phs)
    fit_big = z_big + c_big * cos(fix_phs) + s_big * sin(fix_phs)

    strt   = phsf(0) - period/2.0 + d_phs
    index  = where(abs(fix_phs-strt) LT d_phs, n_index)

    IF ( n_index EQ 1 AND keyword_set(do_sigma) ) then begin
        sigma  = dblarr(n_ec)
        strt   = index(0)
        stop   = strt + n_fitpts - 1
        FOR i    = 0L, n_ec-1 DO BEGIN
            dif      = fit_big(strt:stop)-fix_dat(strt:stop)
            sigma(i) = sqrt( total(dif*dif) / (n_fitpts-3) )
            strt     = strt + n_slide*2
            stop     = stop + n_slide*2
        ENDFOR
    ENDIF

    IF keyword_set(out) then BEGIN

        ; EXPAND FIX_PHS TO COVER ALL DATA
        append_low  = fix_phs(0:2*n_fitpts-1)                 - 2.d*period
        append_high = fix_phs(n_phases-2*n_fitpts:n_phases-1) + 2.d*period
        fix_phs     = [append_low,fix_phs,append_high]

        ; SET UP ARRAYS FOR SIGMA AND OUT
        z_big   = interpol(zero,phsf,fix_phs)
        c_big   = interpol(ec,phsf,fix_phs)
        s_big   = interpol(es,phsf,fix_phs)
        fit_big = z_big + c_big * cos(fix_phs) + s_big * sin(fix_phs)

        fit_big = interpol(fit_big,fix_phs, phs, /spline)
        bad_pts = where(abs(data-fit_big) GT max_error)

    ENDIF
ENDIF ; End of out and do_sigma options.

; FIX PHASES
phs  = phs  - phs0
phsf = phsf - phs0
return

END

