;+
; PROCEDURE:
;         mms_feeps_omni
;
; PURPOSE:
;       Calculates the omni-directional flux for all 24 sensors 
;       
;       (this version re-bins the data due to the different energy channels for each s/c, sensor head and sensor ID)
;       
; INPUT:
;       probe:      spacecraft # (1, 2, 3, or 4)
;
; KEYWORDS:
;
;       datatype:   feeps data types include ['electron', 'electron-bottom', 'electron-top',
;                   'ion', 'ion-bottom', 'ion-top'].
;                   If no value is given the default is 'electron'.
;       data_rate:  instrument data rates for feeps include 'brst' 'srvy'. The
;                   default is 'srvy'
;       tplotnames: names of loaded tplot variables
;       suffix:     suffix used in call to mms_load_data; required to find the correct
;                   variables
;       data_units: specify units for omni-directional calculation
;       quality_flag: Specify the upper level quality index user is comfortable using in calculating omnidirection flux from available eyes. More details in data products guide. 
;                   For BRST or RAW:
;                     Green (Quality Indicator = 0): No contaminated or masked spin sectors (best data quality);
;                     Red (Quality Indicator = 3): Contaminated spin sector is associated with BURST and RAW data
;                     (this data is not recommended for scientific use);
;                     Grey (Quality Indicator = 4): FEEPS calibration data (this data is generated by the instrument
;                     itself and is not recommended for scientific use);
;                   For data in SURVEY mode, there are 5 possible quality indicators:
;                     Green (Quality Indicator = 0): No contaminated or masked spin sectors (best data quality);
;                     Yellow (Quality Indicator = 1): Among the spin sectors used to create a SURVEY aggregate
;                     spin sector, onboard masking is successfully applied to a minimum of 1 and a maximum of 7
;                     spin sectors and there is no contamination in any of the individual spin sectors. Note that
;                     since this indicator applies only when onboard masking is applied, some sectors are not
;                     observed.
;                     Orange (Quality Indicator = 2): There is contamination of between 0 and 50% of the notmasked
;                     sectors. Warning: this data contains contamination.
;                     Red (Quality Indicator = 3): Same as orange except that there is contamination of 50% or
;                     more of the not-masked sectors. IN ADDITION, this indicator applies to the case where
;                     onboard masking is applied to all 8 spin sectors in the aggregate sector.
;                     Grey (Quality Indicator = 4): FEEPS calibration data (this data is generated by the instrument
;                     itself and is not recommended for scientific use);
;                     Added by Christine Gabrielse 10-03-2025
;       get_err:  Gets the error in the omnidirectional FEEPS intensity. 0 value is off, gt 0 sets it on. User MUST run mms_load_feeps for the l1b counts data FIRST to use this function.
;                 Note: the l1b FEEPS data used by this switch is restricted to MMS team members only at present.
;                      
;
; NOTES:
;       New version, 1/26/17 - egrimes
;       Newer version, 1/31/17 - dturner
;       Fixed 2 bugs (3/29/17): 1) off by one bug when setting bottom sensor without data to NaNs, and
;                               2) now initializing output as NaNs, to avoid setting channels with 
;                                 counts=0 to NaN - egrimes
;       New version, 10/03/25 - Christine Gabrielse, added quality flag keyword
;
; $LastChangedBy: jwl $
; $LastChangedDate: 2025-10-08 17:56:45 -0700 (Wed, 08 Oct 2025) $
; $LastChangedRevision: 33719 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/spdsoft/trunk/projects/mms/feeps/mms_feeps_omni.pro $
;-


