;+
; PRO: THM_EFI_CLEAN_EFW
;
; PURPOSE:
;    FIXES HF, REMOVES SPIKES, AND LOW FREQUENCY.
;
;   This routine is designed for general use. It will take EFW data in
;   DSL as input, clean up the data, and transform the data into FAC, and return
;   the cleaned EFP data in DSL/FAC.
;
;   The routine can work on one-day long data, but if only a single particle
;   burst data is interesting, use the keyword TRANGE to specify the time range
;   of that burst can greatly reduce run time.
;
;   NOTE: 1) Spike removal does not always work, so don't be surprised if there
;         still are spikes left on the data. Hopefully the spike issue can be
;         resolved in the future.
;         2) Make sure thx_state_spinper is available! Best to set up for gsm.
;         Make sure timespan is set!!!!
;
; EXAMPLES:
;     For typical use, see thm_crib_cleanefw.pro and thm_crib_cleanefi.pro which
;     are typically located under TDAS_DIR/idl/themis/examples.
;
; INPUT:
;
; KEYWORDS (general):
;    probe          NEEDED: Program does only one sc at a time!
;    Ename          OPTIONAL: Valid TPLOT name for E WAVE, DEFAULT = 'thX_efw'
;                   (will fetch)
;    Bdslname       OPTIONAL, Valid TPLOT name for B, DEFAULT = 'thX_fgh_dsl'
;                   (will fetch)
;    trange         OPTIONAL. Time range.
;                   USE TRANGE TO ISOLATE SINGLE PARTICLE BURST AND GREATLY
;                   REDUCE RUN TIME!!!!
;    talk           OPTIONAL. Plots diagnostics. Can be fun. Slow. DEFAULT = 0
;    EfpName        OPTIONAL: Valid TPLOT name for EFP, DEFAULT = 'thX_efp'
;                   (will fetch)
;                   EFP is used to find spike positions.
;    status         OPTIONAL: A named variable to return the exiting status. If
;                    exit successfully, status = 0, otherwise status = 1
;    EFWHEDACNAME   OPTIONAL: Valid TPLOT name for the AC-coupled E-field header info
;                   DEFAULT = 'thX_efw_hed_ac'
;
; KEYWORDS (Filter):
;    FPole          (OBSOLETE) Kept for backward compatibility.
;
; KEYWORDS (for output tplot variables):
;    Edslname   OPTIONAL: Name for the cleaned efp in DSL.
;               DEFAULT = Ename + '_clean_dsl'
;    Efacname   OPTIONAL: Name for the cleaned efp in DSL.
;               DEFAULT = Ename + '_clean_fac'
;
; KEYWORDS (Remove_Spikes):
;    SpikeRemove    If entered as 0, suppresses spike removal.
;                   DEFAULT = 1 (Spikes are removed.)
;    SpikeNwin      Number of points in spike search window. DFLT = 16
;    SpikeSig       Sigma of spikes. DFLT = 5
;                   (Sometimes Sig = 6 works better)
;    SpikeSigmin    Minimun of sigma. DFLT = 0.01 mV/m.
;                   (Sometimes Sigmin = 0 works better)
;    SpikeNfit      Number of points in the fit window. DFLT = 16
;    SpikeFit       If set, will do a Gaussian fit to Spikes. DFLT = 0
;
; HISTORY:
;   2009-05-06: REE. First Release.
;   2009-06-04: Jianbao Tao (JBT), CU/LASP.
;               The metadata (dlim.data_att) is modified properly.
;   2009-06-17: JBT, CU/LASP.
;               Keyword STATUS added.
;   2010-04-08: JBT, CU/LASP.
;               The high-pass filter is changed from an FFT filter to an
;               FIR filter. Some modification is made to be compatible with
;               two different sample rates 8192 Hz and 16384 Hz.
;   2010-09-13: JBT, CU/LASP.
;               Cleaned up documents.
;               Added keyword dt in calling thm_efi_fix_freq_and_phase.
;   2012-07-26: JBT, SSL, UC Berkeley.
;               1. Incorporated C. Cully's calibration for
;                     1) Plasma-probe coupling
;                     2) Anti-aliasing Bessel filter
;                     3) ADC response
;                     4) DFB response
;   2012-08-30: CMC, University of Calgary.
;               Fixed problems with AC-coupled data sampled at < 16 ksps
;
; VERSION:
;
; $LastChangedBy $
; $LastChangedDate $
; $LastChangedRevision $
; $URL $
;
;-

