;+
;Procedure: thm_part_dist_array
;
;Purpose: Returns and array of particle distributions for use by THM_PART_SLICE2D.
;
;Details: This is basically a wrapper for THM_PART_DIST that will create an
;         array of particle distribution structures sorted by increasing time.
;         
;
;Keywords:
; FORMAT: A string denoting the data that is desired: options are:
;           'tha_peif': Full Esa mode data, ions, probe A
;           'tha_peef': Full Esa mode data, electrons, probe A
;           'tha_peir': Reduced Esa mode data, ions, probe A
;           'tha_peer': Reduced Esa mode data, electrons, probe A
;           'tha_peir': Burst Esa mode data, ions, probe A
;           'tha_peer': Reduced Esa mode data, electrons, probe A
;           'tha_psif': Full Sst mode data, ions, probe A
;           'tha_psef': Full Sst mode data, electrons, probe A
;           'tha_psir': Reduced Sst mode data, ions, probe A
;           'tha_pser': Reduced Sst mode data, electrons, probe A
;         For other probes, just replace 'tha' with the appropriate string,
;         'thb', 'thc', 'thd', 'the'.  If this is not set, then the string is
;         constructed from the type and probe keywords.
; TRANGE: (Optional) Time range of interest (2 element array, string or numerical).
;         If neither TRANGE nor 'timespan' is set the user will be prompted.
; TYPE: Four character string denoting the type of data that you need,
;       e.g., 'peif' for full mode esa data.
; PROBE: The THEMIS probe, 'a','b','c','d','e'.
; MAG_DATA: Tplot variable containing magnetic field data that you want added
;           to the dat structures.
; VEL_DATA: Tplot variable containing velocity data that you want added to dat
;           structure under the tag name 'velocity'.  Overrides VEL_AUTO keyword.
;           The velocity data will be intpolated to the cadence of the particle
;           distributions.
; VEL_AUTO: Automatically calculates velocity at each time sample using V_3D_NEW.PRO
;           and adds it to the dat structure with the tag name 'velocity'.
; SST_CAL: Flag to use new SST calibrations (BETA)
; 
;Examples:  
;           dist_array = thm_part_dist_array(format='thb_peib', trange=trange)
;           dist_array = thm_part_dist_array(probe='b',format='psif', trange=trange, $
;                                            vel_data='tplot_vel', $
;                                            mag_data='tplot_mag')
;
;Contamination Removal:
;
; 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
;
; 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.
; 
; Crib sheets for contamination removal with thm_part_getspec and thm_part_moments
; are listed below in 'See Also'.  The same keywords are valid for this code; a few
; are documented below.
; 
;Contamination/Background Keywords:
;
; The following keywords are passed to the appropriate data retrieval routines. 
; The resulting distributions will have the specified contamination removal applied.
; 
; ESA Keywords:
;   BGND_REMOVE: Flag to turn on ESA background removal.
;   BGND_TYPE: String naming removal type, e.g. 'angle','omni', or 'anode'.
;   BGND_NPOINTS: Number of lowest values points to average over when determining background.
;   BGND_SCALE: Scaling factor that the background will be multiplied by before it is subtracted.
; 
; SST Keywords:
;
;   Check the list of keywords in THM_SST_REMOVE_SUNPULSE for more.
;
;   MASK_REMOVE: Set this keyword to the proportion of values that must be 0 at 
;                all energies to determine that a mask is present. Generally .99 
;                or 1.0 is a good value. The mask is a set of points that are set 
;                to 0 on-board the spacecraft.  By default they will be filled by 
;                linear interpolation across phi.
;
;   METHOD_SUNPULSE_CLEAN:  set this to a string:  Either 'median' or 'spin_fit' or 'z_score_mod'
;              'median':  This will remove all points that are greater 
;                than 2.0 standard deviations from the median.By default they will be filled by a 
;                linear interpolation across the phi angle by default. 
;              'spin_fit':  This will remove all points that are greater
;                than 2.0 standard deviations from a spin fit across phi angle.  The equation used to
;                fit is A+B*sin(phi)+C*cos(phi). By default these points will be filled by a linear
;                interpolation across the phi angle. The fitting is done using the svdfit routine
;                from the idl distribution.
;              'z_score_mod': This will remove all points that have a modified z-score(calculated across phi) greater than 3.5 
;                The modified z-score is a normalized outlier detection test defined as follows:  
;                #1 X_Bar = median(X+1)
;                #2 Sigma = MAD = Median Absolute Deviation = median(abs(X-X_Bar))
;                #3 Z_Score_Mod = .6745*(X - X_Bar)/Sigma
;                This test can often get excellent results because it is insensitive to variation in standard deviation
;                and skew in the distributions.  
;   FILLIN_METHOD: Set this keyword to a string that specifies the method used to fill the points that are
;             removed via the method_sunpulse_clean or the mask_remove keywords.
;             If 'interpolation' is set, this routine will interpolate across the phi angle.  This is the 
;               default behavior. Interpolation is done using the interp_gap routine.
;             If 'spin_fit' is set this routine will perform a spin fit to the data after the points
;               have been removed using the equation A+B*sin(phi)+C*cos(phi).  It will then generate
;               expected values for each removed phi using the equation of fit.   The fitting is done using
;               the svdfit routine from the idl distribution.  Note that if 'spin_fit' is selected for
;               the clean method and the fill method, this routine will perform two spin fits.
;   LIMIT_SUNPULSE_CLEAN: set this equal to a floating point number that will override the default of 2.0 standard
;             deviation tolerance or 3.5 z_score_tolerance, used by the sunpulse cleaning methods by default.
;             This keyword will only have an effect if the method_sunpulse_clean keyword is set.