pro mms_feeps_omni, probe, datatype = datatype, tplotnames = tplotnames, suffix = suffix, $
  data_units = data_units, data_rate = data_rate, level = level, sensor_eyes=sensor_eyes,quality_flag=quality_flag,get_err=get_err
  if undefined(level) then level = 'l2'
  if undefined(probe) then probe = '1' else probe = strcompress(string(probe))
  ; default to electrons
  if undefined(datatype) then datatype = 'electron'
  if undefined(suffix) then suffix = ''
  if undefined(data_rate) then data_rate = 'srvy'
  if undefined(data_units) then data_units = 'cps'
  if (data_units eq 'flux') then data_units = 'intensity'
  if (data_units eq 'cps') then data_units = 'count_rate'
  units_label = data_units eq 'intensity' ? '1/(cm!U2!N-sr-s-keV)' : 'Counts/s'
  if undefined(quality_flag) then quality_flag=3. ;Added by Christine Gabrielse 10-03-2025
  if undefined(get_err) then get_err=0. ;Added by Christine Gabrielse 10-03-2025
  
  probe = strcompress(string(probe), /rem)

  prefix = 'mms'+probe+'_epd_feeps_'

  if datatype eq 'electron' then begin
    energies = [33.200000d, 51.900000d, 70.600000d, 89.400000d, 107.10000d, 125.20000d, 146.50000d, 171.30000d, $
      200.20000d, 234.00000d, 273.40000, 319.40000d, 373.20000d, 436.00000d, 509.20000d]
  endif else energies = [57.900000d, 76.800000d, 95.400000d, 114.10000d, 133.00000d, 153.70000d, 177.60000d, $
       205.10000d, 236.70000d, 273.20000d, 315.40000d, 363.80000d, 419.70000d, 484.20000d,  558.60000d]
  
  ; Added by DLT on 31 Jan 2017: set unique energy and gain correction factors per spacecraft
  eEcorr = [14.0, -1.0, -3.0, -3.0]
  iEcorr = [0.0, 0.0, 0.0, 0.0]
  eGfact = [1.0, 1.0, 1.0, 1.0]
  iGfact = [0.84, 1.0, 1.0, 1.0]
  
  ; Added by DLT on 31 Jan 2017: set unique energy bins per spacecraft     
  ; Electrons:
  if probe eq '1' && datatype eq 'electron' then energies = energies + eEcorr[0]
  if probe eq '2' && datatype eq 'electron' then energies = energies + eEcorr[1]
  if probe eq '3' && datatype eq 'electron' then energies = energies + eEcorr[2]
  if probe eq '4' && datatype eq 'electron' then energies = energies + eEcorr[3]
  ; Ions:
  if probe eq '1' && datatype eq 'ion' then energies = energies + iEcorr[0]
  if probe eq '2' && datatype eq 'ion' then energies = energies + iEcorr[1]
  if probe eq '3' && datatype eq 'ion' then energies = energies + iEcorr[2]
  if probe eq '4' && datatype eq 'ion' then energies = energies + iEcorr[3]
  
 
  ;en_bins = bin_edges(energies) ; commented by DLT on 31 Jan 2017
  n_enbins = n_elements(energies) ; n_elements(en_bins)-1 ; changed by DLT on 31 Jan 2017
  en_label = energies
  en_chk = 0.10  ; percent error around energy bin center to accept data for averaging; anything outside of energies[i] +/- en_chk*energies[i] will be changed to NAN and not averaged 
  
  top_sensors = sensor_eyes['top']
  bot_sensors = sensor_eyes['bottom']
  
  var_name = strcompress(prefix+data_rate+'_'+level+'_'+datatype+'_top_'+data_units+'_sensorid_'+string(top_sensors[0])+'_clean_sun_removed'+suffix, /rem)
  get_data, var_name, data = d, dlimits=dl
  
  if is_struct(d) then begin
    if get_err gt 0. then counts_arr=dblarr(n_elements(d.x), n_elements(d.v), n_elements(top_sensors)+n_elements(bot_sensors))+!values.d_nan
     
    flux_omni = dblarr(n_elements(d.x), n_elements(d.v))
    if level ne 'sitl' then begin
      dalleyes = dblarr(n_elements(d.x), n_elements(d.v), n_elements(top_sensors)+n_elements(bot_sensors))+!values.d_nan
      for isen = 0, n_elements(top_sensors)-1 do begin
        ; Top units:
        var_name = strcompress(prefix+data_rate+'_'+level+'_'+datatype+'_top_'+data_units+'_sensorid_'+string(top_sensors[fix(isen)])+'_clean_sun_removed'+suffix, /rem)
        get_data, var_name, data = d, dlimits=dl
        dalleyes[*,*,isen] = reform(d.y)
        iE = where(abs(energies - d.v) gt en_chk*energies) ; Check for energies beyond en_chk [%] of the corrected energy bin center and replace with NAN
        if iE[0] ne -1 then dalleyes[*,iE,isen] = !values.d_nan
        
        ;;Added 10-03-2025
        if level eq 'l2' then begin
          var_name2 = strcompress(prefix+data_rate+'_'+level+'_'+datatype+'_top_quality_indicator_sensorid_'+string(top_sensors[fix(isen)])+suffix, /rem)   
          get_data, var_name2, data = d2, dlimits=dl2
          bad_quality=where(d2.y ge quality_flag)
          if bad_quality[0] ne -1 then dalleyes[bad_quality,*,isen] = !values.d_nan
          
          if get_err gt 0. then begin
            var_name3 = strcompress(prefix+data_rate+'_l1b_'+datatype+'_top_counts_sensorid_'+string(top_sensors[fix(isen)])+'_clean_sun_removed'+suffix, /rem) 
            get_data, var_name3, data = d3, dlimits=dl3
            counts_arr[*,*,isen] = reform(d3.y)
            if iE[0] ne -1 then counts_arr[*,iE,isen] = !values.d_nan
            if bad_quality[0] ne -1 then counts_arr[bad_quality,*,isen] = !values.d_nan        
          endif
        endif
        ;;;
       ; stop
      endfor
      for isen = 0, n_elements(bot_sensors)-1 do begin
        ; Bottom units:
        var_name = strcompress(prefix+data_rate+'_'+level+'_'+datatype+'_bottom_'+data_units+'_sensorid_'+string(bot_sensors[fix(isen)])+'_clean_sun_removed'+suffix, /rem)
        get_data, var_name, data = d, dlimits=dl
        dalleyes[*,*,isen+n_elements(top_sensors)] = reform(d.y)
        iE = where(abs(energies - d.v) gt en_chk*energies) ; Check for energies beyond en_chk [%] of the corrected energy bin center and replace with NAN
        if iE[0] ne -1 then dalleyes[*,iE,isen+n_elements(top_sensors)] = !values.d_nan
        ;;Added 10-03-2025
        if level eq 'l2' then begin
          var_name2 = strcompress(prefix+data_rate+'_'+level+'_'+datatype+'_bottom_quality_indicator_sensorid_'+string(bot_sensors[fix(isen)])+suffix, /rem)
          get_data, var_name2, data = d2, dlimits=dl2
          bad_quality=where(d2.y ge quality_flag)
          if bad_quality[0] ne -1 then dalleyes[bad_quality,*,isen+n_elements(top_sensors)] = !values.d_nan
          
          if get_err gt 0. then begin
            var_name3 = strcompress(prefix+data_rate+'_l1b_'+datatype+'_bottom_counts_sensorid_'+string(bot_sensors[fix(isen)])+'_clean_sun_removed'+suffix, /rem)
            get_data, var_name3, data = d3, dlimits=dl3
            counts_arr[*,*,isen+n_elements(top_sensors)] = reform(d3.y)
            if iE[0] ne -1 then counts_arr[*,iE,isen+n_elements(top_sensors)] = !values.d_nan
            if bad_quality[0] ne -1 then counts_arr[bad_quality,*,isen+n_elements(top_sensors)] = !values.d_nan

          endif
        endif
        ;;;
      endfor
    endif else begin
      dalleyes = dblarr(n_elements(d.x), n_elements(d.v), n_elements(top_sensors))+!values.d_nan
      for isen = 0, n_elements(top_sensors)-1 do begin
        ; Only Top units in SITL product:
        var_name = strcompress(prefix+data_rate+'_'+level+'_'+datatype+'_top_'+data_units+'_sensorid_'+string(top_sensors[fix(isen)])+'_clean_sun_removed'+suffix, /rem)
        get_data, var_name, data = d, dlimits=dl
        dalleyes[*,*,isen] = reform(d.y)
        iE = where(abs(energies - d.v) gt en_chk*energies) ; Check for energies beyond en_chk [%] of the corrected energy bin center and replace with NAN
        if iE[0] ne -1 then dalleyes[*,iE,isen] = !values.d_nan
        ;;Added 10-03-2025
        if level eq 'l2' then begin
          var_name2 = strcompress(prefix+data_rate+'_'+level+'_'+datatype+'_top_quality_indicator_sensorid_'+string(top_sensors[fix(isen)])+suffix, /rem)
          get_data, var_name2, data = d2, dlimits=dl2
          bad_quality=where(d2.y ge quality_flag)
          if bad_quality[0] ne -1 then dalleyes[bad_quality,*,isen] = !values.d_nan
          
          if get_err gt 0. then begin
            var_name3 = strcompress(prefix+data_rate+'_l1b_'+datatype+'_top_counts_sensorid_'+string(top_sensors[fix(isen)])+'_clean_sun_removed'+suffix, /rem)
            get_data, var_name3, data = d3, dlimits=dl3
            counts_arr[*,*,isen] = reform(d3.y)
            if iE[0] ne -1 then counts_arr[*,iE,isen] = !values.d_nan
            if bad_quality[0] ne -1 then counts_arr[bad_quality,*,isen] = !values.d_nan

          endif
        endif
        ;;;
      endfor
    endelse
  endif
