;+
; 
; Name: thm_part_slice2d_geo.pro
; 
; Purpose:  Helper function for thm_part_slice2d.pro
;           Produces slice using each bin's boundaries.
;
; Arguments:
;   datapoints: N element array of data values
;   velocities: Nx3 element array of cartesian velocities
;   resolution: Single value (R) giving the number of points in each 
;               dimension of the slice
;   dp: N element array of phi ranges
;   dt: N element array of theta ranges
;   dv: N element array of velocity ranges
;   slice orient: 4 Element array, the first 3 elements are the 
;                 slice plane's normal vector, the last is the 
;                 displacement from zero along that vector.
;   
; Input Keywords:
;   ct: Rotation matrix from DSL -> desired coordinates
;   rot: Rotation matrix from given coordinates to specific rotation
;        (e.g. BV, perp_xy)
; 
; Output Keywords:
;   xgrid: R element array of x-axis values for the slice
;   ygrid: R element array of y-axis values for the slice
;   part_slice: RxR element array containing the slice data
; 
; Caveats: This routine will slow as the number of bins (N) increases.
; 
;-
pro thm_part_slice2d_geo, datapoints, velocities, $
                          resolution, dp, dt, dv, $
                          slice_orient, ct=ct, rot=rot, $
                          part_slice=part_slice, $
                          xgrid=xgrid, ygrid=ygrid, $
                          fail=fail

    compile_opt idl2, hidden

  n = resolution
  rd = 180./!dpi
  tolerance = 5d-7


  ;Get slice orienation vectors
  thm_part_slice2d_orientslice, slice_orient, fail=fail, $
                                displacement = displacement, $
                                norm = norm, xvec = xvec
  if keyword_set(fail) then return


  ;Get rotation matrix for custom orienation
  m = thm_part_slice_trans(norm, xvec, fail=fail)
  if keyword_set(fail) then return


  ;Convert bin data to spherical coords
  v = sqrt( total(velocities^2,2) )
  phi = reform(rd * atan(velocities[*,1], velocities[*,0]))
  theta = reform(rd * atan(velocities[*,2], sqrt(total(velocities[*,0:1]^2,2)) ))


  ;Create coordinates for and initialize slice
  vrange = max(abs(v))*[-1,1]
  xgrid = interpol(vrange + [-1,1]*max(dv), n)
  ygrid = interpol(vrange + [-1,1]*max(dv), n)
  part_slice = dblarr(n,n)


  ;Create grid of coordinates for entire slice plane
  u = [ [reform(xgrid # replicate(1.,n), n^2)], $         ;x
        [reform(replicate(1.,n) # ygrid, n^2)], $         ;y
        [reform(replicate(displacement,[n,n]), n^2)]   ]  ;z


  ;Rotate slice coordinates to desired location. The matrices 
  ;are meant to transform INTO the slice plane's coordinatess.  
  ;to find the slice's location in the data's native coordinates 
  ;we combine the rotations then reverse.
  if keyword_set(ct) then begin
    if keyword_set(rot) then m = (ct # rot) # m $
      else m = ct # m
  endif else begin
    if keyword_set(rot) then m = rot # m
  endelse
  u = u # transpose(m)


  ;Convert transformed slice coordinates to spherical
  pcoords = rd * atan(u[*,1],u[*,0])                      ;phi
  tcoords = rd * atan(u[*,2], sqrt(total(u[*,0:1]^2,2)) ) ;theta
  vcoords = sqrt(u[*,0]^2 + u[*,1]^2 + u[*,2]^2)          ;r



  ; Loop over bins to determine what region
  ; each bin covers on the slice plane.
  weight = bytarr(size(pcoords,/dim))
  for i=0, n_elements(datapoints)-1 do begin
    if datapoints[i] eq 0 then continue

    ; Theta--
    tlim = [ theta[i]-0.5*dt[i], theta[i]+0.5*dt[i] ]

    tr = where( abs(tlim - round(tlim)) lt tolerance, ntr)
    if ntr gt 0 then tlim[tr] = round(tlim[tr])
    
    t = (tcoords gt tlim[0]) and (tcoords le tlim[1])

    if total(t) lt 1 then continue

    
    ; Phi--
    plim = [ phi[i]-0.5*dp[i], phi[i]+0.5*dp[i] ]
    
    ; keep limits within [-180,180]
    o = where(plim gt 180, no)
    u = where(plim lt -180, nu)
    if no gt 0 then plim[o] += -360
    if nu gt 0 then plim[u] += 360
    
    pr = where( abs(plim - round(plim)) lt tolerance, npr)
    if npr gt 0 then plim[pr] = round(plim[pr])
    
    if plim[0] gt plim[1] then begin
      p = (pcoords gt plim[0]) or (pcoords le plim[1])
    endif else begin
      p = (pcoords gt plim[0]) and (pcoords le plim[1])
    endelse

    if total(p) lt 1 then continue

    
    ; R (velocity)--
    r = (vcoords ge v[i]-0.5*dv[i]) and (vcoords lt v[i]+0.5*dv[i])

;    if total(r) lt 1 then continue


    ; Combine criteria and assign values--
    bidx = where( p and t and r, nb) 

    if nb gt 0 then begin
      weight[bidx]++ 
      part_slice[bidx] += datapoints[i]
    endif  

  endfor
  
  
  ;testing:
;  a = long(unique(weight))
;  print, '------------'
;  for i=0, n_elements(a)-1 do begin
;    b = where(weight eq a[i], nb)
;    print, strtrim(a[i],2)+':  '+strtrim(nb,2)+'  ('+strtrim(100.*nb/n_elements(weight),2)+'%)'
;  endfor
  
  
  ;Average areas where bins overlapped
  adj = where(weight eq 0, na)
  if na gt 0 then weight[adj] = 1b
  part_slice = part_slice / weight
  
  return

end