;+
;PROCEDURE: THM_SST_ERANGE_BIN_VAL
;
;Purpose:
;   This routine generates the values that will be subtracted to remove electronic noise.  The actual subtraction is done by thm_sst_remove_sunpulse.pro
;   A separate value will be calculated for each combination of bin and energy, but a single value will be calculated for all times(at a particular bin/energy).  
;   This value will be generated using one of several different functions of the values across time.  The functions can be specified by the user. 
;   The user should never call this routine directly, but should instead provide the appropriate arguments to thm_part_moments,thm_part_moments2, & thm_part_getspec
;   These routines will guarantee that this routine is called correctly.
;   
;Arguments(note these arguments are generated correctly in thm_part_moments & thm_part_moments2)
;   thx:  a string representing a probe prefix(ie 'tha')
;
;   instrument:  a string representing the instrument(ie 'psif')  Note that this routine will only perform an operation if
;   it is passed a string representing sst full distribution data.  (ie 'psif', 'psef')
;
;   times:  a list of times for the sst measurements, this list is generated by thm_part_dist
;
;Keywords:
;
; method_clean: If set to 'manual', this keyword disables the enoise code, as another method will be used by thm_sst_remove_sunpulse.pro
;
; 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_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.
;
; 
;
; 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).
;                      
;        
;SEE ALSO:
;  thm_part_moments.pro, thm_part_moments2.pro, thm_part_getspec.pro
;  thm_part_dist.pro, thm_sst_psif.pro, thm_sst_psef.pro,thm_crib_sst_contamination.pro
;  thm_sst_find_masking.pro, thm_sst_remove_sunpulse.pro
;
; $LastChangedBy: pcruce $
; $LastChangedDate: 2014-01-03 12:10:47 -0800 (Fri, 03 Jan 2014) $
; $LastChangedRevision: 13737 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/spdsoft/tags/spedas_1_00/projects/themis/spacecraft/particles/SST/thm_sst_erange_bin_val.pro $
;
;
;-

; Moving to separate file
;;HELPER function
;;take the dat component of a structure and splits it into an array
;;ordered in terms of theta =  energy*angle->energy*theta*phi
;;dimensions 16*64->16*4*16,  phi is guaranteed to be contiguous but
;;not necessarily ascending(some phi may be out of phase by 180 degrees)
;;returns indices to perform this transformation
;function dat2angsplit,dat
;
;  compile_opt idl2,hidden
;  
;  index = indgen(16,64)
;  t_sort = bsort(dat.theta[0,*])
;  index = index[*,t_sort]
;  return,transpose(reform(index,16,16,4),[0,2,1]) ; this magic properly reforms and orders the dimensions  
;
;end

;HELPER function
;performs filling using an svd fit
function do_enoise_spin_fill,dat,ref,phi

  compile_opt idl2,hidden

  dat2 = dat

  good = where(finite(dat2))
  bad = ssl_set_complement(good,indgen(n_elements(dat)))
  
  if good[0] ne -1 && bad[0] ne -1 then begin

    idx = where(dat2[good] ne 0)
    
    if idx[0] eq -1 then begin  ;if the system of equations is not-invertible because it is all zeros then add a constant offset
      dat2[good] += 1
    endif
      
    except_val = !EXCEPT  ;svdfit annoyingly reports a floating point error EVERY TIME IT OCCURS if you don't disable them entirely
    !EXCEPT = 0  
    result = svdfit(phi[good],dat2[good],a=[1,1,1],function_name='svd_func',status=stat) ;fit the data
    !EXCEPT = except_val
                     
    if stat eq 0 then begin
    ; replace values with fit points
      dat2[bad] = result[0] + result[1]*sin(phi[bad]*!DTOR)+result[2]*cos(phi[bad]*!DTOR)         
    endif else begin
      dat2[bad] = ref[bad] ;replace values with ref values to ensure there are no NaNs in the output
    endelse
      
    
    if idx[0] eq -1 then begin ;remove a constant offset if one was added
      dat2 -= 1
    endif
  endif
  
  idx = where(~finite(dat2))
  
  return,dat2

end

