;+
;Function:
;  spd_slice2d
;
;
;Purpose:
;  Return an interpolated 2D slice of 3D particle data for plotting.
;
;
;Interpolation Methods:
;
;   3D Interpolation (default):
;     The entire 3-dimensional distribution is linearly interpolated onto a 
;     regular 3d grid and a slice is extracted from the volume.
;     
;   2D Interpolation:
;     Datapoints within the specified theta or z-axis range are projected onto 
;     the slice plane and linearly interpolated onto a regular 2D grid.      
;    
;   Geometric:
;     Each point on the plot is given the value of the bin it intersects.
;     This allows bin boundaries to be drawn at high resolutions.
;
;
;Calling Sequence:
;  
;  slice = spd_slice2d( data [,data2 [,data3 [,data4]]] $
;                       [,time=time [,window=window | samples=samples]]
;                       [trange=trange] ... )
;
;
;Example Usage:
;  slice = spd_slice2d(data, time=time)            ;get slice from distribution closest to TIME
;  slice = spd_slice2d(data, time=time, samples=4) ;use average of the 4 samples nearest to TIME
;  slice = spd_slice2d(data, time=time, window=10) ;use average of all data within [TIME,TIME+10sec]
;                                                  ;add "/center_time" for [TIME-5sec,TIME+5sec]  
;  slice = spd_slice2d(data, trange=trange)        ;use average of all data within TRANGE
;
;  slice = spd_slice2d(data, time=time, /three_d_interp) ;use 3D interpolation (default) 
;  slice = spd_slice2d(data, time=time, /two_d_interp)   ;use 2D interpolation
;  slice = spd_slice2d(data, time=time, /geometric)      ;use geometric interpolation
;
;  See crib sheets:  
;    THEMIS:  thm_crib_part_slice2d  (called by thm_part_slice2d)
;    MMS:     mms_slice2d_fpi_crib
;             mms_slice2d_hpca_crib
;    
;
;Arguments:
;  DATARR[#]: An array of pointers to 3D data structures.
;            See spd_dist_array.pro for more.
; 
; 
;Basic Keywords:
;  TRANGE: Two-element time range over which data will be averaged. (string or double)
;  TIME: Time at which the slice will be computed. (string or double)
;    SAMPLES: Number of nearest samples to TIME to average. (int/double)
;             If neither SAMPLES nor WINDOW are specified then default=1.
;    WINDOW: Length in seconds from TIME over which data will be averaged. (int/double)  
;      CENTER_TIME: Flag denoting that TIME should be midpoint for window instead of beginning.
;      
;    SUM_SAMPLES: Flag denoting that the data should be summed over the requested trange rather than averaged
;  
;  THREE_D_INTERP: Flag to use 3D interpolation method (described above)      
;  TWO_D_INTERP: Flag to use 2D interpolation method (described above)
;  GEOMETRIC: Flag to use geometric interpolation method (described above)
; 
;
;Orientation Keywords:
;  CUSTOM_ROTATION: Applies a custom rotation matrix to the data.  Input may be a
;                   3x3 rotation matrix or a tplot variable containing matrices.
;                   If the time window covers multiple matrices they will be averaged.
;                   This is applied before other transformations
;  ROTATION: Aligns the data relative to the magnetic field and/or bulk velocity.
;            This is applied after the CUSTOM_ROTATION. (BV and BE are invariant 
;            between coordinate systems)
;
;            Use MAG_DATA keyword to specify magnetic field vector.
;            Use VEL_DATA keyword to specify bulk velocity (optional).
;       
;     'BV':  The x axis is parallel to B field; the bulk velocity defines the x-y plane
;     'BE':  The x axis is parallel to B field; the B x V(bulk) vector defines the x-y plane
;     'xy':  (default) The x axis is along the data's x axis and y is along the data's y axis
;     'xz':  The x axis is along the data's x axis and y is along the data's z axis
;     'yz':  The x axis is along the data's y axis and y is along the data's z axis
;     'xvel':  The x axis is along the data's x axis; the x-y plane is defined by the bulk velocity 
;     'perp':  The x axis is the bulk velocity projected onto the plane normal to the B field; y is B x V(bulk)
;     'perp_xy':  The data's x & y axes are projected onto the plane normal to the B field
;     'perp_xz':  The data's x & z axes are projected onto the plane normal to the B field
;     'perp_yz':  The data's y & z axes are projected onto the plane normal to the B field
;     
;  SLICE_X & SLICE_NORM: These keywords respectively specify the slice plane's 
;                        x-axis and normal within the coordinates specified by  
;                        CUSTOM_ROTATION and ROTATION. Both keywords take 
;                        3-vectors as input. (See note below)
;                       
;                        If SLICE_X is not specified then the given coordinate's 
;                        x-axis will be used. If SLICE_X is not perpendicular to 
;                        the normal it's projection onto the slice plane will be used.
;                        An error will be thrown if no projection exists.
;                        
;                        If SLICE_NORM is not specified then the given coordinate's
;                        z-axis will be used (slice along by x-y plane in those 
;                        coordinates).
;                 
;       examples:
;         Slice along the data's x-z plane: 
;           ROTATION='xz' 
;         
;         Slice plane's x axis is GSM x and y is in the direction of the bulk velocity:
;           CUSTOM_ROTATION='my_gsm_tvar', ROTATION='xvel'
;
;         Slice is perpendicular to "tvar1" and x axis is defined by projection of "tvar2"
;           SLICE_NORM='tvar1', SLICE_X='tvar2'
;           
;       NOTE: Update at 06/04/2018 - The SLICE_X & SLICE_NORM are defined after CUSTOM_ROTATION
;         but before the ROTATION.          
;           
; DISPLACEMENT: Vector. New center of the coordinate system.
;       example:
;         Slice at the point x=0.5, y = 0.5 and z=0.1.
;         DISPLACEMENT = [0.5, 0.5. 0.1]
;
;  MAG_DATA: Name of tplot variable containing magnetic field data or 3-vector.
;            This will be used for slice plane alignment and must be in the 
;            same coordinates as the particle data.
;  VEL_DATA: Name of tplot variable containing the bulk velocity data or 3-vector.
;            This will be used for slice plane alignment and must be in the
;            same coordinates as the particle data.
;            If not set the bulk velocity will be automatically calculated
;            from the distribution (when needed).
;Other Keywords:
;  RESOLUTION: Integer specifying the resolution along each dimension of the
;              slice (defaults:  2D/3D interpolation: 150, geometric: 500) 
;  SMOOTH: An odd integer >=3 specifying the width of a smoothing window in # 
;          of points.  Smoothing is applied to the final plot using a gaussian 
;          convolution. Even entries will be incremented, 0 and 1 are ignored.
;
;  ENERGY: Flag to plot data against energy (in eV) instead of velocity.
;  LOG: Flag to apply logarithmic scaling to the radial measure (i.e. energy/velocity).
;       (on by default if /ENERGY is set)
;
;  ERANGE: Two element array specifying the energy range to be used in eV.
;
;  THETARANGE: (2D interpolation only)
;              Angle range, in degrees [-90,90], used to calculate slice.
;              Default = [-20,20]; will override ZDIRRANGE. 
;  ZDIRRANGE: (2D interpolation only)
;             Z-Axis range, in km/s, used to calculate slice.
;             Ignored if called with THETARANGE.
;
;  AVERAGE_ANGLE: (geometric interpolation only)
;                 Two element array specifying an angle range over which 
;                 averaging will be applied. The angle is measured 
;                 from the slice plane and about the slice's horizontal axis; 
;                 positive in the right handed direction. This will
;                 average over all data within that range. 
;                 
;                 Note: for the default rotation='xy', the angle is measured from the XY 
;                 slice plane and about the x-axis
;                    e.g. rotation='xy', average_angle=[-25,25] will average data within 25 degrees
;                         of the XY slice plane about it's x-axis
;                    or     
;                         rotation='yz', average_angle=[-25,25] will average data within 25 degrees
;                         of the YZ slice plane about it's y-axis
;
;  SUM_ANGLE: (geometric interpolation only)
;                 Two element array specifying an angle range over which 
;                 summing will be applied. The angle is measured 
;                 from the slice plane and about the slice's horizontal axis; 
;                 positive in the right handed direction. This will
;                 sum over all data within that range. 
;                 
;                 Note: for the default rotation='xy', the angle is measured from the XY 
;                 slice plane and about the x-axis
;                    e.g. rotation='xy', sum_angle=[-25,25] will sum data within 25 degrees
;                         of the XY slice plane about it's x-axis
;                    or     
;                         rotation='yz', sum_angle=[-25,25] will sum data within 25 degrees
;                         of the YZ slice plane about it's y-axis
;               
;  MSG_OBJ: Reference to dprint display object.
;  
;  DETERM_TOLERANCE:  tolerance of the determinant of the custom rotation matrix 
;           (maximum acceptable difference from determ(C)=1 where C is the 
;           user's custom rotation matrix); default is 1e-6
; 
;  SUBTRACT_BULK: Subtract bulk velocity vector
;  PERP_SUBTRACT_BULK: Subtract the component of velocity perpendicular to the B field
;Output:
; Returns a structure to be passed to spd_slice2d_plot:
;      {
;       data: two dimensional array (NxN) containing the data to be plotted
;       xgrid: N dimensional array of x-axis values for plotting 
;       ygrid: N dimensional array of y-axis values for plotting
; 
;       project_name: name of project
;       spacecraft: spacecraft designation
;       data_name: string or string array containing the type(s) of distribution used
;
;       mass: particle mass in ev/(km/s)^2
;       units: the data's units
;       xyunits: the x & y axes' units
;       coord: placeholder for coordinate system label
;       rot: the applied rotation option
;       type: flag denoting interpolation type (0,2,3 for geometric, 2D, 3D respectively)
;       rlog: flag denoting radial log scaling
;
;       zrange: two-element array containing the range of the un-interpolated data 
;       rrange: two-element array containing the radial range of the data
;       trange: two-element array containing the numerical time range
;        
;       bulk: 3-vector containing the bulk velocity in the slice plane's coordinates
;       bfield: 3-vector containing the B-field in the slice plane's coordinates
;       sunvec: 3-vector containing the sun direction in the slice plane's coordinates
;       custom_matrix: The applied custom rotation matrix.
;       rotation_matrix: Rotation matrix from the the original or custom coordinates 
;                        to those defined by ROTATION.
;       orient_matrix: Rotation matrix from the coordinates defined by ROTATION to 
;                      the coordinates defined by SLICE_NORM and SLICE_X 
;                      (column matrix of new coord's basis).
;       }
; 
;
;NOTES:
;   - Regions containing no data are assigned zeros instead of NaNs.
;   - Interpolation may occur across data gaps or areas with recorded zeroes
;     when using 3D interpolation (use geometric interpolation to see bins).
;   - The center/midpoint time of a distribution is used as it's timestamp
;     when determining it's inclusion in the requested time range.  The full
;     time range of all included samples is stored in the metadata.
;
;
;CREATED BY:
;  Aaron Flores, based on work by Bryan Kerr, Arjun Raj, and Xuzhi Zhou
;
;
;$LastChangedBy: egrimes $
;$LastChangedDate: 2022-03-08 13:25:16 -0800 (Tue, 08 Mar 2022) $
;$LastChangedRevision: 30661 $
;$URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/spdsoft/trunk/general/science/spd_slice2d/spd_slice2d.pro $
;-

function spd_slice2d, input1, input2, input3, input4, $
                    ; Time options
                      time=time_in, $
                      window=window, $
                      center_time=center_time, $
                      trange=trange_in, $
                      samples=samples, $
                      sum_samples=sum_samples, $ ; sum the data rather than average
                    ; Orientations
                      custom_rotation=custom_rotation, $
                      determ_tolerance=determ_tolerance, $
                      rotation=rotation, $
                      slice_norm=slice_z, $
                      slice_x=slice_x, $
                     ; Shift the origin before interpolation
                      displacement=displacement_in, $                        
                    ; Support Data
                      mag_data=mag_data, $
                      vel_data=vel_data, $
                      sun_data=sun_data, $
                    ; Interpolation options
                      geometric=geometric, $
                      two_d_interp=two_d_interp, $
                      three_d_interp=three_d_interp, $
                    ; 2D interpolation options
                      thetarange=thetarange, $
                      zdirrange=zdirrange, $
                    ; Range limits
                      erange=erange, $
                    ; Other options
                      resolution=resolution, $
                      smooth=smooth,  $
                      log=log, $
                      energy=energy, $
                    ; TBD
                      subtract_bulk=subtract_bulk, $
                      average_angle=average_angle, $
                      sum_angle=sum_angle, $
                    ; Other
                      msg_obj=msg_obj, $
                      fail=fail, $
                    ; Test
                      perp_subtract_bulk=perp_subtract_bulk, $ ;only subtract velocity perp to B
                      perp_vel_data=perp_vel_data, $           ;tplot variable with perp velocity
                      _extra = _extra

    compile_opt idl2


invalid = 0b
fail = ''


if ~ptr_valid(input1[0]) and ~is_struct(input1[0]) then begin
  fail = 'Invalid data.  Input must be pointer or structure array.'
  dprint, dlevel=1, fail
  return, invalid
endif

if undefined(trange_in) then begin
  if undefined(time_in) then begin
    fail = 'Please specify a time or time range over which to compute the slice.  For example: '+ $
           ssl_newline()+'  "TIME=t, WINDOW=w" or "TRANGE=tr" or "TIME=t, SAMPLES=n"'
    dprint, dlevel=1, fail
    return, invalid
  endif else begin
    if undefined(window) && undefined(samples) then begin
      samples = 1 ;use single closest distribution by default
    endif
  endelse
endif

valid_rotations = ['bv', 'be', 'xy', 'xz', 'yz', 'xvel', $
                   'perp', 'perp2', 'perp_xy', 'perp_xz', 'perp_yz', 'b_exb', 'perp1-perp2']
if ~undefined(rotation) then begin
  if ~in_set(strlowcase(rotation),valid_rotations) then begin
    fail = 'Invalid rotation requested.  See spd_crib_slice2d for examples.'
    dprint, dlevel=1, fail
    return, invalid
  endif
endif



; Set default options
;------------------------------------------------------------
if ~undefined(time_in) then time = time_double(time_in[0])
if undefined(rotation) then rotation='xy'
if keyword_set(energy) && undefined(log) then log = 1
if rotation eq 'xyz' then rotation = 'xy'
if rotation eq 'perp_xyz' then rotation = 'perp_xy'

;Interpolation type:
if keyword_set(geometric) then type=0
if keyword_set(two_d_interp) then type = 2
if keyword_set(three_d_interp) then type = 3
if undefined(type) then type = 3


; 2d interp
if type eq 2 then begin 
  if undefined(smooth) then smooth = 7
  if undefined(resolution) then resolution = 150
  if ~keyword_set(thetarange) and ~keyword_set(zdirrange) then begin
    thetarange = [-20.,20]
  endif
endif 

; 3D interp
if type eq 3 then begin 
  if undefined(smooth) then smooth = 7
  if undefined(resolution) then resolution = 150
endif 

; geometric
if ~keyword_set(type) then begin
  if undefined(resolution) then resolution = 500
endif

if type ne 0 && keyword_set(average_angle) then begin
  dprint, dlevel=1, 'Angular averaging is only applicable to the geometric method.'+ $
                    'No averaging will be applied.'
  return, invalid
endif

if type ne 0 && keyword_set(sum_angle) then begin
  dprint, dlevel=1, 'Angular summing is only applicable to the geometric method.'+ $
    'No summing will be applied.'
  return, invalid
endif

if keyword_set(average_angle) && keyword_set(sum_angle) then begin
  dprint, dlevel=1, 'Conflicting keywords specified: sum_angle and average_angle; only one of these keywords can be specified at a time.' + $
    'No summing/averaging will be applied.'
  return, invalid
endif

if undefined(time) then tslice = trange_in[0] else tslice = time
msg = 'Processing slice at ' + time_string(tslice, tformat='YYYY/MM/DD hh:mm:ss.fff') +'... '
dprint, dlevel=2, msg, display_object=msg_obj 



; Aggregate inputs
;  -multiple inputs are allowed to make things easier on the user
;  -pointer arrays are used as a way of supporting dissimilar 
;   structures, which cannot be concatenated
;------------------------------------------------------------

;copy any structs to new pointer var, other inputs are copied and checked next
switch n_params() of 
  4: if is_struct(input4) then p4 = ptr_new(input4) else p4 = input4 
  3: if is_struct(input3) then p3 = ptr_new(input3) else p3 = input3
  2: if is_struct(input2) then p2 = ptr_new(input2) else p2 = input2
  1: if is_struct(input1) then p1 = ptr_new(input1) else p1 = input1
endswitch

;agregate valid pointers
ds = p1
if ptr_valid(p2) then ds = [ds,p2]
if ptr_valid(p3) then ds = [ds,p3]
if ptr_valid(p4) then ds = [ds,p4]



; Get the slice's time range
;------------------------------------------------------------

; get the time range if one was specified
if ~undefined(trange_in) then begin
  trange = minmax(time_double(trange_in))
endif

; get the time range if a time & window were specified instead
if undefined(trange) && keyword_set(window) then begin
  if keyword_set(center_time) then begin
    trange = [time - window/2d, $
              time + window/2d  ]
  endif else begin
    trange = [time, $
              time + window ]
  endelse
endif

; if no time range or window was specified then get a time range 
; from the N closest samples to the specified time
;   (defaults to 1 if SAMPLES is not defined)
if undefined(trange) then begin
  trange = spd_slice2d_nearest(ds, time, samples)
endif
  

; check that there is data in range before proceeding
for i=0, n_elements(ds)-1 do begin
  times_ind = spd_slice2d_intrange(ds[i], trange, n=ndat)
  n_samples = array_concat(ndat,n_samples)
endfor
n_samples = total(n_samples)
if n_samples lt 1 then begin
  fail = 'No particle data in the time window: '+strjoin(time_string(trange),', ')+ $
         '. Time samples may be at low cadence; try adjusting the time window.'
  dprint, dlevel=1, fail 
  return, invalid
endif
dprint, dlevel=3, strtrim(n_samples,2) + ' samples in time window'


; Extract particle data from structures
;  -apply energy limits
;  -average data over time window
;  -output r, phi, theta and dr, dphi, dtheta arrays
;------------------------------------------------------------
spd_slice2d_get_data, ds, trange=trange, erange=erange, energy=energy, fail=fail, $ 
                 data=datapoints, rad=rad, phi=phi, theta=theta, dr=dr, dp=dp, dt=dt, $
                 sum_samples=sum_samples
if keyword_set(fail) then return, invalid


; Known causes of an empty data variable should be caught before this.
if ~keyword_set(datapoints) then begin
  fail = 'Unknown error extracting data from particle distributions.'
  dprint, dlevel=0, fail
  return, invalid
endif


; Get original data and radial ranges for plotting
;  -ignore outliers that may be the result of sanitizations performed outside this routine
;  -this also effectively removes low-significance outliers at large velocities, which
;   can lead to overly large velocity scales for some instruments 
idx = where(datapoints gt 0,n)
if n gt 0 then begin
  dmoms = moment(alog10(datapoints[idx]),maxmom=2)
  min = 10^(dmoms[0] - 2*sqrt(dmoms[1])) ;ignore if < mean - 2*sigma 
  drange = minmax(datapoints[idx],min_value=min)

;*** affects interpolation in some cases ***
;  ;remove points below the minimum 
;  ;  -these points would be hidden on final plot if left and will only 
;  ;   slow the interpolation
;  idx = where(datapoints gt 0, n)
;  if n gt 0 then begin
;    datapoints = datapoints[idx]
;    rad = rad[idx]
;    phi = phi[idx]
;    theta = theta[idx]
;    dr = dr[idx]
;    dp = dp[idx]
;    dt = dt[idx]
;  endif
;
;  drange = [min(datapoints),max(datapoints)]
endif else begin
  drange = [0,0.]
endelse
rrange = [  min( rad - 0.5*dr ), max( rad + 0.5*dr )  ]


; Apply radial log scaling
if keyword_set(log) then begin
  spd_slice2d_rlog, rad, dr
endif


; Convert spherical data to cartesian coordinates for interpolation
if keyword_set(type) then begin
  spd_slice2d_s2c, rad,theta,phi, xyz
endif



; Get/apply coordinate transformations
;  -for 2D/3D interpolation the data is transformed here
;  -for geometric interpolation the data is transformed in spd_slice2d_geo
;  -support data is always rotated at each step
;------------------------------------------------------------

; Get support data for aligning slice
spd_slice2d_get_support, mag_data, trange, output=bfield
spd_slice2d_get_support, vel_data, trange, output=vbulk
spd_slice2d_get_support, sun_data, trange, output=sunvec
spd_slice2d_get_support, slice_x, trange, output=slice_x_vec
spd_slice2d_get_support, slice_z, trange, output=slice_z_vec

; Create perp velocity tplot variable, carefully, by rotating into
; 'bv' system, then set X component of velocity ot zero, and inverse
If(keyword_set(perp_subtract_bulk)) Then Begin
   rotation_tmp = 'bv'
   vbulk_tmp = vbulk
   bfield_tmp = bfield
   spd_slice2d_rotate, rotation=rotation_tmp, fail=fail, $ 
                       vbulk=vbulk_tmp, bfield=bfield_tmp, $
                       matrix=rotation_matrix_tmp
   if keyword_set(fail) then return, invalid
   perp_vbulk = vbulk_tmp & perp_vbulk[0] = 0.00
;now unrotate the perp_vbulk vector
   perp_vbulk = transpose(rotation_matrix_tmp) ## perp_vbulk
Endif

; Custom rotation
spd_slice2d_custom_rotation, custom_rotation=custom_rotation, trange=trange, fail=fail, $
         vectors=xyz, bfield=bfield, vbulk=vbulk, sunvec=sunvec, matrix=custom_matrix, $
         determ_tolerance=determ_tolerance, perp_vbulk=perp_vbulk
if keyword_set(fail) then return, invalid

; Rotation defined by slice normal/x axis options

if (~undefined(slice_x_vec) && ~ARRAY_EQUAL([1,0,0],slice_x_vec)) || $
   (~undefined(slice_z_vec) && ~ARRAY_EQUAL([0,0,1],slice_z_vec)) then begin
   dprint, dlevel=2, 'Warning! behaviour of slice_x_vec and slice_z_vec has been changed (06/04/2018). See the description of the function.', display_object=msg_obj
endif

spd_slice2d_orientslice, slice_x=slice_x_vec, slice_z=slice_z_vec, fail=fail, $
  vectors=xyz, bfield=bfield, vbulk=vbulk, sunvec=sunvec, matrix=orient_matrix, $
  perp_vbulk=perp_vbulk
if keyword_set(fail) then return, invalid

; Built in rotation option
spd_slice2d_rotate, rotation=rotation, fail=fail, $ 
         vectors=xyz, bfield=bfield, vbulk=vbulk, sunvec=sunvec, matrix=rotation_matrix, $
         perp_vbulk=perp_vbulk
if keyword_set(fail) then return, invalid


; Subtract bulk velocity vector
if keyword_set(subtract_bulk) && ~keyword_set(log) then begin
  spd_slice2d_subtract, vectors=xyz, velocity=vbulk, fail=fail
  if keyword_set(fail) then return, invalid
  ;as with rotations, this will be applied later for geometric interp
  geo_shift = vbulk
endif else if keyword_set(perp_subtract_bulk) && ~keyword_set(log) then begin
  spd_slice2d_subtract, vectors=xyz, velocity=perp_vbulk, fail=fail
  if keyword_set(fail) then return, invalid
  ;as with rotations, this will be applied later for geometric interp
  geo_shift = perp_vbulk
endif

; Misc.
;------------------------------------------------------------

; Sort transformed vector grid
;  -this ends up being necessary later 
if ~undefined(xyz) then begin
  sorted = sort(xyz[*,0])
  xyz = xyz[sorted,*]
  datapoints = datapoints[sorted]
endif

minvelinterp = 0
if keyword_set(displacement_in) then begin
  ; simple preventetive measure the code will still work if the inpute is wrong
  if n_elements(displacement_in) ne 3 then begin
    disp_arr = [0., 0., 0.]
    dprint, dlevel=2, 'Displacement must be an array of 3 elements. Displacement is set to zero.', display_object=msg_obj
  endif else disp_arr = displacement_in
   
  d = disp_arr
  case strlowcase(rotation) of
    'xy': disp_arr = [d[0], d[1], d[2]]
    'xz': disp_arr = [d[0], d[2], -1*d[1]] ; -1 because of the rotation see spd_cal_rot([1,0,0],[0,0,1])
    'yz': disp_arr = [d[1], d[2], d[0]]
    else: begin
      disp_arr = [0., 0., 0.]
      dprint, dlevel=2, 'Displacement currently supports only rotation=xy, xz, yz. Displacement is set to zero.', display_object=msg_obj      
    end
  endcase  
  
  if type eq 2 and keyword_set(thetarange) then begin
    dprint, dlevel=2, 'Use zdirrange instead of thetarange with displacement in 2d interpolation.', display_object=msg_obj
  endif
  
  if type eq 1 then begin 
    dprint, dlevel=2, 'Displacement supports 2d and 3d interpolation. Displacement is canceled.', display_object=msg_obj    
  endif else begin 
    spd_slice2d_translate, vectors=xyz, translate=-1*disp_arr, data=datapoints, /truncate
    minvelinterp = 1
  endelse
  
endif


; Create slice:
; 
; TYPE=0   Geometric Interpolation
; TYPE=2   2D Linear Interpolation
; TYPE=3   3D Linear Interpolation
;------------------------------------------------------------

; Linear 2D interpolation from thm_esa_slice2D
if type eq 2 then begin

  dprint, dlevel=4, 'Using 2d linear interpolation'
  spd_slice2d_2di, datapoints, xyz, resolution, $
                   thetarange=thetarange, zdirrange=zdirrange, $
                   slice=slice, xgrid=xgrid, ygrid=ygrid, $
                   minvelinterp=minvelinterp, $
                   fail=fail

; Linear 3D Interpolation
endif else if type eq 3 then begin

  dprint, dlevel=4, 'Using 3d linear interpolation'
  spd_slice2d_3di, datapoints, xyz, resolution, drange=drange, $ 
                   slice=slice, xgrid=xgrid, ygrid=ygrid, $
                   minvelinterp=minvelinterp, $
                   fail=fail

; Geometric Method
endif else begin

  dprint, dlevel=4, 'Using geometric method'
  spd_slice2d_geo, data=datapoints, resolution=resolution, $
                   rad=rad, phi=phi, theta=theta, dr=dr, dp=dp, dt=dt, $
                   custom_matrix=custom_matrix, $
                   rotation_matrix=rotation_matrix, $
                   orient_matrix=orient_matrix, $
                   shift = geo_shift, $
                   average_angle=average_angle, $
                   msg_obj=msg_obj, msg_prefix=msg, $
                   slice=slice, xgrid=xgrid, ygrid=ygrid, $
                   sum_angle=sum_angle, $
                   fail=fail

endelse

if keyword_set(fail) then begin
  dprint, dlevel=1, fail
  return, invalid
endif


; Apply smoothing
if keyword_set(smooth) then begin
  spd_slice2d_smooth, slice, smooth
endif

; Change the axes
if keyword_set(displacement_in) then begin
  if type ne 1 then begin 
    xgrid=xgrid + disp_arr[0]
    ygrid=ygrid + disp_arr[1]
  endif
endif


; Get metadata and return slice structure
;------------------------------------------------------------

; loop over dist structures to get name(s) and time range
for i=0, n_elements(ds)-1 do begin
  
  ;dist type name
  data_name = array_concat( ((*ds[i]).data_name)[0], data_name )

  ;time range of particle data
  times_ind = spd_slice2d_intrange(ds[i], trange, n=ndat)
  if ndat gt 0 then begin
    trc = [ min((*ds[i])[times_ind].time), $
            max((*ds[i])[times_ind].end_time) ]
    tr = minmax( array_concat(trc,tr) ) 
  endif

endfor

; don't repeat distribution types in name
data_name = strjoin( spd_uniq(data_name), '/')


; Output structure
slice = {  $
         ; Data
          data: temporary(slice), $  ;data values
          xgrid: xgrid, $    ;x-axis values
          ygrid: ygrid, $    ;y-axis values

         ; Metadata
          project_name: ((*ds[0]).project_name)[0], $
          spacecraft: ((*ds[0]).spacecraft)[0], $
          data_name: data_name, $            ;names of distributions used (string)
          n_samples: n_samples, $            ;number of distributions used

          mass: ((*ds[0]).mass)[0], $    ;mass in eV/(km/s)^2
          units: ((*ds[0]).units_name)[0], $   ;data units
          xyunits: keyword_set(energy) ? 'eV':'km/s', $ ;x/y/r units (string)
          coord: '', $       ;placeholder for coordinates used (string)
          rot: rotation, $   ;rotation used (string)
          type: type, $      ;integer specifying slice type (0,2,3)
          energy: keyword_set(energy), $ ;flag for energy plots
          rlog: keyword_set(log), $  ;flag for radial log plot

          trange: tr, $      ;total time range of used distributions
          zrange: drange, $  ;data range after averaging but before interpolation
          rrange: rrange , $ ;instrument velocity/energy range

         ; Support Data
          shift: undefined(subtract_bulk) ? 0:vbulk, $ ;for plotting energy limits
          bulk: undefined(vbulk) ? 0:vbulk, $      ;bulk velocity vector
          perp_bulk: undefined(perp_vbulk) ? 0:perp_vbulk, $      ;perp bulk velocity vector
          bfield: undefined(bfield) ? 0:bfield, $ ;bfield vector
          sunvec: undefined(sunvec) ? 0:sunvec, $  ;sun vector
          custom_matrix: custom_matrix, $ ;applied custom rotation matrix
          rotation_matrix: rotation_matrix, $ ;rotation matrix from ROTATION option 
          orient_matrix: orient_matrix $ ;rotation matrix specified by SLICE_NORM/SLICE_X 
          }


msg = 'Finished slice at '+time_string(tslice, format=5)
dprint, dlevel=2, msg, display_object=msg_obj


return, slice

end