;stop
  ; if no data found, just return
  if undefined(dalleyes) then return
  
  flux_omni = reform(average(dalleyes,3,/NAN))
  
  ; Added by DLT on 31 Jan 2017: set unique gain factors per spacecraft
  ; Electrons:
  if probe eq '1' && datatype eq 'electron' then flux_omni = flux_omni*eGfact[0]
  if probe eq '2' && datatype eq 'electron' then flux_omni = flux_omni*eGfact[1]
  if probe eq '3' && datatype eq 'electron' then flux_omni = flux_omni*eGfact[2]
  if probe eq '4' && datatype eq 'electron' then flux_omni = flux_omni*eGfact[3]
  ; Ions:
  if probe eq '1' && datatype eq 'ion' then flux_omni = flux_omni*iGfact[0]
  if probe eq '2' && datatype eq 'ion' then flux_omni = flux_omni*iGfact[1]
  if probe eq '3' && datatype eq 'ion' then flux_omni = flux_omni*iGfact[2]
  if probe eq '4' && datatype eq 'ion' then flux_omni = flux_omni*iGfact[3]

  
  newname = strcompress('mms'+probe+'_epd_feeps_'+data_rate+'_'+level+'_'+datatype+'_'+data_units+'_omni'+suffix, /rem)

  dl.cdf.vatt.catdesc = 'MMS' + probe + ' FEEPS omni-directional ' + datatype + ' ' + data_units
  
  store_data, newname, data={x: d.x, y: flux_omni, v: en_label}, dlimits=dl
  options, newname, spec=1, /ylog, /zlog ;, yticks=3, ystyle=1, zrange=[1., 1.e6];, no_interp=0, y_no_interp=0, x_no_interp=0

  options, newname, ysubtitle='[keV]', ztitle=units_label, ystyle=1, /default,  yrange = minmax(energies)
  append_array, tplotnames, newname
  
  ;added 10-03-2025 cgabrielse
  if get_err gt 0. then begin
    uncertainty_arr=(sqrt(counts_arr)/counts_arr)  ;This is the percent uncertianty for each [time,energy,eye]
    uncertainty_arr_sq=uncertainty_arr^2
    sqrt_uncertainty_arr_sum_sq=sqrt(total(uncertainty_arr_sq,3,/NAN))
    finite_divisor=finite(uncertainty_arr_sq) ;Find out how many eyes were used. For example, if only 2 eyes had counts, you want to average over 2 eyes, not over 18 eyes.
    divisor=total(finite_divisor,3)
    percent_uncertainty=sqrt_uncertainty_arr_sum_sq/divisor
    
   ; percent_uncertainty=reform(average(uncertainty_arr_sq,3,/NAN))
    newname2 = strcompress('mms'+probe+'_epd_feeps_'+data_rate+'_'+level+'_'+datatype+'_'+data_units+'_omni_percent_uncertainty'+suffix, /rem) ;not units of percent. So the value here of 0.1 is 10%. calculated above as sqrt(N)/N.
    store_data,newname2,data={x:d.x,y:percent_uncertainty,v:en_label},dlimits=dl
    options,newname2,spec=1,/ylog,/zlog
    options,newname2,ztitle='[Percent uncertainty/100]',ysubtitle='[keV]',ystyle=1,/default,yrange = minmax(energies)
  endif
end ; pro
