;+
;Procedure:
;  spd_pgs_make_phi_spec
;
;Purpose:
;  Builds phi (longitudinal) spectrogram from simplified particle data structure.
;
;
;Input:
;  data: single sanitized data structure
;  
;
;Input/Output:
;  spec: The spectrogram (ny x ntimes)
;  yaxis: The y axis (ny OR ny x ntimes)
;  resolution: (optional) Specify output resolution
;  
;  -Each time this procedure runs it will concatenate the sample's data
;   to the SPEC variable.
;  -Both variables will be initialized if not set
;  -The y axis will remain a single dimension until a change is detected
;   in the data, at which point it will be expanded to two dimensions.
;
;
;Notes:
;  -Code for original value_locate() based version remains commented out
;   below. It should produce identical spectrograms for regular phi grids.
;
;
;$LastChangedBy: egrimes $
;$LastChangedDate: 2021-04-20 16:08:59 -0700 (Tue, 20 Apr 2021) $
;$LastChangedRevision: 29895 $
;$URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/spdsoft/trunk/general/science/spd_part_products/spd_pgs_make_phi_spec.pro $
;-

pro spd_pgs_make_phi_spec, data, spec=spec, sigma=sigma, yaxis=yaxis, resolution=resolution, _extra=ex

    compile_opt idl2, hidden
  
  
  if ~is_struct(data) then return
  
  
  dr = !dpi/180.
  rd = 1/dr

  ;copy data and zero inactive bins to ensure
  ;areas with no data are represented as NaN
  d = data.data
;  scaling = data.scaling
  idx = where(~data.bins,nd)
  if nd gt 0 then begin
    d[idx] = 0.
  endif
  

  ;determine number of energies
  n_energy = (size(/dim, data.data))[0]

  ;determine number of phi bins at small theta for an energy
  ; -number of phis per theta decreases at higher latitudes for ESA
  ; -using the max number across phi should allow for equal 
  ;  statistics across phi bins
  ; -this assumes the number does not change across energy
  if undefined(resolution) then begin
    dummy = min(abs(data.theta[0,*]),tminidx)
    dummy = where(data.theta[0,*] eq (data.theta[0,*])[tminidx], n_phi)
  endif else begin
    n_phi = resolution
  endelse
  

  ;init this sample's piece of the spectrogram
  ave = replicate(!values.f_nan,n_phi)
  ave_s = ave

  ;form grid specifying the spectrogram's phi bins
  phi_grid = interpol([0,360.],n_phi+1)
  phi_grid_width = median(phi_grid - shift(phi_grid,1))
  
  
;  phi_idx = value_locate(phi_grid, data.phi)
;  phi_half_width = ceil(0.5*data.dphi / phi_grid_width)  
   

  ;get min/max of all data bins
  ;keep phi in [0,360]
  phi_min = (data.phi - 0.5*data.dphi) 
  phi_max = (data.phi + 0.5*data.dphi) mod 360.

  ;algorithm below assumes maximums at 360 not wrapped to 0
  zm = where(phi_max eq 0, nzm)
  if nzm gt 0 then phi_max[zm] = 360
  
  ;keep phi in [0,360]
  ltz = where(phi_min lt 0, nltz)
  if nltz gt 0 then phi_min[ltz] += 360


  ;keep track of bins that span phi=0
  wrapped = phi_min gt phi_max
  
  
  ;When averaging data bins will be weighted by the solid angle of their overlap,
  ;with the given spectrogram bin.  Since each spectrogram bin spans all theta 
  ;values the theta portion of that calculation can be done in advance.  These
  ;values will later be multiplied by the overlap along phi to get the total 
  ;solid angle. 
  omega_part = abs(sin( dr * (data.theta + .5*data.dtheta) ) - $
                   sin( dr * (data.theta - .5*data.dtheta) )    )   
  
  
  ;Loop over each phi bin in the spectrogram and determine which data bins
  ;overlap.  All overlapping bins will be weighted according to the solid 
  ;angle of their intersection and averaged.
  for i=0, n_phi-1 do begin


;    idx = where( (phi_idx eq i) ,n)


    weight = fltarr(size(/dim,phi_min))
    
    ;data bins whose maximum overlaps the current spectrogram bin
    idx_max = where(phi_max gt phi_grid[i] and phi_max lt phi_grid[i+1], nmax)
    if nmax gt 0 then begin
      weight[idx_max] = (phi_max[idx_max] - phi_grid[i]) * omega_part[idx_max]
    endif
    
    ;data bins whose minimum overlaps the current spectrogram bin
    idx_min = where(phi_min gt phi_grid[i] and phi_min lt phi_grid[i+1], nmin)
    if nmin gt 0 then begin
      weight[idx_min] = (phi_grid[i+1] - phi_min[idx_min]) * omega_part[idx_min]
    endif
    
    ;data bins contained withing the current spectrogram bin
    ;check for single phi data (phi_min=phi_max)
    contained = ssl_set_intersection(idx_max,idx_min)
    if contained[0] ne -1 && n_phi gt 1 then begin
      weight[contained] = data.dphi[contained] * omega_part[contained]
    endif
    
    ;data bins which completely cover the current spectrogram bin 
    idx_all = where( phi_min le phi_grid[i] and phi_max ge phi_grid[i+1] $
          or (wrapped and phi_min gt phi_grid[i+1] and phi_max gt phi_grid[i+1]) $
          or (wrapped and phi_min lt phi_grid[i]   and phi_max lt phi_grid[i]), nall)
    if nall gt 0 then begin
      weight[idx_all] = phi_grid_width * omega_part[idx_all]
    endif

    ;combine indices 
    idx = ssl_set_union(idx_min,idx_max)
    idx = ssl_set_union(idx,idx_all)
    
    ;assign a weighted average to this bin
    if (nmax + nmin + nall) gt 0 then begin
    
;      ave[i] = total(d[idx]) / total(data.bins[idx])

      ;indices will contain a -1 if any other the searches failed 
      if idx[0] eq -1 then begin
        idx = idx[1:n_elements(idx)-1]
      endif
      
      ;normalize weighting to selected, active bins
      weight = weight[idx] * data.bins[idx]
      weight = weight / total(weight)
      
      ;average
      ave[i] = total(d[idx] * weight)
      
      ;standard deviation
      ave_s[i] = sqrt(  total(d[idx] * data.scaling[idx] * weight^2)  )
      
    endif else begin
      ;nothing
    endelse
  
  endfor

  ;------------------------------
  
  
  ;get y axis
  y = ( phi_grid + shift(phi_grid,1) ) / 2.
  y = y[1:n_phi]
  
  ;testing kludge for single-phi data
  if n_phi eq 1 then begin
    y = [0.,360]
    ave = [ave,ave]
    ave_s = [ave_s,ave_s]
  endif
  
  ;concatenate y axes
  if undefined(yaxis) then begin
    yaxis = y
  endif else begin
    spd_pgs_concat_yaxis, yaxis, y, ns=dimen2(spec)
  endelse

  
  ;concatenate spectra
  if undefined(spec) then begin
    spec = temporary(ave)
  endif else begin
    spd_pgs_concat_spec, spec, ave
  endelse
  
  
  ;concatenate standard deviation
  if undefined(sigma) then begin
    sigma = temporary(ave_s)
  endif else begin
    spd_pgs_concat_spec, sigma, ave_s
  endelse 
  
end