thm_efi_clean_efp.pro (Documentation for /home/davin/idl/socware/)

;+
; PRO: THM_EFI_CLEAN_EFP
;
; PURPOSE:
;    FIXES AXIAL, REMOVES SPIKES, SPINTONE, SC POTENTIAL, AND OTHER EFP ERRORS.
;
;   This routine is designed for general use. It will take EFP data in DSL as
;   input, clean up the data, and transform the data into FAC and GSM, and
;   return the cleaned EFP data in DSL, FAC, and GSM. If the data is gathered in
;   the dayside, set keyword /SUBSOLAR to use keyword defaults specifically for
;   dayside data.
;   
;   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.
;
;   CAUTION: The input keyword ENAME is optional. When it is not specified,
;       ENAME is defaulted to 'thx_efp', where 'x' is probe name and should be
;       one of these: 'a', 'b', 'c', 'd', and 'e'. In this case, if the existed
;       'thx_efp' is not in DSL, the output will be wrong. So, make sure the
;       existed 'thx_efp' is in DSL, or specify ENAME to another tplot variable
;       which is EFP data in DSL.
;
;   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) Merge can exhibit pathological behavior if B is near the spin plane.
;         spike removal does not always work.
;         3) 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_cleanefp.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, DEFAULT = 'thX_efp' (will
;               fetch)
;    Vname      OPTIONAL: Valid TPLOT name for E, DEFAULT = 'thX_vap' (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!!!! 
;    subsolar   OPTIONAL. If set, defaults for subsolar region are used.
;               DEFAULT = TAIL
;    talk       OPTIONAL. Plots diagnostics. Can be fun. Slow. DEFAULT = 0
;    duration_threshold:  OPTIONAL. If the duration of a burst is less than
;               duration_threshold, it is skipped. DEFAULT = 15. 
;               In units of seconds.
;
; KEYWORDS (for output tplot variables):
;    Edslname   OPTIONAL: Name for the cleaned efp in DSL. 
;               DEFAULT = Ename + '_clean_dsl'
;    Egsmname   OPTIONAL: Name for the cleaned efp in DSL. 
;               DEFAULT = Ename + '_clean_gsm'
;    Efacname   OPTIONAL: Name for the cleaned efp in DSL. 
;               DEFAULT = Ename + '_clean_fac'
;
; KEYWORDS (Remove_SpinTone):
;    SpinRemove     If entered as 0, suppresses spintone removal. DEFAULT = 4
;                   (Spintone removed). The following are the behaviors for
;                   different values of SpinRemove.
;                   1: Only spintone, 2: Spintone and 2nd harm. 3: Spintone and
;                   4th harm.  4: Spintone, 2nd, and 4th harm.
;    SpinNsmooth    Activates median smoothing. # of half periods. Must be odd.
;                   DEFAULT = 19 for subsolar and 39 otherwise.
;    SpinPoly       Activates polyfit smoothing. Can be unstable if >9. DEFAULT
;                   = 0
; 
; KEYWORDS (Remove_Potential):
;    VscRemove      If entered as 0, suppresses potential removal. DEFAULT = 1
;                   (Potential removed.)
;    VscPole        Frequency to median smooth Vsc before applying to Ez.
;                   DEFAULT = 5.0 Hz
;    Vpoly          Actives polyfit of Ez to Vsc.  DEFAULT = 2 (5 /Subsolar)
;    Use_Electrons  (OBSOLETE) Kept for backward compatibility.
;    TeName         (OBSOLETE) Kept for backward compatibility.
;
; 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
;
; KEYWORDS (Remove_Spin_Epoch):
;    EpochRemove    If set, removes all spin-epoch signals. DEFAULT = 0 (Not
;                   used.) Only useful for short stretches or if plasma
;                   condtitions are constant.
;
; KEYWORDS (Fix_Axial)
;    FixAxial       If entered as 0, suppresses axial fix. DEFAULT=1 (Axial is
;                   fixed.)
;    AxPoly         Forces polynomial fit of Ez-Eder of order AxPoly. DEFAULT=2
;                   (9 /SubSolar)
;    Merge          Merges Ederived with Eaxial.  DEFAULT = 0 No merging 
;                   (1 /Subsolar)
;                   1 - Abrupt merging. If Eder is avaliable Eder is used 
;                   (0 - Fmerge). Else Ez.
;                   2 - Soft merging. Bz/|B| is considered when merging.
;                   WARNING: MERGE=1 OR 2 CAN EXHIBIT PATHOLOGICAL BEHAVIOR IF
;                       MAGNETIC FIELD IS NEAR THE SPIN PLANE. PLEASE DOUBLE
;                       CHECK RESULTS IF USED. 
;    Fmerge         Frequency of crossover for merging Ederived with Eaxial. 
;                   DEFAULT = 5.0 Hz (/Subsolar). 
;    MergeRatio     Mimimum value of Bz/|B| to start merge (USED IF MERGE=2).
;                   DEFAULT = 0.05
;    MinBz          Needed to avoid divide by zero in Ederive. DEFAULT = 0.01 nT
;    MinRatio       Mimimum value of Bz/|B| to calculate Ederive. DEFAULT = 0.1
;                   (0.05 /Subsolar)
;
; KEYWORDS (B_Smooth)
;    BspinPoly      Activates polyfit smoothing. Can be unstable if >9. DEFAULT
;                   = 2. -1 suppresses.
;    Bsmooth        Activates B smoothing for fac rotation. DEFAULT = 11. 1 or
;                   less suppresses.
;
; HISTORY:
;   2013-06-20: JBT. Added deconvolution of sheath response.
;   2011-05-21: JBT. Added keyword DURATION_THRESHOLD.
;   2011-05-20: JBT. Fixed a bug that happens when the spin-tone removing fails.
;   2010-09-13: JBT. The default of FixAxial was changed to 0 (not to fix
;               axial E-field component).
;               Updated documents.
;   2010-04-08: JBT. Fixed some bugs.
;   2009-06-04: REE/JBT. Second Release.
;   2009-06-04: Jianbao Tao (JBT). The metadata (dlim.data_att) was improved.
;   2009-05-05: REE. First Release.
;-

function thm_efi_clean_efp_deconvol_inst_resp, E_array, srate, boom_type
; compile_opt idl2, hidden

Exf = E_array
nt_burst = n_elements(Exf)

fsample = srate
kernel_length = 512L
df = fsample / double(kernel_length)
f = dindgen(kernel_length)*df
f[kernel_length/2+1:*] -= double(kernel_length) * df

thm_comp_efi_response, boom_type, f, boom_resp, rsheath=5d6, /complex_response
; thm_comp_efi_response, 'SPB', f, SPB_resp, rsheath=5d6, /complex_response
; thm_comp_efi_response, 'AXB', f, AXB_resp, rsheath=5d6, /complex_response

E12_resp = 1 / boom_resp

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

; 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)
if nindx ne nt_burst then Exf = interpol(Exf[indx], t[indx], t)

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

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

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

