;+ ; ;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 ()" for spinning coordinate systems and ;"Ex ()" for despun coordinate systmes. Y and Z axes are upadated correspondingly. Call after coordinate ;transformations. ; ;SYNTAX: ; thm_update_efi_labels , ; ;Arguments: ; : 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 [, ] ; ;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 . 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