;+
;Purpose:
; Helper function to check support data's existence and time range
;
;Input:
; tvar - string representing tplot variable
; tr - two element numerical time range
;
;Output:
; Returns 1 if time range is covered by the data
;-
function thm_part_getspec_tcheck, tvar, tr
compile_opt idl2, hidden
var = tnames(tvar)
if ~is_string(var) then return, 0
get_data, var, t
;this should be quick compared with total time req
dt = median(t - shift(t,1))
;pads data's time range with dist between last two points
drange = minmax(t) + dt * [-1,1]
return, tr[0] ge drange[0] and tr[1] le drange[1]
end
;+
;PROCEDURE: thm_part_getspec
;PURPOSE:
; Generates tplot variables containing energy-time angular particle spectra
; for any angle range and/or any energy range.
;
;KEYWORDS:
; probe = Probe name. The default is 'all', i.e., load all available probes.
; This can be an array of strings, e.g., ['a', 'b'] or a
; single string delimited by spaces, e.g., 'a b'
; trange = (Optional) Time range of interest (2 element array), if
; this is not set, the default is to prompt the user. Note
; that if the input time range is not a full day, a full
; day's data is loaded.
; phi = Angle range of interest (2 element array) in degrees relative to
; probe-sun direction in the probe's spin plane. Specify angles in
; ascending order (e.g. [90, 180]) to specify the 'daylight'
; hemisphere in DSL coordinates. Default is all 360 degrees (e.g.
; [0, 360]). If the phi range is greater than 360 degrees, then the phi
; bins at the beginning of the phi range are added and wrapped around
; the end of phi range. For example, if phi=[-90,360], the phi bins
; corresponding to 270-360 degrees are appended to the bottom of the plot.
; Also controls which bins are selected for
; pitch/gyrovelocity spectra. Phi is limited to be between -360
; and 720 degrees. Note that, if the ENERGY keyword is set (as
; is true in the default case). This is also true for 'pa'
; and 'gyro' cases.
; theta = Angle range of interest (2 element array) in degrees relative to
; spin plane, e.g. [-90, 0] or [-45, 45] in the probe's spin plane.
; Specify in acending order. Default is all (e.g. [-90, 90]). Also
; controls which bins are selected for pitch/gyrovelocity spectra.
; pitch = Angle range of interest (2 element array) in degrees relative to
; the magnetic field. Default is all (e.g. [0, 180]). Has no effect
; on phi/theta spectra plots.
; gyro = Gyrovelocity angle range of interest (2 element array) in degrees.
; Same rules as specifying phi apply to gyrophse. Default is all
; (e.g. [0, 360]). Has no effect on phi/theta spectra plots.
; erange = Energy range (in eV) of interest (2 element array). Default is all.
; data_type = The type of data to be loaded (e.g. ['peif', 'peef'])
; start_angle = Angle at which to begin plotting the data the bottom of the
; plot along the y-axis. At this time this works for phi only. If
; not set, keyword defaults to first element of phi input.
; suffix = Suffix to add to tplot variable names.
; /AUTOPLOT: Set to activiate a simple plotting routine to display tplot
; variables created by this routine.
; /ENERGY: Set to output tplot variables containing the energy spectrum of the
; energy/angle ranges input above.
; angle = Type of angular spectrum tplot variable created and plotted using the
; energy/angle ranges input above e.g. 'phi'. Default is 'phi'. The
; type of angular spectra is appended to its corresponding tplot
; variable name. All possible options are:
; phi = phi angular spectrum
; theta = theta angular spectrum
; pa = pitch angular spectrum (only works with full angle ranges)
; gyro = gyrovelocity angular spectrum (only works with full angle
; ranges)
; NOTE: If neither of ENERGY and ANGLE keywords are specified, then
; then both are turned on with ANGLE defaulting to 'phi'.
; /GET_SUPPORT_DATA: load support_data variables as well as data variables into
; tplot variables.
; regrid = [m,n] Number of angle bins, used to re-grid (2-element array) along
; phi (gyrovelocity) and theta (pitch), respectively, when calculating
; pitch/gyrovelocity angle spectrum. For example, e.g. [32,16] creates
; a regularly-spaced array of 32 phi angles by 16 theta angles that
; resample the probe(s) native angle bin array. Default is [16,8].
; Keyword has no effect on phi and theta spectragrams. Accuracy of the
; spectragrams increases when more angle bins are used to re-grid.
; Suitable numbers for m and n are numbers like 2^k. So far, k=2-6 has
; been tested and work. Other numbers will work provided 180/n and
; 360/m are rational.
; other_dim = Keyword passed to THM_FAC_MATRIX_MAKE for conversion to field
; aligned coordinates to create pitch angle and gyrovelocity
; spectra. See THM_FAC_MATRIX_MAKE for valid input. Default is
; 'mphigeo'.
; /NORMALIZE: Set to normalize the flux for each time sample to 0-1.
; badbins2mask = A 0-1 array that indicates which SST bins will be masked with
; NaN to eliminate things like sun contamination. The array
; should have the same number of elements as the number of angle
; bins for a given data type. A 0 indicates that will be masked
; with a NaN. This is basically the output from the bins argument
; of EDIT3DBINS.
; datagap = Maximum time gap in seconds over which to interpolate the plot.
; Sets the DATAGAP flag in the dlimits which SPECPLOT uses to
; interpolate pixels in time gaps less than DATAGAP. Use this keyword
; when overlaying spectra plots, allowing the underlying spectra to be
; shown in the data gaps of the overlying spectra. Default is zero.
;
;ESA PEER/PEIR/PEIF Background Removal Keywords:
;
;/bdnd_remove: Turn on ESA background removal.
;
;bgnd_type(Default 'anode'): Set to string naming background removal type:
;'angle','omni', or 'anode'.
;
;bgnd_npoints(Default = 3): Set to the number of lowest values points to average over when determining background.
;
;bgnd_scale(Default=1): Set to a scaling factor that the background will be multiplied by before it is subtracted
;
;GUI-RELATED KEYWORDS:
; /GUI_FLAG: Flag tells code to recognize status bar and history window objects.
; gui_statusBar = Object reference to status bar object in GUI.
; gui_historyWin = Object reference to history window object in GUI.
;
;EXAMPLE: thm_part_getspec, probe='d', trange=['07-06-17','07-06-19']
;
;SEE ALSO:
; THM_CRIB_PART_GETSPEC, THM_PART_MOMENTS2, THM_PART_GETANBINS, THM_LOAD_SST,
; THM_LOAD_ESA_PKT, THM_FAC_MATRIX_MAKE,THM_REMOVE_SUNPULSE,THM_CRIB_SST_CONTAMINATION
;
;NOTES: For documentation on sun contamination correction keywords that
; may be passed in through the _extra keyword please see:
; thm_sst_remove_sunpulse.pro or thm_crib_sst_contamination.pro
;
;BACKGROUND REMOVAL(BGND) Description, Warnings and Caveats(from Vassilis Angelopoulos):
; This code allows for keywords that permit omni-directional or anode-dependent
; background removal from penetrating electrons in the ESA ion and electron
; detectors. Anode-dependent subtraction is used when possible by default,
; i.e., when angle information is available; but user has full control by
; keyword specification. Default bgnd estimates use 3 lowest counts/s values.
; Scaling of the background (artificial scaling) can also allow playing with
; background estimates to account for noise statistics in the background itself.
; The parameters that have worked well for me during high bgnd levels are:
; ,/bgnd_remove, bgnd_type='anode', bgnd_npoints=3, bgnd_scale=1.5
;
; The same keywords when used in thm_part_getspec, and thm_part_moments
; are understood and passed to the data extraction routines, such that
; they will do the removal before computing moments or spectra.
;
; This background subtraction to be used at the inner magnetosphere,
; or when SST electron fluxes indicate presence of significant electron
; fluxes at the satellite (injections). At quiet times the code tends to remove
; real fluxes, so beware.
;
;
;CREATED BY: Bryan Kerr
;
; $LastChangedBy: bckerr $
; $LastChangedDate: 2008-06-13 13:29:12 -0700 (Fri, 13 Jun 2008) $
; $LastChangedRevision: 3204 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/thmsoc/trunk/idl/themis/spacecraft/particles/thm_part_getspec.pro $
;-
pro thm_part_getspec, probe=probe, trange=trange, phi=phi, theta=theta, $
pitch=pitch, erange=erange, data_type=data_type, $
start_angle=start_angle, suffix=suffix, $
autoplot=autoplot, angle=angle, energy=energy, $
nomask=nomask, get_support_data=get_support_data, $
regrid=regrid, gyro=gyro, other_dim=other_dim, $
en_tnames=en_tnames, an_tnames=an_tnames, $
normalize=normalize, bins2mask=bins2mask, $
badbins2mask=badbins2mask, datagap=datagap, $
forceload=forceload, $
; gui-related keywords
gui_flag=gui_flag, gui_statusBar=gui_statusBar, $
gui_historyWin=gui_historyWin, $
sst_cal=sst_cal,$
; misc keywords
_extra=ex, test=test
if n_elements(gui_flag) eq 0 then gui_flag = 0b ;gui_flag for messages
;******************************************************************************
; This section warns when obsolete keywords are used.
;******************************************************************************
if keyword_set(nomask) then begin
print,''
dprint, 'The /NOMASK keyword has been replaced by the /MASK keyword. ', $
'Please see documentation (thm_crib_part_getspec) for usage info.'
print,''
return
endif
if keyword_set(bins2mask) then begin
print,''
dprint, 'The BINS2MASK keyword has been replaced by the BADBINS2MASK ', $
'keyword. Please see documentation (thm_crib_part_getspec) for usage info.'
print,''
return
endif
;******************************************************************************
; Set default input
;******************************************************************************
deftheta = [-90, 90]
defphi = [0, 360]
defpitch = [0, 180]
defgyro = defphi
deferange = [0, 1e7]
defdata_type = 'p??f'
def_angle = 'phi'
defstart_angle = 0
wrapphi = 0
def_other_dim='mphigeo'
vdatatypes = ['peif','peef','psif','psef','peir','peer','psir','pser','peib','peeb','pseb']
datatype = strfilter(vdatatypes,data_type,/fold_case,delimiter=' ',count=ndatatypes)
if ndatatypes eq 0 then begin
dprint,dlevel=1,"No Valid Data Type Selected"
return
endif
;******************************************************************************
; Check for correct angle input
;******************************************************************************
if keyword_set(theta) then begin
if theta[1] lt theta[0] OR (theta[0] lt -90 OR theta[1] gt 90) then begin
theta_mess = 'Error with theta input. Does not satisfy: -90 < theta[0] < theta[1] < 90.'
dprint, ''
dprint, theta_mess
if gui_flag then begin
gui_statusBar -> update, theta_mess
gui_historyWin -> update, theta_mess
endif
return
endif
endif else theta=deftheta
if keyword_set(phi) then begin
phi_mess = ''
if(min(phi) Lt -360.0 Or max(phi) Gt 720.0) then begin
phi_mess = 'Error with phi input. Phi angles must be between -360.0 and 720.0.'
endif
if phi[0] gt phi[1] then begin
phi_mess = 'Error with phi input. Phi angles must be specified in ascending order.'
endif
If(is_string(phi_mess)) Then Begin
dprint, ''
dprint, phi_mess
if gui_flag then begin
gui_statusBar -> update, phi_mess
gui_historyWin -> update, phi_mess
endif
Return
Endif
if (phi[1] - phi[0]) gt 360 then begin
;if angle eq 'phi' || angle eq 'theta' then wrapphi = phi[1] - phi[0] - 360
wrapphi = phi[1] - phi[0] - 360
phi[1] = phi[0] + 360
endif
endif else phi=defphi
if keyword_set(pitch) then begin
if pitch[1] lt pitch[0] OR (pitch[0] lt 0 OR pitch[1] gt 180) then begin
pa_mess = 'Error with pitch input. Does not satisfy: 0 < pitch[0] < pitch[1] < 180.'
dprint, ''
dprint, pa_mess
if gui_flag then begin
gui_statusBar -> update, pa_mess
gui_historyWin -> update, pa_mess
endif
return
endif
endif else pitch=defpitch
if keyword_set(gyro) then begin
if gyro[0] gt gyro[1] then begin
gyro_mess = 'Error with gyro input. Gyrovelocity angles must be specified in ascending order.'
dprint, ''
dprint, gyro_mess
if gui_flag then begin
gui_statusBar -> update, gyro_mess
gui_historyWin -> update, gyro_mess
endif
return
endif
if (gyro[1] - gyro[0]) gt 360 then begin
wrapphi = gyro[1] - gyro[0] - 360
gyro[1] = gyro[0] + 360
endif
endif else begin
gyro=defgyro
if keyword_set(angle) then begin
if angle eq 'gyro' || angle eq 'pa' then wrapphi = 0
endif
endelse
;******************************************************************************
; Check for type of spectra requested (e.g. angle or energy)
;******************************************************************************
if ~ keyword_set(angle) AND ~ keyword_set(energy) then begin
angle = def_angle
dprint, 'setting angle data_type: ', angle
energy=1
endif else if ~ keyword_set(angle) then angle = 0
if keyword_set(erange) eq 0 then erange=deferange
if keyword_set(data_type) eq 0 then begin
data_type = defdata_type
dprint, ''
dprint, 'setting default data_type: ', data_type
dprint, ''
endif
theta1 = min(theta)
theta2 = max(theta)
phi1 = phi[0]
phi2 = phi[1]
if keyword_set(trange) then begin
tr = trange ; copy trange to internal variable
trd = time_double(tr)
ndays = (trd[1] - trd[0]) / 86400
timespan,trd[0],ndays
endif else begin
tr = time_string(timerange())
trd = time_double(tr)
endelse
if keyword_set(badbins2mask) then begin
if array_equal(badbins2mask, -1) then begin
dprint,''
dprint,'WARNING: BADBINS2MASK array is empty. Not masking any SST angle bins.'
dprint,''
endif else begin
dprint,'Masking the following SST angle bins: ', fix(where(badbins2mask eq 0))
endelse
endif
;******************************************************************************
; check for sst data_type requests and load sst data
;******************************************************************************
sst_ind = where(strmid(datatype,1,1) eq 's', sst_count)
if sst_count ne 0 then begin
sst_type = datatype[sst_ind]
if keyword_set(sst_cal) then begin
thm_load_sst2, probe=probe, trange=trd, datatype=sst_type
endif else begin
if ~keyword_set(forceload) && thm_check_part_data(probe, sst_type, trd) then begin
dprint, 'Using previously loaded data...'
endif else begin
thm_load_sst, probe=probe, trange=trd, datatype=sst_type, $
get_support_data=get_support_data
endelse
endelse
endif
;******************************************************************************
; check for esa data_type requests and load esa data
;******************************************************************************
esa_ind = where(strmid(datatype,1,1) eq 'e', esa_count)
if esa_count ne 0 then begin
esa_type = datatype[esa_ind]
if ~keyword_set(forceload) && thm_check_part_data(probe, esa_type, trd) then begin
dprint, 'Using previously loaded data...'
endif else begin
thm_load_esa_pkt, probe=probe, trange=trd, datatype=esa_type, $
get_support_data=get_support_data
endelse
endif
;******************************************************************************
; Setup for pitch angle/gyrovelocity spectra. Load state and mag data.
;******************************************************************************
if strcmp(angle,'pa') || strcmp(angle,'gyro') || keyword_set(energy) then begin
;If the ENERGY keyword is set, this code is needed so that eflux spectra can
;be subject to pitch/gyrovelocity angle restrictions
if ~ keyword_set(regrid) then begin
dprint, 'REGRID keyword not set. Defaulting to [16 phis x 8 thetas].'
regrid=[16,8]
endif
if ~ keyword_set(other_dim) then other_dim=def_other_dim
; move phi minimum to positive value
If(phi1 Lt 0) Then Begin
phi1 = phi1+360
phi2 = phi2+360
Endif
phi = [phi1, phi2]
;The following code block sets phi range to 0, 360 for cases with
;any values outsideof 0, 360 (e.g., -90, 90, which should be set to
;270, 450, is set to 0, 360).
; if ~ ((phi1 le 360 && phi1 ge 0) && (phi2 le 360 && phi2 ge 0)) then begin
; if phi2 lt (phi2 - phi1) then begin
; phi1_1 = 360 + phi1
; phi1_2 = 360
; phi2_1 = 0
; phi2_2 = phi2
; endif
; if phi2 gt (phi2 - phi1) then begin
; phi1_1 = phi1
; phi1_2 = 360
; phi2_1 = 0
; phi2_2 = phi2 - 360
; endif
; endif else begin
; phi1_1 = phi1
; phi1_2 = phi2
; phi2_1 = phi1
; phi2_2 = phi2
; endelse
; phi[0] = min([phi1_1, phi1_2, phi2_1, phi2_2])
; phi[1] = max([phi1_1, phi1_2, phi2_1, phi2_2])
coord = 'dsl'
mag_suffix = '_fgs' ;+coord;+'_'+angle
;Check first for state and FIT/FGM data for each probe, before loading
;similar to thm_ui_check4spin, but without the common block, jmm, 1-9-2009
p1 = thm_valid_input(probe, vinputs = ['a b c d e'], definput = 'a', /include_all)
tr0 = time_double(tr)
For j = 0, n_elements(p1)-1 Do Begin
;Get state support data
var1 = 'th'+p1[j]+'_state_spinper'
var2 = 'th'+p1[j]+'_state_spinphase'
;load data if not already present
if ~thm_part_getspec_tcheck(var1, tr0) || $
~thm_part_getspec_tcheck(var2, tr0) then begin
thm_load_state, probe = p1[j], trange = trd, /get_support_data
if ~thm_part_getspec_tcheck(var1, tr0) || $
~thm_part_getspec_tcheck(var2, tr0) then begin
dprint, 'WARNING: STATE support data does not cover the requested time range.'
endif
endif
;Get lvl1 mag data
;
var = 'th'+p1[j]+'_fgs'
;load lvl 1 FIT data if variable does not already exist
if ~thm_part_getspec_tcheck(var, tr0) then begin
thm_load_fit, probe=p1[j], trange=trd, datatype='fgs', level='l1', $
coord='dsl', get_support_data=get_support_data
if ~thm_part_getspec_tcheck(var, tr0) then begin
dprint, 'WARNING: B-field data (FIT FGS) does not cover the requested time range.'
endif
endif
Endfor
endif
; Calc and create spectra tplot variables
thm_part_moments2, probe=probe, instruments_types=datatype, trange=trd, $
tplotnames=tn, theta=theta, phi=phi, pitch=pitch, $
erange=erange, tplotsuffix=suffix, start_angle=start_angle, $
doangle=angle, doenergy=energy, wrapphi=wrapphi, $
mag_suffix=mag_suffix, regrid=regrid, $
gyro=gyro, other_dim=other_dim, en_tnames=en_tnames, $
an_tnames=an_tnames, normalize=normalize, $
badbins2mask=badbins2mask, datagap=datagap, _extra=ex, $
test=test, gui_flag=gui_flag, gui_statusBar=gui_statusBar, $
gui_historyWin=gui_historyWin,sst_cal=sst_cal
if ~keyword_set(an_tnames) and ~keyword_set(en_tnames) then begin
dprint, 'WARNING: No tplot names returned by thm_part_moments2.'
endif
;******************************************************************************
; Setup angle/energy spectrum tplot variables calculated by thm_part_moments2
;******************************************************************************
if keyword_set(angle) then begin
stheta = strcompress(string(theta))
sphi = strcompress(string(phi))
serang = strcompress(string(erange))
spitch = strcompress(string(pitch))
sgyro = strcompress(string(gyro))
title = strjoin(['theta=', stheta, ', phi=', sphi, ', erange=', serang])
if angle eq 'gyro' || angle eq 'pa' then begin
title = strjoin(['theta=', stheta, ', phi=', sphi, ', pitch=', spitch, $
', gyro=', sgyro, ', erange=', serang])
endif
tplot_options,'title',''
tplot_options,'xmargin',[15,10]
options,an_tnames,'y_no_interp',1,/default
options,an_tnames,'x_no_interp',1,/default
options,an_tnames,'ztitle','Eflux [eV/cm2/s/ster/eV]',/default
options,an_tnames,minzlog=1,/default
zlim,an_tnames,1,1,1,/default
; ylim,'th*an_eflux*',p_start_angle, p_end_angle,log=0,/default
;tdegap, 'th*an_eflux', overwrite=1
endif
if keyword_set(energy) then begin
stheta = strcompress(string(theta))
sphi = strcompress(string(phi))
serang = strcompress(string(erange))
title = strjoin(['theta=', stheta,', phi=',sphi,', erange=', $
serang])
tplot_options,'title',''
tplot_options,'xmargin',[15,10]
options,en_tnames,'y_no_interp',1,/default
options,en_tnames,'x_no_interp',1,/default
options,en_tnames,'ztitle','Eflux [eV/cm2/s/ster/eV]',/default
options,en_tnames,minzlog=1,/default
zlim,en_tnames,1,1,1,/default
options,en_tnames,ystyle=1,/default
endif
;******************************************************************************
; Autoplot the angle/energy spectra
;******************************************************************************
if keyword_set(autoplot) && keyword_set(angle) then begin
tplot_options,'title',title
tplot,an_tnames
endif
if keyword_set(autoplot) && keyword_set(energy) then begin
if keyword_set(angle) then begin
tplot_options,'title',title
tplot, an_tnames
; look for & plot nrg spectra vars w/more than 1 nrg channel
for i=0,size(en_tnames,/n_elements)-1 do begin
get_data,en_tnames[i],x,y,v
if size(v,/n_dimensions) le 1 then begin
dprint, en_tnames[i]," has only 1 energy channel. Can't plot ", $
'spectra unless ERANGE is larger.'
continue
endif else begin
tplot, en_tnames[i], /ADD
endelse
endfor
endif else begin
tplot_options,'title',title
; look for & plot nrg spectra vars w/more than 1 nrg channel
for i=0,size(en_tnames,/n_elements)-1 do begin
get_data,en_tnames[i],x,y,v
if size(v,/n_dimensions) le 1 then begin
dprint, en_tnames[i]," has only 1 energy channel. Can't plot ", $
'spectra unless ERANGE is larger.'
continue
endif else begin
tplot, en_tnames[i], /ADD
endelse
endfor
endelse
endif
end