return, Exf
end


;-------------------------------------------------------------------------------
pro thm_efi_clean_efp, probe=probe, Ename=Ename, Vname=Vname, $ ; GENERAL
        Bdslname=Bdslname, $                            ; GENERAL
        trange=trange, talk=talk, $                     ; GENERAL
        subsolar=subsolar, $                            ; GENERAL
        duration_threshold = duration_threshold, $      ; GENERAL
        Edslname = Edslname, Egsmname = Egsmname, $     ; OUTPUT NAME
        Efacname = Efacname, $                          ; OUTPUT NAME
        SpinNsmooth=SpinNsmooth, SpinPoly=SpinPoly, $   ; REMOVE_SPINTONE
        SpinRemove=SpinRemove, $                        ; REMOVE_SPINTONE
        VscPole=VscPole, Vpoly=Vpoly, $                 ; REMOVE_POTENTIAL
        VscRemove=VscRemove, $                          ; REMOVE_POTENTIAL
        use_electrons=use_electrons, TeName=TeName, $   ; REMOVE_POTENTIAL
        SpikeRemove=SpikeRemove, SpikeNwin=SpikeNwin, $ ; REMOVE_SPIKES
        SpikeSig=SpikeSig, SpikeSigmin=SpikeSigmin, $   ; REMOVE_SPIKES
        SpikeNfit=SpikeNfit, SpikeFit=SpikeFit, $       ; REMOVE_SPIKES
        EpochRemove=EpochRemove, $                      ; REMOVE_SPIN_EPOCH
        FixAxial=FixAxial, AxPoly=AxPoly, $             ; FIX_AXIAL
        Merge=Merge, FMerge=FMerge, $                   ; FIX_AXIAL
        MergeRatio=MergeRatio, $                        ; FIX_AXIAL
        MinBz=MinBz, MinRatio=MinRatio, $               ; FIX_AXIAL (E_DERIVE)
        BspinPoly=BspinPoly, Bsmooth=Bsmooth            ; B_SMOOTH