pro thm_efi_clean_efw, probe=probe, Ename=Ename, Bdslname=Bdslname, $
        trange=trange, talk=talk, EfpName=EfpName, $       ; GENERAL
        FPole=FPole, $                                     ; FILTER
        EFWHEDACNAME=efwhedacname, $                    ; optional name for AC header info
        Edslname = Edslname, $                          ; OUTPUT NAME
        Efacname = Efacname, $                          ; OUTPUT NAME
        SpikeRemove=SpikeRemove, SpikeNwin=SpikeNwin, $    ; REMOVE_SPIKES
        SpikeSig=SpikeSig, SpikeSigmin=SpikeSigmin, $      ; REMOVE_SPIKES
        SpikeNfit=SpikeNfit, SpikeFit=SpikeFit, $          ; REMOVE_SPIKES
        diagnose=diagnose, wt=wt, $                        ; Programmer's use
        status = status

compile_opt idl2

status = 0

; # CHECK INPUTS - GET NEEDED DATA #
IF not keyword_set(probe) then BEGIN
  print, 'THM_EFI_CLEAN_EFW: SC not set. Exiting...'
  status = 1
  return
ENDIF
sc = probe[0]

; LOAD STATE
thm_load_state, probe = probe, /get_support, trange=trange

; SETDEFAULTS
if n_elements(FPole) EQ 0       then FPole       = 5.0D   ; FILTER
if n_elements(SpikeRemove) EQ 0 then SpikeRemove = 1      ; REMOVE_SPIKES
;if n_elements(SpikeNfit) EQ 0   then SpikeNfit   = 400L   ; REMOVE_SPIKES

; CHECK FOR EFW DATA
if not keyword_set(Ename) then Ename = 'th' + sc + '_efw'
IF spd_check_tvar(Ename) then BEGIN
  get_data, ename[0], data=E, dlim=elim
  if ~strcmp(elim.data_att.coord_sys, 'dsl', /fold) then begin
    print, 'THM_EFI_CLEAN_EFW: ' + ename[0] + ' is not in DSL. Exiting...'
    return
  endif
ENDIF ELSE BEGIN
  thm_load_efi, probe=sc, datatype=['efw', 'vaw'], coord='dsl', trange=trange, /get_support_data
  Ename = 'th' + sc + '_efw'
  get_data, ename[0], data=E, dlim=elim
  IF size(/type,E) NE 8 then BEGIN
    print, 'THM_EFI_CLEAN_EFW: Cannot get electric field data. Exiting...'
    status = 1
   return
  ENDIF
ENDELSE




; CHECK FOR MAG DATA
if not keyword_set(Bdslname) then Bdslname = 'th' + sc + '_fgh_dsl'
IF spd_check_tvar(Bdslname) then BEGIN
  get_data, Bdslname[0], data=Bdsl, dlim=blim
ENDIF ELSE BEGIN
  print, 'THM_EFI_CLEAN_EFP: Mag data not stored in dsl. Fetching...'
  thm_load_fgm, probe=sc, datatype = ['fgh'], coord=['dsl'], trange=trange, $
      level = 2
  Bdslname = 'th' + sc + '_fgh_dsl'
  get_data, Bdslname[0], data=Bdsl, dlim=blim
  IF size(/type,Bdsl) NE 8 then BEGIN
    print, 'THM_EFI_CLEAN_EFW: Cannot get MAG data. Exiting...'
    status = 1
    return
  ENDIF
ENDELSE
Bdsl = {x:Bdsl.x, y:Bdsl.y}

; CHECK FOR SPIN DATA
SpinName = 'th' + sc + '_state_spinper'
IF spd_check_tvar(Spinname) then BEGIN
 get_data, SpinName, data=SpinPer
ENDIF ELSE BEGIN
  thm_load_state, probe=sc, datatype='spinper'
  get_data, SpinName, data=SpinPer
  IF size(/type,SpinPer) NE 8 THEN BEGIN
    print, 'THM_EFI_CLEAN_EFW: Cannot get spin period. Exiting...'
    status = 1
    return
  ENDIF
