;+
; 
; 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.
;       The code is essentially copied from that file with 
;       some minor modifications.
;
;-
pro thm_part_slice2d_2di, datapoints, velocities, resolution, $
                            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 err eq -471 then begin
      fail = strjoin(['WARNING: 2D Interpolation cannot be used for this data. ', $
             'Geometric method with low resolution (50-100) and smoothing (3pts) is recommended instead. ', $
             '3D Interpolation will work if the data is free of gaps in look direction and energy.'],string(10b))
    endif else begin
      message, /reissue_last
    endelse
    return
  endif


  xvels = velocities[*,0]
  yvels = temporary(velocities[*,1])

  ; 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, '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