; # CHECK INPUTS - GET NEEDED DATA #
IF not keyword_set(probe) then BEGIN
  dprint, 'SC not set. Exiting...'
  return
ENDIF
sc = probe(0)

; CHECK FOR SUBSOLAR KEYWORD
IF keyword_set(subsolar) then BEGIN
  if n_elements(Vpoly) EQ 0         then Vpoly         = 4    ; REMOVE_POTENTIAL
  if n_elements(use_electrons) EQ 0 then use_electrons = 0    ; REMOVE_POTENTIAL
  if n_elements(Axpoly) EQ 0        then Axpoly        = 7      ; FIX_AXIAL
  if n_elements(Merge) EQ 0         then Merge         = 1      ; FIX_AXIAL
  if n_elements(MinRatio) EQ 0      then MinRatio      = 0.05d  ; FIX_AXIAL
  if n_elements(SpinNsmooth) EQ 0   then SpinNsmooth   = 19   ; REMOVE_SPINTONE
ENDIF

; SET TAIL DEFAULTS
if n_elements(SpinRemove) EQ 0  then SpinRemove  = 4      ; REMOVE_SPINTONE
if n_elements(SpinNsmooth) EQ 0 then SpinNsmooth = 39     ; REMOVE_SPINTONE
if n_elements(VscRemove) EQ 0   then VscRemove   = 1      ; REMOVE_POTENTIAL
if n_elements(VscPole) EQ 0     then VscPole     = 5.0d   ; REMOVE_POTENTIAL
if n_elements(Vpoly) EQ 0       then Vpoly       = 2      ; REMOVE_POTENTIAL
if n_elements(SpikeRemove) EQ 0 then SpikeRemove = 1      ; REMOVE_SPIKES
if n_elements(FixAxial) EQ 0    then FixAxial    = 0      ; FIX_AXIAL
if n_elements(Axpoly) EQ 0      then Axpoly      = 2      ; FIX_AXIAL
if keyword_set(merge)           then Fmerge      = 5.0d   ; FIX_AXIAL
if not keyword_set(merge)       then Fmerge      = 0      ; FIX_AXIAL
IF keyword_set(merge) then BEGIN
  if merge EQ 2 then softmerge=1
ENDIF
if n_elements(BspinPoly) EQ 0   then BspinPoly   = 2      ; B_SMOOTH
if n_elements(Bsmooth) EQ 0     then Bsmooth     = 11     ; B_SMOOTH
if n_elements(duration_threshold) eq 0 then duration_threshold = 15.

; CHECK FOR EFP DATA
if not keyword_set(Ename) then Ename = 'th' + sc + '_efp'
IF thm_check_tvar(Ename) then BEGIN                       ; added by JBT
  get_data, ename(0), data=E, dlim=elim                   ; added by JBT
  ; Check coordinates of Ename.
  if ~strcmp(elim.data_att.coord_sys, 'dsl', /fold) then begin
    dprint, Ename + ' is not in DSL. Exiting...'
    return
  endif
