;+ ; 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 ; ; 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 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 thm_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 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 thm_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 thm_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 thm_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_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 ; Flag for AC/DC coupled (needed for deconvolving instrument response) 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