ENDELSE

; GET EFP DATA FOR SPIKE FINDER
IF keyword_set(SpikeRemove) then BEGIN
  if not keyword_set(EfpName) then EfpName = 'th' + sc + '_efp'
  IF spd_check_tvar(EfpName) then BEGIN
    get_data, EfpName[0], data=Efp, dlim = tmpdlim
    if ~strcmp(tmpdlim.data_att.coord_sys, 'dsl', /fold) then begin
      print, 'THM_EFI_CLEAN_EFW: ' + Efpname[0] + ' is not in DSL. Exiting...'
      return
    endif
  ENDIF ELSE BEGIN
    thm_load_efi, probe=sc, datatype=['efp', 'vap'], coord='dsl', trange=trange, /get_support_data
    get_data, EfpName[0], data=Efp
    IF size(/type,Efp) NE 8 then BEGIN
      print, 'THM_EFI_CLEAN_EFW: Cannot get EFP data. Spikes cannot be removed.'
      SpikeRemove = 0
    ENDIF
  ENDELSE
ENDIF

; CLIP DATA TO RANGE
IF keyword_set(trange) then BEGIN
   trange_clip, E, trange[0], trange[1], /data, BadClip=BadEclip
   trange_clip, Bdsl, trange[0], trange[1], /data, BadClip=BadBclip
   trange_clip, SpinPer, trange[0]-60.d, trange[1]+60.d, /data, BadClip=BadSclip
   IF (keyword_set(BadEclip) OR keyword_set(BadBclip) OR keyword_set(BadSclip) ) THEN BEGIN
     print, 'THM_EFI_CLEAN_EFW: Problem with trange clip. Exiting...'
     print, '0=OK; 1=Problem. E:', BadEclip, 'B:', BadBclip, 'Spin:', BadSclip
     status = 1
     return
   ENDIF
ENDIF

; MAKE ARRAY FOR FAC STORAGE
E = {x:E.x, y:E.y}
Efac = E

; ## IDENTIFY INDIVIDUAL WAVE BURSTS  ##
tE   = E.x
thm_lsp_find_burst, E, istart=bstart, iend=bend, nbursts=nbursts, mdt=mdt
srate = 1. / mdt
if n_elements(SpikeNfit) EQ 0   then begin
   if abs(srate-8192.) lt 10 then SpikeNfit   = 400L
   if abs(srate-16384.) lt 10 then SpikeNfit   = 800L
endif

; CHECK FOR EFW HEADER TPLOT HANDLE KEYWORD SETTING
test  = (size(efwhedacname,/type) ne 7) or (not keyword_set(efwhedacname))
if (test[0]) then efwhacname = 'th'+sc[0]+'_efw_hed_ac' else efwhacname = efwhedacname[0]
; Flag for AC/DC coupled (needed for deconvolving instrument response)
get_data,efwhacname[0],data=efw_ac
if (size(efw_ac,/type) ne 8) then begin
  message,'No efw header info TPLOT handle found --> exiting prematurely',/continue,/informational
  status = 1
  return
endif
;get_data,'th'+sc+'_efw_hed_ac',data=efw_ac
if total(efw_ac.y ne efw_ac.y[0]) ne 0 then begin
  print,'THM_EFI_CLEAN_EFW: Error: EFW data switches coupling (AC/DC) during the requested interval.'
  print,'                   Please split this interval at ',time_string(efw_ac.x[min(where(efw_ac.y ne efw_ac.y[0]))])
  print,'                   and reprocess as separate intervals. Exiting...'
  status = 1
  return
endif
efw_ac=efw_ac.y[0]

; Setup kernel for deconvolving EFI response that includes
;   plasma-probe response
;   anti-aliasing filter response
;   ADC interleaving timing
;   DFB digital filter response
fsample = srate
kernel_length = 1024L
df = fsample / double(kernel_length)
f = dindgen(kernel_length)*df
f[kernel_length/2+1:*] -= double(kernel_length) * df
thm_comp_efi_response, 'SPB', f, SPB_resp, rsheath=5d6, /complex_response
thm_comp_efi_response, 'AXB', f, AXB_resp, rsheath=5d6, /complex_response
if efw_ac then begin
  E12_resp =thm_adc_resp('E12AC',f)