ENDIF ELSE BEGIN
  thm_load_efi, probe=sc, datatype=['efp', 'vap'], coord='dsl', trange=trange
  get_data, ename(0), data=E, dlim=elim
  IF size(/type,E) NE 8 then BEGIN
    dprint, 'Cannot get electric field data. Exiting...'
   return
  ENDIF
ENDELSE

; CHECK FOR VOLTAGES
if ~keyword_set(Vname) then Vname = 'th' + sc + '_vap'
if ~keyword_set(Vscname) then Vscname = 'th' + sc + '_vsc'
IF thm_check_tvar(Vscname) then BEGIN                          
  get_data, Vscname[0], data=Vsc                                  
ENDIF ELSE BEGIN                                               
  thm_efi_get_potential, Vname(0), trange=trange
  get_data, Vscname[0], data=Vsc
    IF size(/type,Vsc) NE 8 then BEGIN
      dprint, 'Cannot get SC potential. Exiting...'
    return
  ENDIF
ENDELSE

; CHECK FOR MAG DATA
if not keyword_set(Bdslname) then Bdslname = 'th' + sc + '_fgh_dsl'
IF thm_check_tvar(Bdslname) then BEGIN                         
  get_data, Bdslname(0), data=Bdsl, dlim=blim                  
ENDIF ELSE BEGIN                                               
  dprint, '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
    dprint, 'Cannot get MAG data. Exiting...'
    return
  ENDIF
ENDELSE

; CHECK FOR SPIN DATA
SpinName = 'th' + sc + '_state_spinper'                    
IF thm_check_tvar(Spinname) then BEGIN  
 get_data, SpinName[0], data=SpinPer                          
ENDIF ELSE BEGIN                                               
  thm_load_state, probe=sc, datatype='spinper'           
  get_data, SpinName[0], data=SpinPer   
  IF size(/type,SpinPer) NE 8 THEN BEGIN
    dprint, 'Cannot get spin period. Exiting...'
    return
  ENDIF
ENDELSE                                                  

; GET ELECTRON TEMPERATURE DATA 
;IF keyword_set(use_electrons) then BEGIN
;  if not keyword_set(TeName) then TeName = 'th' + sc + '_peeb_avgtemp'
;  IF thm_check_tvar(TeName) then BEGIN                          
;    get_data, TeName[0], data=peeb_t       
;  ENDIF ELSE BEGIN                                               
;    print, 'THM_EFI_CLEAN_EFP: Electron temperature data not stored. Fetching...'
;    thm_load_esa, probe=sc, datatype = 'peeb_avgtemp', level = 'l2', $
;      /get_support, trange = trange
;    TeName = 'th' + sc + '_peeb_avgtemp'
;    get_data, TeName[0], data=peeb_t
;    IF size(/type,peeb_t) NE 8 then BEGIN
;      print, 'THM_EFI_CLEAN_EFP: Cannot get electron temperature. Not using.'
;      use_electrons = 0
;    ENDIF
;  ENDELSE                                                  
;ENDIF
use_electrons = 0

; For making trange_clip work.
E = {x:E.x, y:E.y}
vsc = {x:vsc.x, y:vsc.y}
Bdsl = {x:Bdsl.x, y:Bdsl.y}

; CLIP DATA TO RANGE
IF keyword_set(trange) then BEGIN
   trange_clip, E, trange(0), trange(1), /data, BadClip=BadEclip
   trange_clip, Vsc, trange(0), trange(1), /data, BadClip=BadVclip
   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 size(/type,peeb_t) eq 8 then $
;     trange_clip, peeb_t, trange(0), trange(1), /data, BadClip=BadPclip
   IF (keyword_set(BadEclip) OR keyword_set(BadVclip) OR $
       keyword_set(BadBclip) OR $
       keyword_set(BadSclip) OR keyword_set(BadPclip) ) THEN BEGIN
     dprint, 'Problem with trange clip. Exiting...'
     dprint, '0=OK; 1=Problem. E:', BadEclip, 'V:', BadVclip, 'B:', $
        BadBclip, 'Spin:', BadSclip
     return
   ENDIF
ENDIF

