;+ ; PROCEDURE: ; rbsp_rbspice_pad ; ; PURPOSE: ; Calculate pitch angle distributions using data from the ; RBSP Radiation Belt Storm Probes Ion Composition Experiment (RBSPICE) ; ; KEYWORDS: ; ; probe: RBSP spacecraft indicator [Options: 'a' (default), 'b'] ; trange: Time range of interest (2 element array); if not set, the default is to prompt the user ; datatype: desired data type [Options: 'TOFxEH' (default), 'TOFxEnonH'] ; level: data level ['l1','l2','l3' (default),'l3pap'] ; energy: user-defined energy range to include in the calculation in keV [default = [0,1000]] ; bin_size: desired size of the pitch angle bins in degrees [default = 15] ; scopes: string array of telescopes to be included in PAD [0-5, default is all] ; ; EXAMPLES: ; ; ; OUTPUT: ; ; ; NOTES: ; This was written by Brian Walsh; minor modifications by egrimes@igpp and Ian Cohen (APL) ; ; REVISION HISTORY: ; + 2017-02-22, I. Cohen : created based on mms_eis_pad.pro, edited to handle omni-directional or single telescope ; ;$LastChangedBy: egrimes $ ;$LastChangedDate: 2017-03-03 08:08:58 -0800 (Fri, 03 Mar 2017) $ ;$LastChangedRevision: 22902 $ ;$URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/spdsoft/tags/spedas_6_0/general/missions/rbsp/rbspice/rbsp_rbspice_pad.pro $ ;- pro rbsp_rbspice_pad,probe = probe, trange = trange, datatype = datatype, level = level, $ energy = energy, bin_size = bin_size, scopes = scopes if not KEYWORD_SET(probe) then probe = 'a' if not KEYWORD_SET(datatype) then datatype = 'TOFxEH' case datatype of 'TOFxEH': species = 'proton' 'TOFxEnonH': species = ['helium','oxygen'] 'TOFxPHHHELT': species = ['proton','oxygen'] endcase if not KEYWORD_SET(level) then level = 'l3' if (level ne 'l1') then units_label = '1/(cm!U2!N-sr-s-keV)' else units_label = 'counts/s' if not KEYWORD_SET(energy) then energy = [0,1000] if not KEYWORD_SET(bin_size) then bin_size = 15. if not KEYWORD_SET(scopes) then scopes = [0,1,2,3,4,5] prefix = 'rbsp'+probe+'_rbspice_'+level+'_'+datatype+'_' if energy[0] gt energy[1] then begin print, 'Low energy must be given first, then high energy in "energy" keyword' return endif ; set up the number of pa bins to create bin_size = float(bin_size) n_pabins = 180./bin_size pa_bins = 180.*indgen(n_pabins+1)/n_pabins pa_label = 180.*indgen(n_pabins)/n_pabins+bin_size/2. dprint, dlevel=0, 'Num PA bins: ', string(n_pabins) dprint, dlevel=0, 'PA bins: ', string(pa_bins) status = 1 ; check to make sure the data exist get_data, prefix + 'Alpha', data=d, index = index if index eq 0 then begin print, 'No '+datatype+' data is currently loaded for probe rbsp-'+probe+' for the selected time period' status = 0 return endif var_data = tnames(prefix+species) ; if data exists continue if status ne 0 then begin for ion_type_idx = 0, n_elements(species)-1 do begin ; get pitch angle data (all telescopes in single variable) get_data, prefix + 'Alpha', data = d_pa pa_file = fltarr(n_elements(d_pa.x),n_elements(scopes)) ; time steps, look direction for aa=0,n_elements(scopes)-1 do pa_file[*,scopes[aa]] = d_pa.y[*,scopes[aa]] pa_flux = fltarr(n_elements(d_pa.x),n_pabins,n_elements(scopes)) pa_flux[where(pa_flux eq 0)] = !values.f_nan pa_num_in_bin = fltarr(n_elements(d_pa.x), n_pabins, n_elements(scopes)) for qq=0,n_elements(species)-1 do begin ; get flux data (all telescopes in single variable) get_data, prefix + species[qq], data = d_flux dprint, dlevel=1, prefix + species[qq] flux_file = fltarr(n_elements(d_flux.x),n_elements(scopes)) ; time steps, look direction flux_file[where(flux_file eq 0)] = !values.f_nan new_pa_flux = fltarr(n_elements(d_flux.x),n_pabins,n_elements(scopes)) ; the average for each bin ; get energy range of interest e = d_flux.v indx = where((e lt energy[1]) and (e gt energy[0]), energy_count) if energy_count eq 0 then begin print, 'Energy range selected is not covered by the detector for ' + datatype + ' ' + species[ion_type_idx] continue endif for t=0,n_elements(scopes)-1 do begin ; Loop through each time step and get: ; 1. the total flux for the energy range of interest for each detector ; 2. flux in each pa bin for i=0, n_elements(d_flux.x)-1 do begin ; loop through time flux_file[i,t] = total(reform(d_flux.y[i,indx,scopes[t]]), /nan) ; start with lowest energy for j=0, n_pabins-1 do begin ; loop through pa bins if (pa_file[i,t] gt pa_bins[j]) and (pa_file[i,t] lt pa_bins[j+1]) then begin if ~finite(pa_flux[i,j,t]) then begin pa_flux[i,j,t] = flux_file[i,t] endif else begin pa_flux[i,j,t] = pa_flux[i,j,t] + flux_file[i,t] endelse pa_num_in_bin[i,j,t] += 1.0 endif endfor endfor ; endfor ; loop over time for i=0, n_elements(pa_flux[*,0,0])-1 do begin ; loop over bins for bin_idx = 0, n_elements(pa_flux[i,*,0])-1 do begin if pa_num_in_bin[i,bin_idx,t] ne 0.0 then begin new_pa_flux[i,bin_idx,t] = pa_flux[i,bin_idx,t]/pa_num_in_bin[i,bin_idx,t] endif else begin new_pa_flux[i,bin_idx,t] = !values.f_nan endelse endfor endfor endfor en_range_string = strcompress(string(energy[0]), /rem) + '-' + strcompress(string(energy[1]), /remove_all) + 'keV' if (n_elements(scopes) eq 6) then begin new_name = prefix+species(qq)+'_omni_'+en_range_string+'_pad' new_omni_pa_flux = reform(new_pa_flux[*,*,0]) for ii=0,n_elements(new_pa_flux[*,0,0])-1 do for jj=0,n_elements(new_pa_flux[0,*,0])-1 do new_omni_pa_flux(ii,jj)= mean(new_pa_flux(ii,jj,*),/NAN) store_data, new_name, data={x:d_flux.x, y:new_omni_pa_flux, v:pa_label} options, new_name, yrange = [0,180], ystyle=1, spec = 1, no_interp=1 , $ ytitle = 'rbsp-'+probe+'!Crbspice!C'+species[ion_type_idx]+'!Comni', ysubtitle=en_range_string+'!CPA [Deg]', ztitle=units_label, minzlog=.01 zlim, new_name, 0, 0, 1 endif else begin new_name = strarr(n_elements(scopes)) for ii=0,n_elements(scopes)-1 do begin new_name = prefix+species(qq)+'_T'+strcompress(string(scopes[ii]),/remove_all)+'_'+en_range_string+'_pad' store_data, new_name, data={x:d_flux.x, y:reform(new_pa_flux[*,*,ii]), v:pa_label} options, new_name, yrange = [0,180], ystyle=1, spec = 1, no_interp=1 , $ ytitle = 'rbsp-'+probe+'!Crbspice!C' +species[ion_type_idx]+'!CT'+strcompress(scopes[t],/remove_all), ysubtitle=en_range_string+'!CPA [Deg]', ztitle=units_label, minzlog=.01 zlim, new_name, 0, 0, 1 endfor endelse ; now do the spin average rbsp_rbspice_pad_spinavg, probe=probe, datatype=datatype, species=species[ion_type_idx], energy=energy, bin_size=bin_size, scopes=scopes endfor endfor endif end