function thm_sst_erange_bin_val,thx,instrument,times,enoise_bins=enoise_bins,enoise_bgnd_times=enoise_bgnd_time,enoise_remove_method=enoise_remove_method,enoise_remove_fit_method=enoise_remove_fit_method,method_clean=method_clean,limit_sunpulse_clean=limit_sunpulse_clean,_extra=ex

  compile_opt idl2 
 
  ;Validate inputs
  
  if keyword_set(method_clean) && strlowcase(method_clean) eq 'manual' then begin
    return,-1
  endif else if (keyword_set(method_clean) && strlowcase(method_clean) eq 'automatic') then begin
   
    if ~strmatch(instrument,'ps?f') && ~strmatch(instrument,'pseb') then begin
      dprint,'Automatic sun contamination  not enabled for instrument ' + instrument
      return,-1
    endif
   
    trange = times ; note that this is not nec the user specified time range. Is at least a full day.
    n_tm = n_elements(trange)
    format = thx+'_'+instrument
   
    if keyword_set(limit_sunpulse_clean) && is_num(limit_sunpulse_clean) then begin
      zscore_max = limit_sunpulse_clean
    endif else begin
      zscore_max = 10 ;In standard statistical tests, this is an extremely large limit, making for a weak outlier test.  
                      ;But for sun contamination outliers generally have gigantic scores, so it shouldn't be a problem.
                      ;What tricks the test isn't the threshold, but the heterogenous nature of the data across time.
                      ;A better method in the future would test attenuator open and closed data separately
    endelse
     
    totalzscore = dblarr(64)
    ; The modified z score for each angular bin (summed over energy) is calculated for each time.
    ; The truncated mean of these z scores over time is then used to find outliers.
    ; The test is applied for all energies summed and for bin 9 separately.  
    ; For some reason, there is sometimes a pattern of contamination the hits bin 9 only

    data_arr = dblarr(n_tm,16,64)
    
    for j = 0, n_tm-1 do begin ; loop over times
      dat = thm_part_dist(format,index=j,_extra=ex,/no_sun)  ;get data for this time
      if is_struct(dat) then begin
        data_arr[j,*,*]=dat.data/thm_sst_atten_scale(dat.atten,dimen(dat.data))
      endif
    endfor
    
    totaloverenergy = total(data_arr,2)
    ;totaloverenergy = reform(data_arr,n_tm,16*64)
    zscr_med = median(totaloverenergy+1,dimension=2); Adding 1 here prevents infinities in the case that all totaloverenergy=0
    ;zscr_sigma = median(abs(totaloverenergy-(zscr_med#(dblarr(64)+1))),dimension=2)
    zscr_sigma = median(abs(totaloverenergy-(zscr_med#(dblarr(64)+1))),dimension=2)
    zscr_arr = .6745*(totaloverenergy-zscr_med#(dblarr(64)+1))/(zscr_sigma#(dblarr(64)+1))
    ;totalzscore = total(zscr_arr,1)
    
    lim1=0.1 ;all energies
    lim2=0.1 ;bin 9
   
    ;sort so that we can do a truncated mean
    ;This removes any outliers in the z-score days due to glitches
    ;Glitches can create z-scores that are big enough that they throw off an entire day's average even though they represent only a few points out of 10-20k points
    for j = 0,64-1 do begin
      sort_idx=bsort(zscr_arr[*,j])
      zscr_arr[*,j]=zscr_arr[sort_idx,j]  
    endfor
    
    low_ten=floor(n_tm*lim1)
    high_ten=floor(n_tm*(1-lim1))
    
    ;truncated mean removes outliers due to glitches.
    averagezscore = total(zscr_arr[low_ten:high_ten,*],1,/nan)/(high_ten-low_ten)
    ;averagezscore=total(zscr_arr,1,/nan)/n_tm
        
    ;averagezscore = totalzscore/n_tm
    zscore_bad = where(averagezscore gt zscore_max,c); not using abs value as we are really interested in values too high, not too low (masking is removed separately)
    out = dblarr(16,64)
    if c gt 0 then begin
      dprint,'Sun Contaminated Bins Flagged: ' + strjoin(strtrim(zscore_bad,2),' '),dlevel=2
      out[*,zscore_bad] = !VALUES.D_NAN ;store result
    endif
    
    ;We do a separate test on energy bin 9.  For a reason that is not understood, there is a pattern of sun contamination that only affects bin 9 on occassion.
    ;Possible explanations include idiosyncrasies of the electronics & glint off of the side of the collimator.
    totaloverenergy=reform(data_arr[*,9,*])
    ;totaloverenergy = reform(data_arr,n_tm,16*64)
    zscr_med = median(totaloverenergy+1,dimension=2); Adding 1 here prevents infinities in the case that all totaloverenergy=0
    zscr_sigma = median(abs(totaloverenergy-(zscr_med#(dblarr(64)+1))),dimension=2)
    zscr_arr = .6745*(totaloverenergy-zscr_med#(dblarr(64)+1))/(zscr_sigma#(dblarr(64)+1))
    ;totalzscore = total(zscr_arr,1)
    
    ;sort so that we can do a truncated mean
    for j = 0,64-1 do begin
      sort_idx=bsort(zscr_arr[*,j])
      zscr_arr[*,j]=zscr_arr[sort_idx,j]  
    endfor
    
    low_ten=floor(n_tm*lim2)
    high_ten=floor(n_tm*(1-lim2))
    
    ;truncated mean removes outliers due to glitches.
    averagezscore = total(zscr_arr[low_ten:high_ten,*],1,/nan)/(high_ten-low_ten)
    ;averagezscore=total(zscr_arr,1,/nan)/n_tm
        
    ;averagezscore = totalzscore/n_tm
    zscore_bad = where(averagezscore gt zscore_max,c); not using abs value as we are really interested in values too high, not too low (masking is removed separately)
    if c gt 0 then begin
      dprint,'Earth Shine Contaminated Bins Flagged: ' + strjoin(strtrim(zscore_bad,2),' '),dlevel=2
      out[*,zscore_bad] = !VALUES.D_NAN ;store result
    endif

    return,out 
   
  endif else if keyword_set(method_clean) then begin
    dprint,'Unrecognized value for keyword method_clean: ' + method_clean
    dprint,'Did you mean to use method_sunpulse_clean?"
  endif 
  
  
; enoise_bins rather than automatic clean

  if ~keyword_set(enoise_bins) then begin
    return,-1
  endif

  if ~strmatch(instrument,'ps?f') && ~strmatch(instrument,'pseb') then begin
    dprint,'Electronic removal not enabled for instrument ' + instrument
    return,-1
  endif

  ;fix to use 0-1 array rather than indexes
  enoise_bins2 = where(enoise_bins eq 0,bin_cnt)
  
  if enoise_bins2[0] eq -1 then begin
    dprint,'No bins selected for ' + thx + '/' + instrument + '  Enoise removal not applied'
    return,-1
  endif
  
  if keyword_set(enoise_remove_method) && strlowcase(enoise_remove_method) eq 'fill' then begin
 
    dprint,'Enoise Filling ' + strtrim(bin_cnt,2) + ' angle bins for ' + thx + '/' + instrument
    out = dblarr(16,64) ;output array
    out[*,enoise_bins2] = !VALUES.D_NAN
    return,out
  
  endif
  
  if ~keyword_set(enoise_bgnd_time) then begin
    dprint,'Cannot apply selected enoise remove method for '+ thx + '/' + instrument + '.  Enoise_bgnd_time not set'
  endif
  
  d = dimen(enoise_bgnd_time)
  
  if d[0] ne 2 then begin
    dprint,'enoise_bgnd_time has illegal dimensions, electronic noise removal will not be performed'
    return,-1
  endif
 
  n_bin = n_elements(enoise_bins2)
  
  if n_elements(d) eq 1 then begin 
    e_times = make_array(2,n_bin,type=size(enoise_bgnd_time,/type))
    e_times[0,*] = enoise_bgnd_time[0]
    e_times[1,*] = enoise_bgnd_time[1]
  endif else begin
    e_times = enoise_bgnd_time
  
    if d[1] ne n_bin then begin
      dprint,'number of enoise_bins does not match number of enoise_times, electronic noise removal will not be performed'
      return,-1
    endif
    
  endelse
    

  
  ;convert to numeric format if input is not numeric
  if is_string(e_times) then begin
    e_times = time_double(e_times)
  endif
  
  ;Identify method
  if ~keyword_set(enoise_remove_method) then begin
    enoise_remove_method = 'fit_median'
  endif
  
  if strlowcase(enoise_remove_method) eq 'fit_median' then begin
    fit_median = 1
  endif else if strlowcase(enoise_remove_method) eq 'fit_average' then begin
    fit_average = 1
  endif else if strlowcase(enoise_remove_method) eq 'median' then begin
    median = 1
  endif else if strlowcase(enoise_remove_method) eq 'average' then begin
    average = 1
  endif else if strlowcase(enoise_remove_method) eq 'min' then begin
    min = 1
  endif else begin
    dprint,'enoise remove method "' + enoise_remove_method + '" unrecognized, electronic noise removal will not be performed'
    return,-1
  endelse
  
  if (keyword_set(min) || keyword_set(average) || keyword_set(min)) && $
     keyword_set(enoise_remove_fit_method) then begin
     dprint,'Enoise fit method does not have an effect if remove_method is not fit_average,fit_median, or fill'
  endif 
  
  if ~keyword_set(enoise_remove_fit_method) then begin
    enoise_remove_fit_method = 'interpolation'
  endif
  
  if strlowcase(enoise_remove_fit_method) eq 'interpolation' then begin
    interpol_fit = 1
  endif else if strlowcase(enoise_remove_fit_method) eq 'spin_fit' then begin
    spin_fit = 1
  endif else begin
    dprint,'enoise fit method "' + enoise_remove_fit_method + '" unrecognized, electronic noise removal will not be performed'
    return,-1
  endelse
 
  
  format = thx+'_'+instrument
  
  out = dblarr(16,64) ;output array
  
  for i = 0,n_bin-1 do begin ; loop over angle bins
  
    trange = e_times[*,i]  ; get times and bins
    ind = where(times ge trange[0] and times le trange[1],cnt)      
    bin = enoise_bins2[i]
      
    if cnt eq 0 then begin
      dprint,'No valid times found for bin: ' + strcompress(string(bin))
      continue
    endif
      
    n_tm = n_elements(ind)  ;loop upper bound
    totd = dblarr(16,n_tm)   ;temporary storage of processed values
    totm = dblarr(16,n_tm)
      
;    print,'TIME!!! ',n_tm
;    print,ind
;    print,time_string(times[ind])
      
    for j = 0, n_tm-1 do begin ; loop over times
      
      dat = thm_part_dist(format,index=ind[j],_extra=ex)  ;get data for this time
      
      ;copy of Jim M's quick fix, to guarantee time matchup
      If(size(dat, /type) Eq 8) Then Begin
        mid_time = (dat.time+dat.end_time)/2.0
        trange2 = [dat.time, dat.end_time]             
        dat.time = mid_time
        str_element, dat, 'end_time', /delete
        str_element, dat, 'trange', trange2, /add_replace
      Endif
      
      if is_struct(dat) then begin
           
        if keyword_set(fit_median) || keyword_set(fit_average) then begin
       
          ang_split = thm_sst_dat2angsplit(dat)   ;indexes to turn the data from 16x64 to 16x4x64
          edat = dat
          edat.data[*,enoise_bins2] = !VALUES.F_NAN ;these bins are going to be interpolated over 
          data = dat.data[ang_split]
          edata = edat.data[ang_split]
          phi = dat.phi[ang_split]
          
          thetas = dat.theta[uniq(dat.theta,sort(dat.theta))]    
          theta = where(dat.theta[0,bin] eq thetas)    ;which theta is the bin in question in 
          
          phi[*,[0,2],0:7] -= 360  ; make phi's monotonic      
  
          for k = 0, 16-1 do begin  ;loop over energy
          
             phis = reform(phi[k,theta,*])
             edatas = reform(edata[k,theta,*])
             ref_datas = reform(data[k,theta,*])
             
             triple_phi = [phis-360,phis,phis+360]
             triple_dat = [edatas,edatas,edatas]
             triple_ref = [ref_datas,ref_datas,ref_datas]
             
             if keyword_set(interpol_fit) then begin
               interp_gap, triple_phi, triple_dat
             endif else begin
               triple_dat = do_enoise_spin_fill(triple_dat,triple_ref,triple_phi)
             endelse
             ;triple_dat = do_spin_fill(triple_dat,triple_ref,triple_phi)          
       
             edata[k,theta,*] = triple_dat[16:31]
             ;data[k,theta,*] = edatas
         endfor
         
         ;replace data
         inv_ang_split = bsort(ang_split)
         edat.data = edata[inv_ang_split]        
         totm[*,j] = edat.data[*,bin]
         
        endif
       
        totd[*,j] = dat.data[*,bin]  ;store the data for this part of the time interval
     
      endif
    endfor
    
    if keyword_set(fit_median) then begin  ;perform requested processing operation  
      out[*,bin] = median(totd,dimension=2) - median(totm,dimension=2)
    endif else if keyword_set(fit_average) then begin
      out[*,bin] = average(totd,2,/nan) - average(totm,2,/nan)
    endif else if keyword_set(median) then begin
      out[*,bin] = median(totd,dimension=2)
    endif else if keyword_set(average) then begin
      out[*,bin] = average(totd,2,/nan)
    endif else if keyword_set(min) then begin
      out[*,bin] = min(totd,dimension=2,/nan)
    endif
    
  endfor
  
;  idx = where(out lt 0,c)
;  
;  if c gt 0 then begin
;    out[idx] = 0
;  endif
   
  return,out
   
end