; CREATE ARRAYS FOR OUTPUT
EderSave  = E.y(*,0) * 0

; ## IDENTIFY INDIVIDUAL PARTICLE BURSTS  ##
tE   = E.x
thm_lsp_find_burst, E, istart=bstart, iend=bend, nbursts=nbursts, mdt=mdt

; sample rate
srate = round(1d / median(tE[1:*] - tE)) + 0d
; print, 'srate = ', srate
; return
 
; START LOOP OVER INDIVIDUAL BURSTS
FOR ib=0L, nbursts-1 DO BEGIN
  dprint, 'BURST: ', ib+1, ' out of: ', nbursts

  ; 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]  
  
  ; CHECK THE TEMPORAL LENGTH OF THE BURST. IF IT IS SHORTER THAN 15 S, SKIP IT.
  blen  = tE[bend[ib]] - tE[bstart[ib]]
  dprint, 'burst temporal length = ', blen
  if blen lt duration_threshold then begin
      dprint, 'Burst #' + string(ib + 1, format='(I0)') + $
          ' is shorter than ' + string(duration_threshold, form = '(I0)') + $
          ' second ' + $
         'and is skipped from cleaning.'
      EderSave(bstart(ib):bend(ib))  = !values.d_nan
      continue
  endif

  ; Remove EFI boom response
  Ex = thm_efi_clean_efp_deconvol_inst_resp(Ex, srate, 'SPB')
  Ey = thm_efi_clean_efp_deconvol_inst_resp(Ey, srate, 'SPB')
  Ez = thm_efi_clean_efp_deconvol_inst_resp(Ez, srate, 'AXB')

  ; 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
    dprint, 'Spin period missing during burst. Exiting...'
    return
  ENDIF
  per = median(spinper.y(ind))
   
  ; GET SC POTENTIAL
  VscTemp = Vsc
  trange_clip, VscTemp, Ttemp(0), Ttemp(1), /data

  ; Remove spin tone in Vsc. If the removing fails, skip the removing.
  VscF   = thm_lsp_remove_spintone(VscTemp.x, VscTemp.y, per, $
      talk=talk, ns=SpinNsmooth, sp = spinpoly, fail = fail)
  if fail gt 0 then VscF = VscTemp.y

  VscF_old = VscF
  VscF   = thm_lsp_remove_spintone(VscTemp.x, VscF, per/2,    $
      talk=talk, ns=SpinNsmooth, sp = spinpoly, fail = fail)
  if fail gt 0 then VscF = VscF_old

  VscF_old = VscF
  VscF   = thm_lsp_remove_spintone(VscTemp.x, VscF, per/4,    $
      talk=talk, ns=SpinNsmooth, sp = spinpoly, fail = fail)
  if fail gt 0 then VscF = VscF_old

  ; REMOVE SPIN TONES FROM E
  IF keyword_set(SpinRemove) then BEGIN
    Exf = thm_lsp_remove_spintone(t, Ex, per,    $
        talk=talk, ns=SpinNsmooth, sp=spinpoly, fail = fail)
    if fail gt 0 then Exf = Ex

    Eyf = thm_lsp_remove_spintone(t, Ey, per,    $
        talk=talk, ns=SpinNsmooth, sp=spinpoly, fail = fail)
    if fail gt 0 then Eyf = Ey

    Ezf = thm_lsp_remove_spintone(t, Ez, per,    $
        talk=talk, ns=SpinNsmooth, sp=spinpoly, fail = fail)
    if fail gt 0 then Ezf = Ez

    IF ( (SpinRemove EQ 2) OR (SpinRemove GE 4) ) then BEGIN
      Exf_old = Exf
      Exf = thm_lsp_remove_spintone(t, Exf, per/2, $
          talk=talk, ns=SpinNsmooth, sp=spinpoly, fail = fail)
      if fail gt 0 then Exf = Exf_old
      
      Eyf_old = Eyf
      Eyf = thm_lsp_remove_spintone(t, Eyf, per/2, $
          talk=talk, ns=SpinNsmooth, sp=spinpoly, fail = fail)
      if fail gt 0 then Eyf = Eyf_old

      Ezf_old = Ezf
      Ezf = thm_lsp_remove_spintone(t, Ezf, per/2, $
          talk=talk, ns=SpinNsmooth, sp=spinpoly, fail = fail)
      if fail gt 0 then Ezf = Ezf_old
    ENDIF

    IF (SpinRemove GE 3) then BEGIN
      Exf_old = Exf
      Exf = thm_lsp_remove_spintone(t, Exf, per/4, $
          talk=talk, ns=SpinNsmooth, sp=spinpoly, fail = fail)
      if fail gt 0 then Exf = Exf_old

      Eyf_old = Eyf
      Eyf = thm_lsp_remove_spintone(t, Eyf, per/4, $
          talk=talk, ns=SpinNsmooth, sp=spinpoly, fail = fail)
      if fail gt 0 then Eyf = Eyf_old

      Ezf_old = Ezf
      Ezf = thm_lsp_remove_spintone(t, Ezf, per/4, $
          talk=talk, ns=SpinNsmooth, sp=spinpoly, fail = fail)
      if fail gt 0 then Ezf = Ezf_old
    ENDIF
  ENDIF

  ; REMOVE POTENTIAL FROM Ez
  if keyword_set(VscRemove) then $
    Ezf = thm_lsp_remove_potential(t, Ezf, VscTemp.x, VscF, peeb_t=peeb_t, $
                                   VscPole=VscPole, Vpoly=Vpoly, talk=talk)

  ; REMOVE SPIKES
  IF keyword_set(SpikeRemove) then BEGIN
    thm_lsp_remove_spikes, t, Exf, Eyf, Ezf, per, Nwin=SpikeNwin, $
        SpikeSig=SpikeSig, Sigmin=SpikeSigmin, Nfit=SpikeNfit, Fit=SpikeFit, $
        talk=talk, diagnose=diagnose, wt=wt
  ENDIF
  
  ; REMOVE SPIN TONES FROM E
  IF keyword_set(SpinRemove) then BEGIN
    Exf_old = Exf
    Exf = thm_lsp_remove_spintone(t, Exf, per,    $
        talk=talk, ns=SpinNsmooth, sp=spinpoly, fail = fail)
    if fail gt 0 then Exf = Exf_old

    Eyf_old = Eyf
    Eyf = thm_lsp_remove_spintone(t, Eyf, per,    $
        talk=talk, ns=SpinNsmooth, sp=spinpoly, fail = fail)
    if fail gt 0 then Eyf = Eyf_old

    Ezf_old = Ezf
    Ezf = thm_lsp_remove_spintone(t, Ezf, per,    $
        talk=talk, ns=SpinNsmooth, sp=spinpoly, fail = fail)
    if fail gt 0 then Ezf = Ezf_old

    IF ( (SpinRemove EQ 2) OR (SpinRemove GE 4) ) then BEGIN
      Exf_old = Exf
      Exf = thm_lsp_remove_spintone(t, Exf, per/2,    $
          talk=talk, ns=SpinNsmooth, sp=spinpoly, fail = fail)
      if fail gt 0 then Exf = Exf_old

      Eyf_old = Eyf
      Eyf = thm_lsp_remove_spintone(t, Eyf, per/2,    $
          talk=talk, ns=SpinNsmooth, sp=spinpoly, fail = fail)
      if fail gt 0 then Eyf = Eyf_old

      Ezf_old = Ezf
      Ezf = thm_lsp_remove_spintone(t, Ezf, per/2,    $
          talk=talk, ns=SpinNsmooth, sp=spinpoly, fail = fail)
      if fail gt 0 then Ezf = Ezf_old