;   
;   ENOISE_BINS: A 0-1 array that indicates which bins should be used to calculate 
;               electronic noise.  A 0 indicates that the bin should be used for 
;               electronic noise calculations.  This is basically the output from 
;               the bins argument of edit3dbins. It should have dimensions 16x4.
;   ENOISE_REMOVE_METHOD: (default: 'fit_median') set the keyword to a string specifying the method you want to use to calculate the electronic noise that will be subtracted 
;               This function combines values across time.  The allowable options are:
;                 'min':  Use the minimum value in the time interval for each bin/energy combination.
;                 'average': Use the average value in the time interval for each bin/energy combination.
;                 'median': Use the median value in the time interval for each bin/energy combination.
;                 'fit_average': Fill in selected bins with a value that is interpolated across phi then subtracts the average of the difference
;                                between the interpolated value and the actual value from each selected bin/energy combination.
;                 'fit_median': Fill in selected bins with a value that is interpolated across phi then subtracts the median of the difference
;                               between the interpolated value and the actual value from each selected bin/energy combination.
;                 'fill': Fill the selected bins using the technique specified by enoise_remove_method_fit, or interpolation by default.
;                         (note that if this method is used, enoise_bgnd_time is not required)
;   ENOISE_REMOVE_FIT_METHOD: (default:'interpolation'):  Set this keyword to control the method used in 'fit_average' & 'fit_median' to
;                    fit across phi. Options are 'interpolation' & 'spin_fit' By default, missing bins are interpolated across phi.  Setting
;                    enoise_remove_fit_method='spin_fit' will instead try to fill by fitting to a curve of the form A+B*sin(phi)+C*cos(phi).
;   ENOISE_BGND_TIMES: This should be either a 2 element array or a 2xN element array(where n is the number of elements in enoise_bins).  
;                    The arguments represents the start and end times over which the electronic background will be calculated for each
;                    bin.  If you pass a 2 element array the same start and end times can be used for each bin.  If you pass a 2xN element
;                    array, then the Ith bin in enoise_bins will use the time enoise_bgnd_time[0,I] as the start time and enoise_bgnd_time[1,I] as
;                    the end time for the calculation of the background for that bin.  If this keyword is not set then electronic noise will not 
;                    be subtracted.
;
;
;See Also: thm_part_slice2d, thm_clib_part_slice2d, 
;          thm_crib_esa_bgnd_remove, thm_esa_bgnd_remove, 
;          thm_crib_sst_contamination, thm_sst_remove_sunpulse
;
;
;Created by Bryan Kerr
;Modified by A. Flores
;-

function thm_part_dist_array, format=format, trange=trange, type=type, $
                              probe=probe, mag_data=mag_data, vel_data=vel_data, $
                              vel_auto=vel_auto, suffix=suffix, err_msg=err_msg, $
                              gettimes=gettimes, sst_cal=sst_cal, $
                              _extra = _extra 

  compile_opt idl2


if keyword_set(format) then begin
  format = strlowcase(format)
  probe = strmid(format,2,1)
  dtype  = strmid(format,4,4)