;   E12_resp = 1 / (SPB_resp * thm_eac_filter_resp(f) $
;                            * thm_adc_resp('E12AC',f) $
;                            * thm_dfb_dig_filter_resp(f, fsample,/EAC))
;   E34_resp = 1 / (SPB_resp * thm_eac_filter_resp(f) $
;                            * thm_adc_resp('E34AC',f) $
;                            * thm_dfb_dig_filter_resp(f, fsample,/EAC))
;   E56_resp = 1 / (AXB_resp * thm_eac_filter_resp(f) $
;                            * thm_adc_resp('E56AC',f) $
;                            * thm_dfb_dig_filter_resp(f, fsample,/EAC))
  E12_resp = 1 / (SPB_resp * thm_eac_filter_resp(f) $
                           * thm_adc_resp('E12AC',f) $
                           * thm_dfb_dig_filter_resp(f, fsample))
  E34_resp = 1 / (SPB_resp * thm_eac_filter_resp(f) $
                           * thm_adc_resp('E34AC',f) $
                           * thm_dfb_dig_filter_resp(f, fsample))
  E56_resp = 1 / (AXB_resp * thm_eac_filter_resp(f) $
                           * thm_adc_resp('E56AC',f) $
                           * thm_dfb_dig_filter_resp(f, fsample))
  E12_resp[0] = 0
  E34_resp[0] = 0
  E56_resp[0] = 0
endif else begin
  E12_resp = 1 / (SPB_resp * bessel_filter_resp(f,4096,4) $
    * thm_adc_resp('E12DC',f) * thm_dfb_dig_filter_resp(f, fsample))
  E34_resp = 1 / (SPB_resp * bessel_filter_resp(f,4096,4) $
    * thm_adc_resp('E34DC',f) * thm_dfb_dig_filter_resp(f, fsample))
  E56_resp = 1 / (AXB_resp * bessel_filter_resp(f,4096,4) $
    * thm_adc_resp('E56DC',f) * thm_dfb_dig_filter_resp(f, fsample))
endelse

; Transfer kernel into time domain: take inverse FFT and center
E12_resp = shift((fft(E12_resp,1)), kernel_length/2) / kernel_length
E34_resp = shift((fft(E34_resp,1)), kernel_length/2) / kernel_length
E56_resp = shift((fft(E56_resp,1)), kernel_length/2) / kernel_length