;       Exf = thm_lsp_remove_spintone(t, Exf, per/2, $
;           talk=talk, ns=SpinNsmooth, sp=spinpoly)
;       Eyf = thm_lsp_remove_spintone(t, Eyf, per/2, $
;           talk=talk, ns=SpinNsmooth, sp=spinpoly)
;       Ezf = thm_lsp_remove_spintone(t, Ezf, per/2, $
;           talk=talk, ns=SpinNsmooth, sp=spinpoly)
    ENDIF
    IF (SpinRemove GE 3) then BEGIN
      Exf_old = Exf
      Exf = thm_lsp_remove_spintone(t, Exf, per/4,    $
          talk=talk, ns=SpinNsmooth, sp=spinpoly, fail = fail)
      if fail gt 0 then Exf = Exf_old

      Eyf_old = Eyf
      Eyf = thm_lsp_remove_spintone(t, Eyf, per/4,    $
          talk=talk, ns=SpinNsmooth, sp=spinpoly, fail = fail)
      if fail gt 0 then Eyf = Eyf_old

      Ezf_old = Ezf
      Ezf = thm_lsp_remove_spintone(t, Ezf, per/4,    $
          talk=talk, ns=SpinNsmooth, sp=spinpoly, fail = fail)
      if fail gt 0 then Ezf = Ezf_old
