;+
; 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 spd_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 spd_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 spd_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 spd_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 spd_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