; START LOOP OVER INDIVIDUAL BURSTS
FOR ib=0L, nbursts-1 DO BEGIN

  ; BREAK OUT DATA
  t  = tE[bstart[ib]:bend[ib]]
  Ex = E.y[bstart[ib]:bend[ib],0]
  Ey = E.y[bstart[ib]:bend[ib],1]
  Ez = E.y[bstart[ib]:bend[ib],2]
  Ttemp = [min(t)-mdt/2, max(t)+mdt/2]
  nt_burst = n_elements(t)

  print, ''
  print, 'BURST: ', string(ib+1,form='(I0)'), ' out of: ', $
    string(nbursts,form='(I0)'), '; sample rate: ', $
    string(srate, form='(I0)'), '; burst data points: ', $
    string(nt_burst, form='(I0)')
  print, ''

  if nt_burst lt kernel_length then begin
     Efac.y[bstart[ib]:bend[ib],0]    = !values.d_nan
     Efac.y[bstart[ib]:bend[ib],1]    = !values.d_nan
     Efac.y[bstart[ib]:bend[ib],2]    = !values.d_nan
     print, 'Burst too short. Skipping...'
     print, ''
     continue
   endif

  ; Check NaNs
  indx = where(finite(Ex), nindx)
  indy = where(finite(Ey), nindy)
  indz = where(finite(Ez), nindz)
  if nindx le nt_burst/2. or nindy le nt_burst/2. or nindz le nt_burst/2. $
    then begin
    Efac.y[bstart[ib]:bend[ib],0]    = !values.d_nan
    Efac.y[bstart[ib]:bend[ib],1]    = !values.d_nan
    Efac.y[bstart[ib]:bend[ib],2]    = !values.d_nan
    print, 'Too many NaNs in the burst. Skipping...'
    continue
  endif
  if nindx le nt_burst then Ex = interpol(Ex[indx], t[indx], t)
  if nindy le nt_burst then Ey = interpol(Ey[indy], t[indy], t)
  if nindz le nt_burst then Ez = interpol(Ez[indz], t[indz], t)

  ; Minimize offset   - by JBT
  Ex = Ex - median(Ex)
  Ey = Ey - median(Ey)
  Ez = Ez - median(Ez)


  ; CALCULATE SPIN PERIOD
  ind = where( (SpinPer.x GE (Ttemp[0]-60.d)) AND $
    (SpinPer.x LE (Ttemp[1]+60.d)), nind)
  IF nind EQ 0 then BEGIN
    print, 'THM_EFI_CLEAN_EFW: Spin period missing during burst. Exiting...'
    status = 1
    return
  ENDIF
  per = median(spinper.y[ind])

  ; DO HIGH-PASS ON E
  if keyword_set(talk) then print, 'HIGH-PASS FILTER'
  Exf = thm_lsp_filter_highpass(Ex, mdt, freqlow=flow)
  Eyf = thm_lsp_filter_highpass(Ey, mdt, freqlow=flow)
  Ezf = thm_lsp_filter_highpass(Ez, mdt, freqlow=flow)

  ; REMOVE SPIKES
  IF keyword_set(SpikeRemove) then BEGIN
    thm_lsp_remove_spikes, t, Exf, Eyf, Ezf, per, Efp=Efp, Nwin=SpikeNwin, $
        SpikeSig=SpikeSig, Sigmin=SpikeSigmin, Nfit=SpikeNfit, Fit=SpikeFit, $
        talk=talk, diagnose=diagnose, wt=wt
  ENDIF

;   ; FIX FREQUENCY AND PAHSE
;   Exf = thm_efi_fix_freq_and_phase(Exf, dt = mdt)
;   Eyf = thm_efi_fix_freq_and_phase(Eyf, dt = mdt)
;   Ezf = thm_efi_fix_freq_and_phase(Ezf, dt = mdt,/ax)

  ; De-convolve instrument responses.
  b_length = 8 * kernel_length
  while b_length gt nt_burst do b_length /= 2
  print, 'b_length = ', b_length

  ; Remove NaNs
  indx = where(finite(Exf), nindx)
  indy = where(finite(Eyf), nindy)
  indz = where(finite(Ezf), nindz)
  if nindx ne nt_burst then Exf = interpol(Exf[indx], t[indx], t)
  if nindy ne nt_burst then Eyf = interpol(Eyf[indy], t[indy], t)
  if nindz ne nt_burst then Ezf = interpol(Ezf[indz], t[indz], t)

  ;-- Zero-pad data to account for edge wrap
  Exf = [Exf, fltarr(kernel_length/2)]
  Eyf = [Eyf, fltarr(kernel_length/2)]
  Ezf = [Ezf, fltarr(kernel_length/2)]

  ;-- Deconvolve transfer function
  Exf = shift(blk_con(E12_resp, Exf, b_length=b_length),-kernel_length/2)
  Eyf = shift(blk_con(E34_resp, Eyf, b_length=b_length),-kernel_length/2)
  Ezf = shift(blk_con(E56_resp, Ezf, b_length=b_length),-kernel_length/2)

  ;-- Remove the padding
  Exf = Exf[0:nt_burst-1]
  Eyf = Eyf[0:nt_burst-1]
  Ezf = Ezf[0:nt_burst-1]

  ; SAVE E DATA
  E.y[bstart[ib]:bend[ib],0]    = Exf
  E.y[bstart[ib]:bend[ib],1]    = Eyf
  E.y[bstart[ib]:bend[ib],2]    = Ezf

  ; Transform data into FAC.
  Eclip = E
  trange_clip, Eclip, Ttemp[0], Ttemp[1], /data
  store_data, 'Eclip', data=Eclip, dlim=elim

  ; CLIP MAGNETOMETER DATA
  Bclip = Bdsl
  trange_clip, Bclip, Ttemp[0]-3.d, Ttemp[1]+3.d, /data, badclip=badclip
  if keyword_set(badclip) then begin
     Efac.y[bstart[ib]:bend[ib],0]    = !values.d_nan
     Efac.y[bstart[ib]:bend[ib],1]    = !values.d_nan
     Efac.y[bstart[ib]:bend[ib],2]    = !values.d_nan
     continue
  endif

  store_data, 'Bclip', data=Bclip, dlim=blim
  nsmpts = 11  ; smooth points
  if n_elements(Bclip.x) le nsmpts then nsmpts = n_elements(data.x)/2
  nsmpts >= 1
  tsmooth2, 'Bclip', nsmpts, newname='Bclip'

  ; GO TO FAC COORDINATES (JIANBAO)
  thm_lsp_clean_timestamp, 'Bclip'
  thm_lsp_clean_timestamp, 'Eclip'
  thm_fac_matrix_make, 'Bclip', other_dim='zdsl', $
             newname='th'+sc+'_fgh_fac_mat'
  tvector_rotate, 'th'+sc+'_fgh_fac_mat', 'Eclip', $
          newname='Eclip', error=error

  ; GET ECLIP AND SAVE
  get_data, 'Eclip', data = Eclip
  Efac.y[bstart[ib]:bend[ib],0]    = Eclip.y[*,0]
  Efac.y[bstart[ib]:bend[ib],1]    = Eclip.y[*,1]
  Efac.y[bstart[ib]:bend[ib],2]    = Eclip.y[*,2]