;       Exf = thm_lsp_remove_spintone(t, Exf, per/4, $
;           talk=talk, ns=SpinNsmooth, sp=spinpoly)
;       Eyf = thm_lsp_remove_spintone(t, Eyf, per/4, $
;           talk=talk, ns=SpinNsmooth, sp=spinpoly)
;       Ezf = thm_lsp_remove_spintone(t, Ezf, per/4, $
;           talk=talk, ns=SpinNsmooth, sp=spinpoly)
    ENDIF
  ENDIF

  ; REMOVE EPOCH
  IF keyword_set(EpochRemove) then BEGIN
    Exf = thm_lsp_remove_spin_epoch(t, Exf, per, talk=talk)
    Eyf = thm_lsp_remove_spin_epoch(t, Eyf, per, talk=talk)
    Ezf = thm_lsp_remove_spin_epoch(t, Ezf, per, talk=talk)
  ENDIF

  Btemp = Bdsl

  ; FIT AXIAL TO EDERIVED
  IF keyword_set(FixAxial) then BEGIN
    Etemp = E
    trange_clip, Etemp, Ttemp(0), Ttemp(1), /data
    trange_clip, Btemp, Ttemp(0), Ttemp(1), /data
    Etemp.y(*,0) = Exf
    Etemp.y(*,1) = Eyf
    Etemp.y(*,2) = Ezf
    Eder = thm_lsp_derive_Ez(Etemp, Btemp, minBz=minBz, $
        minRat=MinRatio, ratio=ratio)
    Ezf  = thm_lsp_fix_axial(t, Ezf, Eder, talk=talk, $
        AxPoly=Axpoly, soft=softmerge,  $
        Fmerge=Fmerge, MergeRatio=MergeRatio, ratio=ratio)
    EderSave(bstart(ib):bend(ib))  = Eder
  ENDIF 
  
  ; 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

  ; FIX B SPINTONE
  magstart = value_locate(Bdsl.x, Btemp.x(0)+1.d-8)
  magstop  = magstart + n_elements(Btemp.x) - 1L
  IF BspinPoly GE 0 then BEGIN
;     Btemp.y(*,0) 
    tmp = thm_lsp_remove_spintone(Btemp.x, Btemp.y(*,0), per, $
        talk=talk, spinpoly=Bspinpoly, fail = fail)
    if fail eq 0 then Btemp.y[*, 0] = tmp

    tmp = thm_lsp_remove_spintone(Btemp.x, Btemp.y(*,1), per, $
        talk=talk, spinpoly=Bspinpoly, fail = fail)
    if fail eq 0 then Btemp.y[*, 1] = tmp

    tmp = thm_lsp_remove_spintone(Btemp.x, Btemp.y(*,2), per, $
        talk=talk, spinpoly=Bspinpoly, fail = fail)
    if fail eq 0 then Btemp.y[*, 2] = tmp