endif else if keyword_set(probe) && keyword_set(type) then begin
  probe = strlowcase(probe)
  dtype = strlowcase(type)
  format = 'th'+probe+'_'+dtype
endif else begin
  dprint, 'Either FORMAT or PROBE & TYPE keywords must be specified."
  return, -1
endelse

inst = strmid(dtype,1,1)
species = strmid(dtype,2,1)

; check requested probe
dummy = where(probe eq ['a','b','c','d','e'], yes_probe)
if yes_probe lt 1 then begin
  dprint, 'Invalid probe: ' + probe
  return, -1
endif

; check requested instrument type
if inst ne 'e' && inst ne 's' then begin
  dprint, 'Invalid instrument type: ' + inst
  return, -1
endif

; check requested species
if species ne 'e' && species ne 'i' then begin
  dprint, 'Invalid species: ' + species
  return, -1
endif

if ~keyword_set(vel_data) then vel_auto = 1b

; copy time range to new variable(s) and make sure timespan gets set
if keyword_set(trange) then begin
   trd = minmax(time_double(trange))
   tr = time_string(trd)
   ndays = (trd[1] - trd[0]) / 86400
   timespan,trd[0],ndays
endif else begin
   trd = timerange()
   tr = time_string(trd)
endelse


; load L1 data
if inst eq 'e' then begin
  thm_load_esa_pkt, probe=probe, trange=tr, datatype=dtype, suffix=suffix
endif else begin
  if keyword_set(sst_cal) then begin
    thm_load_sst2, probe=probe, trange=tr, datatype=dtype, suffix=suffix
  endif else begin
    thm_load_sst, probe=probe, trange=tr, datatype=dtype, suffix=suffix
  endelse
endelse


; get time indexes of data in requested time range
times = thm_part_dist(format, /times, sst_cal=sst_cal)
if size(times,/type) eq 8 then begin
  err_msg = 'Unable to retrieve times for th'+probe+'_'+inst+ $
            ' between '+tr[0]+ $
            ' and '+tr[1]+'.' 
  dprint, err_msg
  return, -1
endif


;time correction to point at bin center is applied for ESA, but not for SST
if strmid(inst,1,1) eq 's' then begin
  times += 1.5
endif

if keyword_set(gettimes) then return, times


; check that data exists in requested range
time_ind = where(times ge trd[0] and times le trd[1], n_times) 
if (size(times,/type) ne 5) or (n_times lt 1) then begin
  err_msg = 'No '+format+' data for time range '+tr[0]+ $
             ' to '+tr[1]+'.' 
  dprint, err_msg
  return, -1
endif 


;interpolate velocity data
if keyword_set(vel_data) then begin

  tinterpol_mxn, vel_data, times[time_ind], /nan_extrapolate, error=success

  if success then begin
    get_data, vel_data+'_interp', data=d
    vel = d.y
  endif else begin
    err_msg = 'Error interpolating velocity data. Not adding velocity data to dat structure.'
    dprint, err_msg
    vel_data = 0
  endelse
endif


enoise_tot = thm_sst_erange_bin_val('th'+probe, inst, times, _extra=_extra)
mask_tot = thm_sst_find_masking('th'+probe, inst, time_ind, _extra=_extra)


; Loop to create array
for i=0L,n_times-1 do begin


  dat = thm_part_dist(format, index=time_ind[i], sst_cal=sst_cal, _extra=_extra)


  ; add mag data to dat structure
  if keyword_set(mag_data) then add_magf, dat, mag_data
  
  ; add velocity data to dat structure
  if keyword_set(vel_data) then begin
    ; add user-input velocity data to data structure
    add_str_element, dat, 'velocity', reform(vel[i,*])    
  endif else if keyword_set(vel_auto) then begin
    ; calculate velocity in km/s and add to dat stucture
    vel=v_3d_new(dat);*1000.
    add_str_element, dat, 'velocity', vel
  endif

  
  ; append current dat structure to array of dat structures
  if keyword_set(dist_arr) then begin
    if ~array_equal( size(dist_arr[0].bins,/dim), size(dat.bins,/dim) ) then begin
      err_msg = 'Mode change encountered at '+time_string(dat.time)+'. Data after this point must be loaded separately.'
      dprint, err_msg
      return, dist_arr
    endif
    dist_arr = [dist_arr, dat]
  endif else begin
    dist_arr = dat
  endelse

endfor

return, dist_arr
end