ENDFOR
; ## END OF LOOP

; STORE E DATA
; add BAND to data_att -JBT
flow = floor((flow + 2.5)/5.) * 5.
if abs(srate-8192.) lt 10 then bandmsg = '~' + string(flow, format='(I0)') + $
         ' Hz -- ~3.3 kHz'
if abs(srate-16384.) lt 10 then bandmsg = '~' + string(flow, format='(I0)') + $
         ' Hz -- ~6.6 kHz'
data_att = {DATA_TYPE: elim.data_att.DATA_TYPE, $
            COORD_SYS: elim.data_att.COORD_SYS, $
            UNITS: elim.data_att.UNITS, $
            CAL_PAR_TIME: elim.data_att.CAL_PAR_TIME, $
            OFFSET: elim.data_att.OFFSET, $
            EDC_GAIN: elim.data_att.EDC_GAIN, $
            EAC_GAIN: elim.data_att.EAC_GAIN, $
            BOOM_LENGTH: elim.data_att.BOOM_LENGTH, $
            BOOM_SHORTING_FACTOR: elim.data_att.BOOM_SHORTING_FACTOR, $
            DSC_OFFSET: elim.data_att.DSC_OFFSET, $
            BAND: bandmsg}
                     ; BAND - the freq band of the data

if ~keyword_set(Edslname) then Edslname = Ename + '_clean_dsl'
ename2 = Edslname
dlim = {CDF: elim.cdf, SPEC: 0b, LOG: 0b, YSUBTITLE: '[mV/m]', $
  DATA_ATT: data_att, COLORS: elim.colors, $
  LABELS: ['Ex', 'Ey', 'Ez'], LABFLAG: 1, $
  YTITLE: 'E_DSL (th' + sc +')'}
store_data, ename2[0], data=E, dlim=dlim

; STORE FAC DATA
if ~keyword_set(Efacname) then Efacname = Ename + '_clean_fac'
perp1 = 'E!DSP!N'
perp2 = 'E!Dperp!N'
para = 'E!D||!N'
data_att.coord_sys = 'fac: x in spin-plane'
dlim = {CDF: elim.cdf, SPEC: 0b, LOG: 0b, YSUBTITLE: '[mV/m]', $
  DATA_ATT: data_att, COLORS: elim.colors, $
  LABELS: [perp1, perp2, para], LABFLAG: 1, $
  YTITLE: 'E_FAC (th' + sc +')'}
store_data, efacname[0], data=Efac, dlim=dlim

; Clean up.
thx = 'th' + sc + '_'
store_data, ['Eclip', 'Bclip', thx + 'fgh_fac_mat'], /delete
end