;     Btemp.y(*,0) = thm_lsp_remove_spintone(Btemp.x, Btemp.y(*,0), per, $
;         talk=talk, spinpoly=Bspinpoly)
;     Btemp.y(*,1) = thm_lsp_remove_spintone(Btemp.x, Btemp.y(*,1), per, $
;         talk=talk, spinpoly=Bspinpoly)
;     Btemp.y(*,2) = thm_lsp_remove_spintone(Btemp.x, Btemp.y(*,2), per, $
;         talk=talk, spinpoly=Bspinpoly)
  ENDIF
  
  ; SMOOTH B
  IF Bsmooth GT 1 THEN BEGIN
    Btemp.y(*,0) = smooth(Btemp.y(*,0), Bsmooth, /nan, /edge_trunc)
    Btemp.y(*,1) = smooth(Btemp.y(*,1), Bsmooth, /nan, /edge_trunc)
    Btemp.y(*,2) = smooth(Btemp.y(*,2), Bsmooth, /nan, /edge_trunc)
  ENDIF
  
  ; SAVE B DATA
  Bdsl.y(magstart:magstop, 0) =  Btemp.y(*,0)
  Bdsl.y(magstart:magstop, 1) =  Btemp.y(*,1)
  Bdsl.y(magstart:magstop, 2) =  Btemp.y(*,2)

ENDFOR
; ## END OF LOOP

; ### STORE DATA AS TPLOT VARIABLES

; add BAND to data_att -JBT
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: 'DC - ~50 Hz'}   ; BAND - the freq band of the data

; STORE E 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: elim.labflag, $
  YTITLE: 'E_DSL (th' + sc +')'}
store_data, ename2[0], data=E, dlim=dlim

; STORE B DATA
bname2 = Bdslname + '_smooth'
store_data, bname2[0], data=Bdsl, dlim=blim

; STORE Ederived DATA
if keyword_set(FixAxial) then begin
  Edername = ename + '_derived'
  Etemp = dblarr(n_elements(E.x),4)
  Etemp(*,0) = E.y(*,0)
  Etemp(*,1) = E.y(*,1)
  Etemp(*,2) = E.y(*,2)
  Etemp(*,3) = EderSave
  dlim = {CDF: elim.cdf, SPEC: 0b, LOG: 0b, YSUBTITLE: '(mV/m)', $
    DATA_ATT: data_att, COLORS: [elim.colors, 0], $
    LABELS: ['Ex', 'Ey', 'Ez', 'E!Dzder!N'], LABFLAG: elim.labflag, $
    YTITLE: 'E - ' + data_att.coord_sys}
  store_data, Edername[0], data={X:E.x, Y: Etemp, V: [1,2,3,4]}, dlim=dlim
endif

; #### COORDINATE TRANSFORMATIONS

; TRANSFORM TO GSM
if ~keyword_set(Egsmname) then Egsmname = Ename + '_clean_gsm'
thm_cotrans, ename2, egsmname, in_coord='dsl', out_coord='gsm'
get_data, egsmname[0], data = data
data_att.coord_sys = 'gsm'
dlim = {CDF: elim.cdf, SPEC: 0b, LOG: 0b, YSUBTITLE: '(mV/m)', $
  DATA_ATT: data_att, COLORS: elim.colors, $
  LABELS: ['Ex', 'Ey', 'Ez'], LABFLAG: elim.labflag, $
  YTITLE: 'E_GSM (th' + sc +')'}
store_data, egsmname[0], data = data, dlim = dlim

; GO TO FAC COORDINATES (JIANBAO)
if ~keyword_set(Efacname) then Efacname = Ename + '_clean_fac'
thm_lsp_clean_timestamp, bname2
thm_lsp_clean_timestamp, ename2
thm_fac_matrix_make, bname2, other_dim='zdsl', newname='th'+sc+'_fgh_fac_mat'
tvector_rotate, 'th'+sc+'_fgh_fac_mat', ename2, $
          newname=efacname, error=error
get_data, efacname[0], data = data
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: elim.labflag, $
  YTITLE: 'E_FAC (th' + sc +')'}
store_data, efacname[0], data = data, dlim = dlim

; Clean up.
thx = 'th' + sc + '_'
store_data, [bname2[0], thx + 'fgh_fac_mat'], /del

end