;+
;
;function THM_GET_EFI_EDC_OFFSET()
;
;Purpose:
; Estimate the EDC offset by low-pass filtering the ADC data with a moving estimation window. Returns -1 if
;
;Syntax:
; result = thm_get_efi_edc_offset( Y, Sample_Rate, N_Spins, Spin_Period [, Offset_Estimation_Window_Truncated ] [, NEW_NSPINS ] )
;
; where
;
; Y (n-element array) Input The ADC data needing the EDC offset calculation.
; Sample_Rate Input The data rate.
; N_Spins Input The width of the estimation window in # of spins
; Spin_Period Input The period of a spin [s].
; Edge_Truncate Input Passed to SMOOTH, if SMOOTH is the estimation function.
; Offset_Estimation_Window_Truncated Output Indicates that the data interval is shorter than the estimation window requested, so the estimation window
; has been truncated to the greatest integral number of spins that fit into the data interval.
; New_Nspins Output If OFFSET_ESTIMATION_WINDOW_TRUNCATED is set, then this will contain the adjusted number of spins.
;
;-
function thm_get_efi_edc_offset, y, samp_rate, nspins, spin_period, edge_truncate, offset_estimation_window_truncated, nspins_temp
compile_opt idl2
; Initialize outputs:
;
undefine, nspins_temp
offset_estimation_window_truncated = 0b
npts = ceil( float(nspins)*spin_period*samp_rate ) ; width of smoothing window in # of points.
n_y = n_elements( y )
if npts ge n_y then begin
offset_estimation_window_truncated = 1b
npts = (n_y mod 2) ? n_y-2 : n_y-1 ; if chunk is less than NSPINS wide, then make NPTS just fit.
nspins_temp = floor(npts/spin_period/samp_rate) ; recalculate NSPINS (NSPINS_TEMP) as least whole # w/in interval.
npts = ceil( float(nspins_temp)*spin_period*samp_rate ) ; recalculate NPTS.
endif
if npts gt 0 then return, smooth( y, npts, /nan, edge_truncate = edge_truncate ) else begin
undefine, offset_estimation_window_truncated
undefine, nspins_temp
return, -1
endelse
end ;thm_get_edc_offset_efi
;+
;
;THM_UPDATE_EFI_LABELS.PRO
;
;PURPOSE:
;To be used on THEMIS EFI 3-vectors. Labels x-axis as "E12 (<coordinate system in DLIMITS.LABLES>)" for spinning coordinate systems and
;"Ex (<coord. sys. in DLIMITS.LABELS>)" for despun coordinate systmes. Y and Z axes are upadated correspondingly. Call after coordinate
;transformations.
;
;SYNTAX:
; thm_update_efi_labels ,<String>
;
;Arguments:
; <String>: A TPLOT variable name referring to a 3-vector TPLOT variable.
;
;Code: W. Michael Feuerstein, 7/2008.
;
;-
pro thm_update_efi_labels ,tvarname
compile_opt idl2, hidden
coord=cotrans_get_coord(tvarname)
if where(coord eq ['spg','ssl']) ge 0 then options,/default,tvarname,'labels',['E12','E34','E56']+' ('+coord+')' else $
options,/default,tvarname,'labels',['Ex','Ey','Ez']+' ('+coord+')'
end
;+
;Procedure: THM_CAL_EFI
;
;Purpose: Converts raw EFI (V, EDC, and EAC waveform) data into
; physical quantities.
;
;Syntax: THM_CAL_EFI [, <optional keywords below> ]
;
;keywords:
; /VERBOSE or VERBOSE=n ; set to enable diagnostic message output.
; higher values of n produce more and lower-level diagnostic
; messages.
; datatype - Default setting is to calibrate all raw quantites and also
; produce all _0 and _dot0 quantities that are passed in
; DATATYPE kw.
; probe - Defaults to all probes.
; valid_names - Return valid datatypes, print them, and return.
; coord - Coordinate system of output.
; TEST: Disables selected /CONTINUE to MESSAGE. For QA testing only.
; STORED_TNAMES: An OUTPUT. Returns a string array containing each TPLOT variable name
; invoked in a STORE_DATA operation (chronological order).
; (Not sorted or uniqued.)
; ONTHEFLYEDCOFFSET: Return the EDC offset array calculated on-the-fly *** FOR THE LAST EF? DATATYPE PROCESSED ONLY ***.
; GAP_TRIGGER_VALUE: An INPUT and OUTPUT. For on-the-fly EDC offset calculation, consider anything greater than or equal to GAP_TRIGGER_VALUE to be
; a gap in the data. Default: 0.5 s.
; NUMBER_OF_SPINS: I/O. Specify the number of spins for the on-the-fly EDC offset
; calculation estimation window, or read out the default (20 spins).
; OFFSET_ESTIMATION_FUNC: OUTPUT <string scalar>. The name of the function used to
; estimate the EDC offset for the on-the-fly window.
; /EDGE_TRUNCATE: I/O. Set to 0 to disable edge truncation in SMOOTH (for the
; on-the-fly offset calculation). Assign to a variable to read the
; default (= 1). Undefined, if on-the-fly offset not done.
;
;keyword parameters for _dot0 computation:
; max_angle - Maximum angle of B field to spin axis to calculate _dot0.
; Typical = 80 degrees. No default.
; min_bz - Minimum value of Bz. Typical value is 2.0 nT.
; Default= 1.0 nT. Not compatible with max_angle keyword.
; max_bxy_bz- Maximum value of abs(bx/bz) or abs(by/bz). Typical value is
; 5. ~= tan(79 degrees) (think of Bx/Bz).
; Default is not to use this method (no default value).
; bz_offset - Offset in nT that will be added to Z axis measurement of B.
; Defaults to 0.0.
; fgm_datatype - 'fgl', 'fgh', 'fgs' or 'fge'
;
;Example:
; thm_cal_efi, probe='c'
;
;Modifications:
; Added boom_shorting_factor_e12, boom_shorting_factor_e34,
; offset_e12, offset_e34, offset_dsc_x, offset_dsc_y
; fields to data_att structure in the default limits structure.
; Also, added mechanism to fill new field from
; THM_GET_EFI_CAL_PARS.PRO, W.M.Feuerstein, 2/26-27/2008.
; Changed "History" to "Modifications" in header, made "1000" a float ("."),
; switched all internal calibration parameters over to those read through
; THM_GET_EFI_CAL_PARS.PRO, WMF, 3/12/2008.
; Updated DATA_ATT fields to match revised cal. files, WMF, 3/12/2008.
; Updated error handling, updated doc'n, WMF, 3/12/2008.
; Updated doc'n, made datatypes more consistent (speed), WMF, 3/14/2008.
; Changed "no default" for MIN_BZ kw to default=1.0 nT, fixed MAX_ANGLE
; and MIN_BZ kwd's "both set" handling, WMF, 3/14/2008.
; Changed an NaN multiplication to an assignment (speed), simplified MIN_BZ
; default assignment, turned on good practice compile options,
; reorganized opening statements (flow), turned off ARG_PRESENT(MAX_BXY_BZ)
; bug (no default for this kw), updated doc'n, WMF, 3/17/2008.
; Now calculates Ez (_dot0 section) AFTER implementing MAX_BXY_BZ kw as per
; J. Bonnell's request (it also happens to be faster), updated doc'n,
; WMF, 3/18/2008.
; Removed redundant datatype conversions, WMF, 3/20/2008.
; Using BOOM_SHORTING_FACTOR in field calibration, WMF, 3/20/2008 (Th).
; Fixed potential logical vs. bitwise conditional bug, WMF, 3/21/2008 (F).
; Updated doc'n, WMF, 3/27/2008.
; Removed "TWEAK_GAINS" kw to THM_EFI_DESPIN.PRO per J.Bonnell's req.,
; WMF, 4/4/2008 (F).
; Added TEST kw to disable certain /CONTINUE to MESSAGE, and passed TEST
; through to THM_GET_EFI_CAL_PARS.PRO, WMF, 4/7/2008 (M).
; Implemented time-dependent EAD/EDC gain conditional, WMF, 4/22/2008 (Tu).
; Implemented time-dependent calibration parameters for v?? and
; e?? datatypes, WMF, 5/21 - 6/5/2008.
; Reconciled last three changes from non-time-dependent version: COORD kw is now not case
; sensitive, made sure 'E12','E34','E56' labels go with SSL and SPG coordinates and
; 'Ex','Ey','Ez' go with all other coordinate systems (THM_UPDATE_EFI_LABELS.PRO),
; fixed COORD='gse' crash, WMF, 8/21/2008.
; Added error message for hed_ac = 255, WMF, 8/26/2008.
; Renamed from "thm_cal_efi_td.pro" to "thm_cal_efi.pro", WMF, 9/9/2008.
; Put in cases for 2 and 4 spin-dependent offsets in the boom plane, WMF, 3/4/2009.
; Insert on NaN on FGM degap, memory management, WMF, 5/5/2009.
;
;Notes:
; -- fixed, nominal calibration pars used (gains and
; frequency responses), rather than proper time-dependent parameters.
;
; $LastChangedBy: michf $
; $LastChangedDate: 2008-04-22 17:50:30 -0700 (Tue, 22 Apr 2008) $
; $LastChangedRevision: 2804 $
; $URL $
;-
pro thm_cal_efi, probe=probe, datatype=datatype, files=files,trange=trange,$
valid_names=valid_names, verbose=verbose, coord=coord, $
in_suffix=in_suf,out_suffix=out_suf, max_angle=max_angle, $
min_bz=min_bz, max_bxy_bz=max_bxy_bz, bz_offset=bz_offset, $
fgm_datatype=fgm_datatype, test=test, $
stored_tnames = stored_tnames, $
ontheflyedcoffset = d_smooth, $
gap_trigger_value = gap_trigger_value, $
number_of_spins = nspins, $
offset_estimation_func = offset_estimation_func, $
edge_truncate = edge_truncate
compile_opt idl2, strictarrsubs ;Bring this routine up to date.
if ~keyword_set(test) then test = 0
;Check for "max_angle" and "min_bz" kwd's both being set:
;========================================================
;if keyword_set(max_angle) and keyword_set(min_bz) then begin ;Fails on var=0.
if size(max_angle,/type) ne 0 && size(min_bz,/type) ne 0 then begin ;This way.
message,/cont, $
'MIN_BZ and MAX_ANGLE keywords are mutually exclusive. Returning...'
return
endif
thm_init
vb = size(verbose, /type) ne 0 ? verbose : !themis.verbose
;Define valid EFI "datatypes":
;=============================
primary_datatypes = [ 'vaf', 'vap', 'vaw', 'vbf', 'vbp', 'vbw', 'eff', 'efp', 'efw'] ; note that efi_dq is deleted.
u0s=['eff_0','efp_0','efw_0']
udot0s=['eff_dot0','efp_dot0','efw_dot0']
efi_valid_names = [primary_datatypes, u0s, udot0s]
;Return valid datatypes (and return), if they are called for:
;============================================================
if arg_present( valid_names) then begin
valid_names = efi_valid_names
message, /info, string( strjoin( efi_valid_names, ','), format='( "Valid names:",X,A,".")')
return
endif
;Set COORD default to 'dsl':
;============================
if not keyword_set(coord) then coord='dsl' else coord=strlowcase(coord)
;Set in_suf/out_suf variables to null, if not set:
;=================================================
if not keyword_set(in_suf) then in_suf=''
if not keyword_set(out_suf) then out_suf=''
;Define "probes" variable ("f" is not a valid probe unless it is the only
;element). Return if "probes" does not get defined:
;===================================================
vprobes = ['a','b','c','d','e']
if n_elements(probe) eq 1 then if probe eq 'f' then vprobes = ['f']
if not keyword_set(probe) then probes = vprobes else probes = thm_check_valid_name(strlowcase(probe), vprobes, /include_all)
if not keyword_set(probes) then return
if keyword_set(vb) then printdat, probes, /value, varname='Probes'
;Define "dts" (datatypes) variable. Return if "dts" does not get defined:
;=========================================================================
if not keyword_set(datatype) then dts = efi_valid_names else dts = thm_check_valid_name(strlowcase(datatype), efi_valid_names, /include_all)
if not keyword_set(dts) then return
if keyword_set(vb) then printdat, dts, /value, varname='Datatypes'
;Make min_bz = 1.0 nT the default (as long as MAX_ANGLE is not defined!):
;========================================================================
if ~size(min_bz,/type) && ~size(max_angle,/type) then min_bz = 1.0 ;nT
;Set BZ_OFFSET default to 0.0:
;=============================
if ~size(bz_offset,/type) then bz_offset=0.0
;Loop on "probes":
;=================
for s=0L,n_elements(probes)-1L do begin
sc = probes[s]
;Loop on "dts":
;==============
for n=0L, n_elements(dts)-1L do begin
name = dts[ n]
;NAMERAW= 'aaa', NAME= (e.g.) 'aaa_dot0':
;========================================
if where(name eq u0s or name eq udot0s) ne -1 then nameraw=strmid(name,0,3) else nameraw=name
;*****************************************
;Define tplot var. names and get raw data:
;*****************************************
tplot_var_raw = thm_tplot_var( sc, nameraw)+in_suf
tplot_var = thm_tplot_var( sc, name)+out_suf
;
;Test for the presence of "veryunusualprefixtempfoo_" prefixed TPLOT vars. If not present, assume
;this call is not from THM_LOAD_EFI and look for normal (non-temporary) TPLOT vars to calibrate:
;***********************************************************************************************
if (tnames('veryunusualprefixtempfoo_*'))[0] ne '' then begin
get_data, 'veryunusualprefixtempfoo_'+tplot_var_raw, data=d, limit=l, dlim=dl
endif else begin
get_data, tplot_var_raw, data=d, limit=l, dlim=dl
endelse
if keyword_set( vb) then message, /cont, string( tplot_var, format='("Working on TPLOT variable",X,A)')
;Check to see if the data has already been calibrated:
;=====================================================
if( thm_data_calibrated( dl)) then begin
;If DATATYPE = first three letters of DATATYPE, then the data has
;been calibrated already:
;========================
if nameraw eq name then begin
if keyword_set( vb) then message, /cont, string( tplot_var, format= '(A,X,"has already been calibrated.")')
;Else it's gotta be either "_0", or "_dot0", so get cal_pars ("cp") and
;transform to 'dsl' (via THM_COTRANS for non-'dsl' and non-'spg', else via THM_EFI_DESPIN.PRO):
;==============================================================================================
endif else begin
thm_get_efi_cal_pars, d.x, nameraw, probes, cal_pars=cp,test=test
if where(cotrans_get_coord(dl) eq ['dsl','spg']) lt 0 then begin
thm_cotrans, tplot_var_raw,tplot_var,out_coord='dsl'
endif else begin
thm_efi_despin, sc, nameraw, cp.dsc_offset[0,*], [1,1,1], $ ;Pass [1,1,1] in order to not duplicate boom_shorting factor.
tplot_name=tplot_var_raw, newname=tplot_var, stored_tnames = foo
if ~(~size(stored_tnames,/type)) then stored_tnames = [ stored_tnames, foo ] else stored_tnames = foo
endelse
endelse
;Otherwise, proceed with the RAW->PHYS transformation:
;=====================================================
endif else begin
hed_tplot_var=thm_tplot_var(sc,nameraw)
tplot_var_hed = hed_tplot_var + '_hed'
get_data, tplot_var_hed, data=d_hed, limit=l_hed, dlim=dl_hed
;Check that returned data and hed structures are structures
;(get_data returns 0 if no TPLOT variable exists):
;=================================================
if (size( d, /type) eq 8) and (size( d_hed, /type) eq 8) then begin
hed_data = thm_unpack_hed( name, d_hed.y) ;Leaving this line for refernce, but product not used. 1/27/09.
;Match datatype by a switch statement, get corresponding
;calibration parameters, apply calibration, and update
;default limits structure:
;=========================
switch strlowcase( name) of
'vaf':
'vbf':
'vap':
'vbp':
'vaw':
'vbw': begin
;Grab calibration params, and initialize result:
;***********************************************
thm_get_efi_cal_pars, d.x, name, probes, cal_pars=cp, test=test
res = fltarr( size( d.y, /dimensions))
;Remember time-doubles and how many time-dependent intervals span data:
;**********************************************************************
cptd = time_double(cp.cal_par_time) ;"Calibration Parameter Time Double."
ntdi=n_elements(cp.cal_par_time) ;"# Time Dependent Indicies".
;Loop on time-dependent calibration intervals calibrating result for each interval:
;**********************************************************************************
for i=0,ntdi-1 do begin
;Get data time indicies for interval:
;************************************
if i ne ntdi-1 then begin
w = where(d.x ge cptd[i] and d.x lt cptd[i+1])
endif else begin
w = where(d.x ge cptd[i])
endelse
;Calibrate, looping over EFI booms:
;**********************************
for icomp=0L,5L do begin
if ntdi eq 1 then begin
res[w, icomp] = cp.gain[icomp]*(d.y[w, icomp] - cp.offset[icomp])
endif else begin
res[w, icomp] = cp.gain[i,icomp]*(d.y[w, icomp] - cp.offset[i,icomp]) ;This will work for both cases, now.
endelse
endfor
endfor
;; update the DLIMIT elements to reflect RAW->PHYS
;; transformation, coordinate system, etc.
units = cp.units ;Now taken from cal. file.
str_element, dl, 'data_att', data_att, success=has_data_att
if has_data_att then begin
str_element, data_att, 'data_type', 'calibrated', /add
endif else data_att = { data_type: 'calibrated' }
str_element, data_att, 'coord_sys', 'efi_sensor', /add
str_element, data_att, 'units', units[0], /add ;We cheat here. Only units of first interval are shown!
str_element, data_att, 'cal_par_time', cp.cal_par_time, /add
;"Documented calibration factors" to match expanded
;cal. file params. (3/12/2008):
;==============================
str_element,data_att,'offset',cp.offset,/add
str_element,data_att,'gain',cp.gain,/add
str_element,data_att,'boom_length',cp.boom_length,/add
str_element,data_att,'boom_shorting_factor', $
cp.boom_shorting_factor,/add
str_element, dl, 'data_att', data_att, /add
str_element, dl, 'ytitle', string(tplot_var, units[0], format= $ ;We cheat here. Only units of first interval are shown!
'(A,"!C!C[",A,"]")'), /add
str_element, dl, 'units', units[0], /add ;We cheat here. Only units of first interval are shown!
str_element, dl, 'labels', [ 'V1', 'V2', 'V3', 'V4', 'V5', 'V6'], /add
str_element, dl, 'labflag', 1, /add
str_element, dl, 'colors', [ 1, 2, 3, 4, 5, 6]
;; store the transformed spectra back into the original
;; TPLOT variable.
store_data, tplot_var, data = { x:d.x, y: temporary(res), v:d.v }, lim=l, dlim=dl
if ~(~size(stored_tnames,/type)) then stored_tnames = [ stored_tnames, tplot_var ] else stored_tnames = [ tplot_var ]
break
end ; end of { VAF, VBF, VAP, VBP, VAW, VBW} calibration clause.
'eff':
'efp':
'efw':
'efp_0':
'eff_0':
'efw_0':
'efp_dot0':
'efw_dot0':
'eff_dot0': begin ; { EFF, EFP, EFW} calibration clause.
; For Metadata/settings:
;
offset_estimation_func = 'SMOOTH' ; This may be user-selectable later.
if undefined(edge_truncate) then edge_truncate = 1b
offset_estimation_window_truncated = 0b ; Initialize to 0 for each datatype.
;These are boom lengths in meters:
;=================================
;e12=49.6 ;These are now taken from the cal. file.
;e34=40.4
;e56=5.63
;
;Grab calibration params, and initialize result:
;***********************************************
thm_get_efi_cal_pars, d.x, nameraw, probes, cal_pars=cp, test=test
res = fltarr( size( d.y, /dimensions))
;Remember time-doubles and how many time-dependent intervals span data:
;**********************************************************************
cptd = time_double(cp.cal_par_time) ;"Calibration Parameter Time Double."
ntdi=n_elements(cp.cal_par_time) ;"# Time Dependent Indicies".
;Loop on time-dependent calibration intervals calibrating result for each interval:
;**********************************************************************************
for i=0,ntdi-1 do begin
exx=cp.boom_length[i,*]*cp.boom_shorting_factor[i,*]
;Get data time indicies for interval:
;************************************
if i ne ntdi-1 then begin
w = where(d.x ge cptd[i] and d.x lt cptd[i+1])
endif else begin
w = where(d.x ge cptd[i])
endelse
;
;==================
;Calculate E field:
;==================
;
;================================================================
;Test for header_ac data. If not found, default to EAC
;gain and print warning. If found, set accordingly:
;(This step must be done w/in NTDI loop b/c the gain may change.)
;================================================================
tpname_hed_ac = 'th'+probes[0]+'_'+nameraw+'_hed_ac'
tpnames = tnames(tpname_hed_ac) ;Avoid known name length bug of TPLOT_NAMES.
if ~(~size(tpnames,/type)) && tpnames[0] ne '' then begin
get_data,tpnames[0],data=data_hed_ac
;
;Index AC vs. DC-coupled points in HED_AC:
;=========================================
n_hed_ac = n_elements(data_hed_ac.y)
n_hed_ac_2=n_hed_ac-2
if (where(data_hed_ac.y eq 255))[0] ge 0 then message, $
'Spurious "hed_ac" flag value of 255 encounterd. Contact THEMIS software team with this error message!'
wac = where(data_hed_ac.y,n_wac,complement = wdc)
;
;
case 1 of
n_wac eq n_hed_ac: gain = cp.eac_gain[i,*]
n_wac eq 0: gain = cp.edc_gain[i,*]
else: begin
;
;GAIN will be set to EAC_GAIN (if all points are 1), EDC_GAIN
;(if all points are 0), or interpolated to TIMES otherwise:
;==========================================================
aci=bytarr(n_elements(w)) ;"AC Interpolated".
acx = data_hed_ac.x ;Do not make structure ref. inside loop.
acy = data_hed_ac.y ;Do not make structure ref. inside loop.
;
;Interpolation of HED_AC data:
;=============================
dxw=(d.x)[w] ;Do not repeat array reference inside loop.
;j=0 ;Index to HED_AC data.
w0=where(acx le dxw[0])
j=0>w0[n_elements(w0)-1] ;Index of largest HED_AC index where acx le dxw[0] (get the loop started at the right spot).
for k=0,n_elements(w)-1 do begin
if dxw[k]-acx[j] ge acx[j+1]-dxw[k] then j = j ne n_hed_ac_2 ? ++j : j
aci[k] = acy[j]
endfor
;
;Multiply and merge:
;===================
gain = dblarr(n_elements(w),3)
gain[*,0] = aci*cp.eac_gain[i,0] + (~aci)*cp.edc_gain[i,0]
gain[*,1] = aci*cp.eac_gain[i,1] + (~aci)*cp.edc_gain[i,1]
gain[*,2] = aci*cp.eac_gain[i,2] + (~aci)*cp.edc_gain[i,2]
end
endcase
endif else begin
message,'*** WARNING!: TPLOT variable "'+ $
tpname_hed_ac+ $
'" not found. Defaulting to EDC gain.', $
continue=~test
gain=cp.edc_gain[i,*]
endelse
case (size(cp.dsc_offset,/dimensions))[1] of
5: begin
;***************************************
;4 spin-dependent offsets in boom plane:
;***************************************
;
message,/cont,'Using 4 spin-dependent offsets in boom plane.'
;
;Interpolate spin phase:
;=======================
;
model=spinmodel_get_ptr(probes[0])
spinmodel_interp_t,model=model,time=(d.x)[w], spinphase=phase,spinper=spinper ;a la J. L.
phase*=!dtor
;Calculate corrected (SPG) field:
;********************************
;
ex012=cp.dsc_offset[i,0]
ex034=cp.dsc_offset[i,1]
ey012=cp.dsc_offset[i,2]
ey034=cp.dsc_offset[i,3]
e012=cp.offset[i,0]
e034=cp.offset[i,1]
;GAIN is time-independent case:
;==============================
if (size(gain,/dimensions))[0] le 1 then begin
res[w,0]= -1000.*gain[0]/exx[0]*(d.y[w,0] - e012) - ex012*cos(phase) - ey012*sin(phase)
res[w,1]= -1000.*gain[1]/exx[1]*(d.y[w,1] - e034) + ex034*sin(phase) - ey034*cos(phase)
res[w,2]= -1000.*gain[2] * (d.y[w,2] - cp.offset[i,2])/exx[2] ;Z [mV/m]
;
;GAIN is an array case:
;======================
endif else begin
res[w,0]= -1000.*gain[*,0]/exx[0]*(d.y[w,0] - e012) - ex012*cos(phase) - ey012*sin(phase)
res[w,1]= -1000.*gain[*,1]/exx[1]*(d.y[w,1] - e034) + ex034*sin(phase) - ey034*cos(phase)
res[w,2]= -1000.*gain[*,2] * (d.y[w,2] - cp.offset[i,2])/exx[2] ;Z [mV/m]
endelse
end
3: begin
;***************************************
;2 spin-dependent offsets in boom plane:
;***************************************
;
message,/cont,'Using 2 spin-dependent offsets in boom plane.'
;Calibrate, looping over EFI booms:
;**********************************
;On-the-fly EDC offset calculation (smoothing or convolution w/ hanning window):
;*******************************************************************************
;
;
; Define # of spins and get spinmodel info:
;
if undefined(nspins) then begin
nspins = 20L ;width of smoothing window in spins.
endif else begin
case 1 of
is_num(nspins): nspins = nspins[0]
else: message, 'NUMBER_OF_SPINS keyword input must be scalar numeric. Please correct input and retry.'
endcase
endelse
;
; *** Need to insert J. L.'s code for checking the up-to-dateness of the state information here (when ready)! ***
;
model=spinmodel_get_ptr(probes[0])
spinmodel_interp_t,model=model,time=(d.x)[w], spinper=spin_period ;a la J. L.
spin_period = median(spin_period)
; Inform user of on-the-fly offset method:
;
message, /info, 'Using on-the-fly EDC offset removal with '+strtrim(nspins, 2)+ ' spin estimation window assuming a '+ $ ;smooth
strtrim(spin_period, 2) +' [s] spin period (median of loaded data spin periods).'
message, /info, 'Using ' + strtrim(offset_estimation_func,2)+' for offset estimation function. ' + $
'Edge truncation is '+(edge_truncate ? 'ON.':'OFF.')
; Define array for EDC offset calculation:
;
d_smooth = dblarr( n_elements(w), 2 ) ; Smoothing E12 and E34 only -- not E56.
; Find any gaps in the data (for E12 and E34 separately):
;
if undefined(gap_trigger_value) then gap_trigger_value = 0.5
xdegap, gap_trigger_value, 0., d.x[w], float(d.y[w,0:1]), x_out, y_out, iindices = ii, /onenanpergap, /nowarning
; Smooth each chunk of data individually (if there are N gaps, then there are N+1 chunks):
; *** Truncate NSPINS, if chunk is too short, and send warning MESSAGE. ***
; ****************************************************************************************
;
for j = 0,1 do begin ;0 for E12, 1 for E34.
if y_out[0] ne -1 then q = where( finite( y_out[*, j], /nan ) ) else q = -1
if q[0] ne -1 then begin
imap = lindgen(n_elements(ii)) ;Maps indices in x/y_out into d.y[w,*].
n_q = n_elements( q )
;
; 1st chunk:
;
samp_rate = 1./median( x_out[ 1:q[0]-1 ] - x_out[ 0:q[0]-2 ] ) ; samp/s.
wdii = where( ii lt q[0] )
wd = imap[ wdii ]
result = thm_get_efi_edc_offset( d.y[wd, j], samp_rate, nspins, spin_period, edge_truncate, window_truncated, new_nspins )
if result[0] ne -1 then begin
d_smooth[wd, j] = temporary(result)
if window_truncated then begin
offset_estimation_window_truncated = 1b
message, /info, $
' *** WARNING: smoothing window truncated to '+string(new_nspins,format='(i0)')+ ' spins to fit data sub-interval.'
endif
endif else message, /info, ' *** WARNING: Interval less than 1 spin wide. EDC offset not calculated and not subtracted!'
;
; Nth chunk:
;
if n_q ge 2 then begin ; Only do this, if there are at least two gaps.
for k = 0, n_elements(q)-2 do begin
samp_rate = 1./median( x_out[ q[k]+2:q[k+1]-1 ] - x_out[ q[k]+1:q[k+1]-2 ] )
wdii = where( ii gt q[k] and ii lt q[k+1] )
wd = imap[ wdii ]
result = thm_get_efi_edc_offset( d.y[wd, j], samp_rate, nspins, spin_period, edge_truncate, window_truncated, new_nspins )
if result[0] ne -1 then begin
d_smooth[wd, j] = temporary(result)
if window_truncated then begin
offset_estimation_window_truncated = 1b
message, /info, $
' *** WARNING: smoothing window truncated to '+string(new_nspins,format='(i0)')+ ' spins to fit data sub-interval.'
endif
endif else message, /info, ' *** WARNING: Interval less than 1 spin wide. EDC offset not calculated and not subtracted!'
endfor ;k
endif
;
; Last chunk:
;
samp_rate = 1./median( x_out[ q[n_q-1]+2:* ] - x_out[ q[n_q-1]+1: n_elements(x_out)-2 ] )
wdii = where( ii gt q[n_q-1] )
wd = imap[ wdii ]
result = thm_get_efi_edc_offset( d.y[wd, j], samp_rate, nspins, spin_period, edge_truncate, window_truncated, new_nspins )
if result[0] ne -1 then begin
d_smooth[wd, j] = temporary(result)
if window_truncated then begin
offset_estimation_window_truncated = 1b
message, /info, $
' *** WARNING: smoothing window truncated to '+string(new_nspins,format='(i0)')+ ' spins to fit data sub-interval.'
endif
endif else message, /info, ' *** WARNING: Interval less than 1 spin wide. EDC offset not calculated and not subtracted!'
endif else begin
;
; No gaps detected: No need to loop on chunks -- smooth entire period.
;
samp_rate = 1./median( d.x[ 1:* ] - d.x[ 0: n_elements( d.x )-2 ] )
result = thm_get_efi_edc_offset( d.y[*, j], samp_rate, nspins, spin_period, edge_truncate, window_truncated, new_nspins )
if result[0] ne -1 then begin
d_smooth[*, j] = temporary(result)
if window_truncated then begin
offset_estimation_window_truncated = 1b
message, /info, $
' *** WARNING: smoothing window truncated to '+string(new_nspins,format='(i0)')+ ' spins to fit data sub-interval.'
endif
endif else message, /info, ' *** WARNING: Interval less than 1 spin wide. EDC offset not calculated and not subtracted!'
endelse
endfor ;j
; ; Notes on hanning method (would have to replace SMOOTH everywhere -- and also would need to adjust warning messages):
; ;
; kernel = hanning(npts)
; for icomp = 0, 1l do d_smooth[*,icomp] = $
; fix( round( convol( float(d.y[w, icomp]), kernel, total(kernel), /nan ) ) ) ;convol in float, return in int.
;GAIN is time-independent case:
;==============================
if (size(gain,/dimensions))[0] le 1 then begin
; for icomp=0L,2L do res[ w, icomp] = -1000.*gain[icomp] * (d.y[w,icomp] - cp.offset[i,icomp])/exx[icomp] ;EDC offset from calib. file.
for icomp=0L,1L do res[ w, icomp] = -1000.*gain[icomp] * (d.y[w,icomp] - d_smooth[*,icomp] )/exx[icomp] ;EDC removal on-the-fly.
for icomp=2L,2L do res[ w, icomp] = -1000.*gain[icomp] * (d.y[w,icomp] - cp.offset[i,icomp])/exx[icomp] ;EDC offset from calib. file.
;
;GAIN is an array case:
;======================
endif else begin
for icomp=0L,2L do res[ w, icomp] = -1000.*gain[*,icomp] * (d.y[w,icomp] - cp.offset[i,icomp])/exx[icomp] ;milivolts/meter [mV/m]
endelse
end
endcase
endfor
;; update the DLIMIT elements to reflect RAW->PHYS
;; transformation, coordinate system, etc.
;; ***********************************************
;
units = cp.units
;
str_element, dl, 'data_att', data_att, success=has_data_att
if has_data_att then begin
str_element, data_att, 'data_type', 'calibrated', /add
endif else data_att = { data_type: 'calibrated' }
str_element, data_att, 'coord_sys', 'spg', /add
str_element, data_att, 'units', units[0], /add ;We cheat here. Only units of first interval are shown!
str_element, data_att, 'cal_par_time', cp.cal_par_time, /add ;Different w/ T-D calibration.
;"Documented calibration factors" to match expanded
;cal. file params. (3/12/2008):
;*** All of these need a time dimensions w/ T-D calibration ***
;==============================================================
str_element,data_att,'offset',cp.offset,/add
str_element,data_att,'edc_gain',cp.edc_gain,/add
str_element,data_att,'eac_gain',cp.eac_gain,/add
str_element,data_att,'boom_length',cp.boom_length,/add
str_element,data_att,'boom_shorting_factor', cp.boom_shorting_factor,/add
str_element,data_att,'dsc_offset',cp.dsc_offset,/add
; Store metadata for the on-the-fly EDC offset calculation:
;
str_element,data_att,'number_of_spins',nspins,/add
str_element,data_att,'gap_trigger_value',gap_trigger_value,/add
str_element,data_att,'edge_truncate',edge_truncate,/add
str_element,data_att,'offset_estimation_window_truncated',offset_estimation_window_truncated,/add
str_element,data_att,'offset_estimation_func',offset_estimation_func ,/add
str_element, dl, 'data_att', data_att, /add
str_element, dl, 'ytitle', string(tplot_var, units[0], format= $ ;We cheat here. Only units of first interval are shown!
'(A,"!C!C[",A,"]")'), /add
str_element, dl, 'labels', [ 'E12', 'E34', 'E56'], /add
str_element, dl, 'labflag', 1, /add
str_element, dl, 'colors', [ 2, 4,6]
;Store the transformed spectra into the output
;TPLOT variable:
;===============
store_data, tplot_var, data = { x:d.x, y: temporary(res), v:d.v }, lim=l, dlim=dl
if ~(~size(stored_tnames,/type)) then stored_tnames = [ stored_tnames, tplot_var ] else stored_tnames = [ tplot_var ]
;If these are dot0 or 0 datatypes, then transform to dsl
;for further processing, else transform according to coord
;kw (if not 'spg'), but transform to dsl first!
;(labels: spun coord. systems get "E12,34,56"; despun systems get "Ex,y,z"):
;===========================================================================
if where(name eq [udot0s, u0s]) ge 0 then begin
thm_efi_despin, sc, nameraw, cp.dsc_offset[0,*], [1,1,1], $ ;Pass [1,1,1] in order to not duplicate boom_shorting factor.
tplot_name=tplot_var, newname=tplot_var, stored_tnames = foo
if ~(~size(stored_tnames,/type)) then stored_tnames = [ stored_tnames, foo ] else stored_tnames = foo
thm_update_efi_labels,tplot_var
endif else if keyword_set(coord) && coord ne 'spg' then begin
thm_efi_despin, sc, nameraw, cp.dsc_offset[0,*], [1,1,1], $ ;Pass [1,1,1] in order to not duplicate boom_shorting factor.
tplot_name=tplot_var, newname=tplot_var, stored_tnames = foo
if ~(~size(stored_tnames,/type)) then stored_tnames = [ stored_tnames, foo ] else stored_tnames = foo
if coord ne 'dsl' then thm_cotrans,tplot_var,out_coord=coord
endif
thm_update_efi_labels,tplot_var
break
end ; of { EFF, EFP, EFW} calibration clause.
else: begin ; improperly defined name of EFI quantity.
message,/continue, $
string( tplot_var_raw, tplot_var_hed, format= $
'("Improperly defined name of EFI quantity ("'+ $
',A,X,A,").")')
end
endswitch
endif else begin ; necessary TPLOT variables not present.
if keyword_set(vb) then message, /continue, $
string( tplot_var_raw, tplot_var_hed, format= '("necessary TPLOT variables (",A,X,A,") ' + 'not present for RAW->PHYS transformation.")')
endelse
endelse ; done with RAW->PHYS transformation
;stop
;=======================================
;Finish processing of _0 quantity, if present
;and transform coordinates as requested:
;=======================================
if where(name eq u0s) ne -1 then begin
get_data, tplot_var, data=d, limit=l, dlim=dl
if size(d, /type) ne 8 then continue
res = d.y
res[*,2]=0.0
;
store_data, tplot_var, data = { x:d.x, y: temporary(res), v:d.v }, lim=l, dlim=dl
if ~(~size(stored_tnames,/type)) then stored_tnames = [ stored_tnames, tplot_var ] else stored_tnames = [ tplot_var ]
if keyword_set(coord) && coord ne 'dsl' then begin
if coord ne 'spg' and coord ne 'ssl' then begin
thm_update_efi_labels,tplot_var
message, /continue, $
'Warning: coord keyword ignored: ' + $
'only ssl, dsl, or spg are physically meaningful ' + $
'coordinate systems for _0 quantites.'
endif else begin
thm_cotrans,tplot_var,out_coord=coord
thm_update_efi_labels,tplot_var
endelse
endif
endif
;================================================
;Finish processing of _dot0 quantity, if present:
;================================================
if where(name eq udot0s) ne -1 then begin
get_data, tplot_var, data=d, limit=l, dlim=dl
if size(d, /type) ne 8 then continue
; res = d.y
;==========================================
;Temporarily load FGM data as a TPLOT var.:
;==========================================
suff = '_thm_cal_efi_priv'
thm_load_fgm,probe=sc,level='l2',coord='dsl', suffix = suff
if keyword_set(fgm_datatype) then ftype = fgm_datatype else begin
if nameraw eq 'efp' || nameraw eq 'efw' $
then ftype='fgh' $
else ftype='fgl'
endelse
if vb ge 2 then print, 'Using '+ftype
;===========
;Degap data:
;===========
; tdegap, 'th'+sc+'_'+ftype+'_dsl'+suff, dt = cp.cal_par_time, margin = cp.cal_par_time * 0.2, /overwrite ;These inputs don't make sense.
tdegap, 'th'+sc+'_'+ftype+'_dsl'+suff, /overwrite, /onenanpergap, /nowarning ;Try defaults.
;=======================================================
;Get relevant data and delete temporary TPLOT variables:
;=======================================================
get_data,'th'+sc+'_'+ftype+'_dsl'+suff, data=thx_fgx_dsl, limit=fl, dlim=fdl
del_data, '*'+suff
;==================================
;Interpolate FGM data to EFI times:
;==================================
Bx = interpol(thx_fgx_dsl.y[*,0],thx_fgx_dsl.x, d.x)
By = interpol(thx_fgx_dsl.y[*,1],thx_fgx_dsl.x, d.x)
Bz = interpol(thx_fgx_dsl.y[*,2],thx_fgx_dsl.x, d.x) + bz_offset
;
undefine, thx_fgx_dsl, fl, fdl
;=====================
;Initialize Ez to NaN:
;=====================
; res[*,2] = !values.f_nan ;Assign.
d.y[*,2] = !values.f_nan ;Assign.
;==================================================
;Apply max_angle kw, if present. If not, apply
;min_bz kw, if present. If not, let all B's be good.
;If there are no good elements, then continue loop:
;==================================================
if keyword_set(max_angle) then begin
angle=acos(Bz/sqrt(Bx^2+By^2+Bz^2))*180d/!dpi
good = where(abs(angle) le max_angle, ngood)
if ngood eq 0 then continue
endif else if keyword_set(min_bz) then begin
good = where(abs(Bz) ge min_bz, ngood)
if ngood eq 0 then continue
endif else good = lindgen(n_elements(d.x))
;Calculate Bx over Bz, By over Bz:
;=================================
bx_bz = Bx[good]/Bz[good]
by_bz = By[good]/Bz[good]
;===================================
;Apply max_bxy_bz kw, if present. If not, keep all.
;If no good elements, continue loop:
;===================================
if keyword_set(max_bxy_bz) then begin
keep = where(abs(bx_bz) lt max_bxy_bz and abs(by_bz) lt max_bxy_bz, nkeep)
if nkeep eq 0 then continue
endif else keep = lindgen(ngood)
goodkeep=good[keep] ;Otherwise this gets eval'td 3X.
;Calculate Ez and assign directly to result:
;===========================================
; res[goodkeep,2] = -(bx_bz[keep]*res[goodkeep,0] + $
; by_bz[keep]*res[goodkeep,1])
d.y[goodkeep,2] = -(bx_bz[keep]*d.y[goodkeep,0] + $
by_bz[keep]*d.y[goodkeep,1])
;Store result:
;=============
; store_data, tplot_var, data = { x:d.x, y: temporary(res), v:d.v }, lim=l, dlim=dl
store_data, tplot_var, data = { x:d.x, y: d.y, v:d.v }, lim=l, dlim=dl
if ~(~size(stored_tnames,/type)) then stored_tnames = [ stored_tnames, tplot_var ] else stored_tnames = [ tplot_var ]
undefine, d
;At this point, TPLOT_VAR is in dsl. If COORD ne 'dsl', then transform via THM_COTRANS.PRO
;and set lables (spun coord. systems get "E12,34,56"; despun systems get "Ex,y,z"):
;==================================================================================
if keyword_set(coord) && coord ne 'dsl' then begin
thm_cotrans,tplot_var,out_coord=coord
thm_update_efi_labels, tplot_var
endif
endif
endfor ; loop over names.
endfor ; loop over spacecraft.
end