;+
; 
; Name: thm_part_slice2d_2di.pro
; 
; Purpose:  Helper function for thm_part_slice2d.pro
;           Produces slice using 2D linear interpolation
;           
; Note: This code is meant to preserve the functionality of
;       thm_esa_slice2d for the new set of slices routines.
;
;-
pro thm_part_slice2d_2di, datapoints, velocities, resolution, $
                          thetarange=thetarange, zdirrange=zdirrange, $ 
                          part_slice=part_slice, $
                          xgrid=xgrid, ygrid=ygrid, $
                          fail=fail

    compile_opt idl2, hidden


  ;catch common "points colinear" error from triangulate
  catch, err
  if err ne 0 then begin
    catch, /cancel
    if strmessage(err,/name) eq 'IDL_M_POINTS_COLINEAR' then begin
      fail = strjoin(['WARNING: 2D Interpolation cannot be used for this data. ', $
             ' 3D Interpolation is recommended if the data is free of large gaps in look direction and energy.', $
             ' Otherwise use the Geometric method with low resolution (50-150) and smoothing (3-5pts). '],string(10b))
    endif else begin
      message, /reissue_last
    endelse
    return
  endif


  ; Cut by theta value
  if keyword_set(thetarange) then begin
  
    thetarange = minmax(thetarange)
    
    r = sqrt(velocities[*,0]^2 + velocities[*,1]^2 + velocities[*,2]^2)
    eachangle = asin(velocities[*,2]/r)/!dtor
    
    index = where(eachangle le thetarange[1] and eachangle ge thetarange[0], $
                  count, ncomplement=ncount)
    if count ne 0 then begin
      if ncount ne 0 then begin
        velocities = velocities[index,*]
        datapoints = datapoints[index]
      endif
    endif else begin
      fail = 'No data points in given theta range'
      return
    endelse
    
  ;Cut by z-axis value
  endif else if keyword_set(zdirrange) then begin
    
    zdirrange = minmax(zdirrange)
    
    index = where(velocities[*,2] ge zdirrange[0] and velocities[*,2] le zdirrange[1], $
                  count, ncomplement=ncount)
    if count ne 0 then begin
      if ncount ne 0 then begin
        velocities = velocities[index,*]
        datapoints = datapoints[index]
      endif
    endif else begin
      fail = 'No data points in given Z direction range'
      return
    endelse  
  endif

  ;copy vars, ignoring z comp
  xvels = velocities[*,0]
  yvels = temporary(velocities[*,1]) ;save mem

  ; Average duplicate points (the long way)---------
  uni2=uniq(xvels)
  uni1=[0,uni2[0:n_elements(uni2)-2]+1]
  
  kk=0
  for i=0,n_elements(uni2)-1 do begin
      yvelsi=yvels[uni1[i]:uni2[i]]
      xvelsi=xvels[uni1[i]:uni2[i]]
      datapointsi=datapoints[uni1[i]:uni2[i]]
  
      xvelsi=xvelsi[sort(yvelsi)]
      datapointsi=datapointsi[sort(yvelsi)]
      yvelsi=yvelsi[sort(yvelsi)]
  
    index2=uniq(yvelsi)
    if n_elements(index2) eq 1 then begin
        index1=0
    endif else begin
        index1=[0,index2[0:n_elements(index2)-2]+1]
    endelse
  
      for j=0,n_elements(index2)-1 do begin
          yvels[kk]=yvelsi[index1[j]]
          xvels[kk]=xvelsi[index1[j]]
          if index1[j] eq index2[j] then begin
              datapoints[kk]=datapointsi[index1[j]]
          endif else begin
              datapoints_moment=moment(datapointsi[index1[j]:index2[j]])
              datapoints[kk]=datapoints_moment[0]
          endelse
          kk=kk+1
      endfor
  endfor

  yvels=yvels[0:kk-1]
  xvels=xvels[0:kk-1]
  datapoints=datapoints[0:kk-1]
  ; ---------------


  ; Create triangulation
  triangulate,xvels,yvels,tr,b
  
  ; Remove triangles whose total x-y plane velocity is less than
  ; minimum velocity from distribution (prevents interpolation over
  ; lower energy limits)
  index = where( 1./9 * total( xvels[tr[0:2,*]] ,1 )^2 + $
                 1./9 * total( yvels[tr[0:2,*]] ,1 )^2 $
                  gt min(xvels^2+yvels^2), $
                  count, ncomplement=ncomp)
  if count gt 0 then begin
    if ncomp gt 0 then begin
      tr=tr[*,index]
    endif
  endif else dprint,dlevel=4, 'Possible error in triangulation' 
  
  ; Set spacing
  xmax = max(abs([yvels,xvels]))
  xrange = [-1*xmax,xmax]
  spacing = (xrange[1]-xrange[0])/(resolution-1)
  
  ; Interpolate to regular grid (slice data)
  part_slice = trigrid(xvels,yvels,datapoints,tr,[spacing,spacing], $
                       [xrange[0],xrange[0],xrange[1],xrange[1]],   $
                       xgrid = xgrid,ygrid = ygrid )


end