; Returns vector with mag data or -1 on error
;
function thm_part_slice2d_getmag, thedata, datstructs, magdata=magdata, $
rotation=rotation, fail=fail
compile_opt idl2, hidden
if keyword_set(magdata) then begin
; average mag data over time window
bfield = dat_avg(magdata, datStructs[0].time, thedata.end_time)
endif else begin
; average mag data over time window
bfield = average(datStructs.magf, 2)
; check for NaNs that could corrupt rotation matrix
if in_set(finite(bfield),0) then begin
if ~in_set(['xy','xz','yz','xvel'], rotation) then begin
fail = 'Magnetic field data incomplete for desired time range.'
dprint, fail
return, -1
endif
endif
endelse
return, bfield
end
;Returns bulk velocity vector or -1 on error
;
;Note: thm_esa_slice2d Divides total flux (over entire time range) by
; the total density; here the flux/density was previously
; calculated for each distribution and is now averaged.
function thm_part_slice2d_getvel, thedata, datstructs, veldata=veldata, $
rotation=rotation, fail=fail
compile_opt idl2, hidden
if keyword_set(veldata) then begin
;user specified velocitys can currently be added to the
;dat structure in thm_part_dist_array
vbulk = dat_avg(veldata, datStructs[0].time, thedata.end_time)
endif else begin
vbulk = average(datStructs.velocity, 2) ;vel in km/s
; Check for NaNs that could corrupt rotation matrix
if in_set(finite(vbulk), 0) then begin
if in_set(['BV','BE','xvel','perp'], rotation) then begin
fail='Bulk velocity data incomplete for desired time range.'+ $
' Unable to perform corrdinate rotation to '+rotation+'.'
dprint, fail
return, -1
endif
endif
endelse
; convert to m/s
return, vbulk * 1000.
end
;Performs rotations from DSL
; *Assumes the transformation from the data's native DSL to
; the requested coords (GSM, GSE) is invariant over the
; time range of the slice.*
pro thm_part_slice2d_cotrans, vectors, coord, bfield, vbulk, type, $
probe=probe, twin=twin, ct=ct, $
fail=fail
compile_opt idl2, hidden
ctn = '2dslice_temp_cotrans'
if stregex(coord, '^dsl$', /bool, /fold_case) then return
if ~keyword_set(probe) then begin
fail = 'Cannot determine spacecraft. Coordinate transform canceled.'
dprint, fail
return
endif
if ~stregex(coord, '(^gsm$)|(^gse$)', /bool, /fold_case) then begin
fail = 'Invalid coordinate keyword: '+coord
dprint, fail
return
endif
probe = strlowcase(probe)
coord = strlowcase(coord)
dprint, dlevel=4, 'Rotating to '+strupcase(coord)+' coordinates'
;Load necessary support data (2 min padding at each end)
thm_load_state, probe=probe, trange = twin + 120*[-1,1], suffix='_'+ctn, $
/get_support_data, /keep_spin_data, verbose=0
;Transform x,y,z unit vectors to new coordinates
store_data, ctn, data = {x:replicate(mean(twin),3), $
y: [ [1.,0,0], [0,1,0], [0,0,1] ] }, verbose=0
thm_cotrans, ctn, probe=probe, in_coord='dsl', out_coord=coord, $
support_suffix='_'+ctn, out_suff='_'+coord
;Form transformation matrix from transformed basis
get_data, ctn+'_'+coord, data=ctd
if size(ctd,/type) ne 8 then begin
fail = 'Could not obtain coordinate transform from DSL -> '+strupcase(coord)+ $
', STATE data may be absent for the requested period.'
dprint, fail
return
endif
ct = transpose(ctd.y)
;Apply to velocity data and required rotation vectors
if keyword_set(type) then vectors = ct ## temporary(vectors)
if keyword_set(vbulk) then vbulk = ct ## vbulk
if keyword_set(bfield) then bfield = ct ## bfield
;Delete temporary tplot data
store_data, '*'+ctn+'*', /delete, verbose=0
end
;Perform rotation on cartesian velocity bins
pro thm_part_slice2d_rotate, velocities, rotation, bfield, vbulk, type, $
rot=rot, fail=fail
compile_opt idl2, hidden
; Create rotation matrix
case strlowcase(rotation) of
'bv': rot=cal_rot(bfield,vbulk)
'be': rot=cal_rot(bfield,crossp(bfield,vbulk))
'xy': rot=cal_rot([1,0,0],[0,1,0])
'xz': rot=cal_rot([1,0,0],[0,0,1])
'yz': rot=cal_rot([0,1,0],[0,0,1])
'xvel': rot=cal_rot([1,0,0],vbulk)
'perp': rot=cal_rot(crossp(crossp(bfield,vbulk),bfield),crossp(bfield,vbulk))
'perp_yz': rot=cal_rot(crossp(crossp(bfield,[0,1,0]),bfield),crossp(crossp(bfield,[0,0,1]),bfield))
'perp_xy': rot=cal_rot(crossp(crossp(bfield,[1,0,0]),bfield),crossp(crossp(bfield,[0,1,0]),bfield))
'perp_xz': rot=cal_rot(crossp(crossp(bfield,[1,0,0]),bfield),crossp(crossp(bfield,[0,0,1]),bfield))
else: begin
fail = 'Unrecognized rotation: "'+rotation+'".'
dprint, fail
return
end
endcase
; Prevent data from being mutated to doubles
rot = float(rot)
if rotation ne 'xy' then begin
dprint, dlevel=4, 'Alligning slice plane to ' +rotation
endif else return
; Rotate data
if keyword_set(type) then velocities = rot ## temporary(velocities)
if keyword_set(vbulk) then vbulk = rot ## vbulk
if keyword_set(bfield) then bfield = rot ## bfield
end
; Produces the slice's normal and x-vector
;
pro thm_part_slice2d_orientslice, slice_orient, fail=fail, $
matrix = matrix, $
displacement = displacement, $
norm = unit_norm, $
xvec = xvec
compile_opt idl2, hidden
; Create unit vector for slice orientation
slice_norm = float(slice_orient[0:2])
displacement = float(slice_orient[3])
unit_norm = slice_norm / sqrt(total(slice_norm^2)) ; convert slice normal to a unit vec
if in_set(finite(unit_norm), 0) then begin
fail = 'Invalid orientation vector.'
return
endif
; Project current x axis into slice plane or use y axis if orthogonal
a = crossp(unit_norm,crossp([1,0,0], unit_norm))
if total(abs(a)) eq 0. then begin
if unit_norm[0] lt 0 then xvec=[0, -1, 0] else xvec=[0, 1, 0]
endif else begin
xvec = a / sqrt(total(a^2))
endelse
; matrix = cal_rot(xvec, crossp(unit_norm,xvec))
matrix = thm_part_slice_trans(unit_norm, xvec, fail=fail)
end
; Subtracts bulk velocity vector from velocities array
; todo: fix for geo
pro thm_part_slice2d_subtract, velocities, vbulk, veldata=veldata
compile_opt idl2, hidden
if keyword_set(veldata) then begin
dprint, dlevel=4,'Velocity used for subtraction is '+veldata
endif else begin
dprint, dlevel=4, 'Velocity used for subtraction is V_3D'
endelse
velocities[*,0] = velocities[*,0] - vbulk[0]
velocities[*,1] = velocities[*,1] - vbulk[1]
velocities[*,2] = velocities[*,2] - vbulk[2]
end
; Set slicing limits for 2d interpolation method (from thm_esa_slice2d)
;
pro thm_part_slice2d_2dlimits, velocities, datapoints, $
thetarange=thetarange, zdirrange=zdirrange, $
fail = fail
compile_opt idl2, hidden
; 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'
dprint, fail
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'
dprint, fail
return
endelse
endif
end
;+
;Procedure: thm_part_slice2d
;
;Purpose: Returns a 2-D slice of the 3-D THEMIS ESA/SST ion or electron distribution
; function. The 2-D slice is returned via the PART_SLICE, XGRID, YGRID, and
; SLICE_INFO keywords.
;
; This procedure works in conjunction with thm_part_dist_array.pro and
; plot_part_slice2d.pro. The corresponding crib is thm_crib_part_slice2d.pro
;
; There are currently four methods of generating slices for plotting:
;
; Geomtric:
; Exact method showing the bounds of all bins intersecting the slice plane.
; Values are averaged in cases of overlapping bins or when the slice plane
; is coplanar with a bin's boundary.
;
; 2D Interpolation (from thm_esa_slice2d):
; Datapoints within the specified theta or z-axis range are projected onto
; the slice plane and linearly interpolated onto a regular 2d grid.
; X,Y-range values determined slightly differently for 2D Interpolation
;
; 3D Interpolation:
; The entire 3-dimensional distribution is linearly interpolated onto a
; regular 3d grid and a this slice is extracted. The slice will consist of
; up to resolution^2 datapoints from within 0.866/resolution of the total
; velocity range along the normal. This method works best with dense data.
; Note: - Not all points within the aforemention range will be used.
; - Interpolation may occur across data gaps or areas with recorded zeroes
; - Higher resolution regridding will severely slow this method
;
; 2D Nearest Neighbor (mainly for testing/not reccomended):
; A slice of data with user defined width is projected onto the slice plane
; and interpolated onto a regular grid using the nearest neighbor method.
;
;
;
;Callins Sequence:
; thm_part_slice2d, datArr, [datArr2, [datArr3, [datArr4]]], $
; timewin = timewin, slice_time=slice_time, $
; part_slice=part_slice, xgrid=xgrid, ygrid=ygrid, $
; slice_info=slice_info, $
; fail=fail
;
;
;Arguments:
; DATARR[#]: An array of data structures as returned by one of the get_th?_p???
; routines (or thm_part_dist_array.pro). The structure must have
; a magf field containing a three-component magnetic field vector.
;
;Input Keywords:
; SLICE_TIME: Beginning of time window in seconds since Jan. 1, 1970. If
; CENTER_TIME keyword set, then TIME is the center of the time widow
; specified by the TIMEWIN keyword.
; CENTER_TIME: Flag that, when set, centers the time window around the time
; specified by SLICE_TIME keyword.
; TIMEWIN: Length in seconds over which to compute the slice.
; TWO_D_INTERP: Flag to use 2D interpolation method (described above)
; THREE_D_INTERP: Flag to use 3D interpolation method (described above)
; COORD: A string designating what coordinate system the slice will be based off.
; Default is 'DSL', but can also be 'GSE' or 'GSM'
; SLICE_ORIENT: (ignored for 2D interpolation)
; [x, y, z, vel] A four-element array that specifies the
; orientation of the slice plane. The first three elements
; specify the direction of the plane's normal vector. The fourth
; element is the distance of the slice plane's center from the
; origin along the plane's normal vector in units of meters/sec.
; Default = [0,0,1,0]
; UNITS: A string designating the desired units ('Counts','DF','Rate','CRate',
; 'Flux','EFlux','E2Flux','E3Flux')
; Default = DF
; ONECOUNT: Flag to remove bins that fall below the once count average after averaging.
; RESOLUTION: A single integer specfying the resolution of each dimension of the
; slice (as well as the interpolated volume in the the 3d case).
; Default = 250
; REGRID: A three element array specifying regrid dimensions in phi, theta, and
; energy respectively. If set, all distributions' data will first be
; spherically interpolated using the nearest neighbor method. For example
; the default [64,32,64] will yield distributions with 63x32 points per
; energy, and 64 energy bins). Increasing the regridding resolution will
; slow slice production, particularly if using 3D interpolation.
; Default = [64,32,64]
; ROTATION: The rotation keyword can specify the orienation of the slice plane
; within the specified coordinates.
; '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'/'xyz': (Default) The x axis is V_x and the y axis is V_y.
; 'xz': The x axis is V_x and the y axis is V_z.
; 'yz': The x axis is V_y and the y axis is V_z.
; 'xvel': The x axis is V_x; the y axis 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'/'perp_xyz': Geometric xy coordinates projected onto the plane normal to the B field.
; 'perp_xz': Geometric xz coordinates projected onto the plane normal to the B field.
; 'perp_yz': Geometric yz coordinates projected onto the plane normal to the B field.
; ERANGE: Two element array specifying the energy range to be used
; SUBTRACT: (buggy, check again later)
; Flag to subtract the bulk velocity before plot. The exact quantity
; to use can be specified in the call to thm_part_dist_array.
; SLICE_WIDTH: (2d nearest neighbor only)
; Specifies width of the slice in % of the total range
; Default =
; SMOOTH: An odd integer >=3 specifing the width of the smoothing window in #
; of points. Even entries will be incrimented, 0 and 1 are ignored.
; 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.
;
;
;Output variables for plotting:
; PART_SLICE: The 2-D data array to be plotted (the slice).
; XGRID: Array of x-locations for the slice.
; YGRID: Array of y-locations for the slice.
; SLICE_INFO: Structure to be passed to plot_part_slice2d to determine plot labels.
; {
; probe: string containing the probe
; dist: string (or array of) containing the type of distribution used
; coord: string describing the coordinate system used for the slice
; rot: string describing the user specified rotation (N/A for 2D interp)
; units: string describing the units
; twin: string describing the time window of the slice
; range: array containing the range of the un-interpolated data
; vrange: array containing the velocity range of the instrument(s)
; trange: array containing the numerical time range
; shift: 3-vector containing any translations made in addition to
; requested rotations (e.g. subtracted bulk velocity)
; bulk: 3-vector containing the bulk velocity in the slice plane's
; coordinates
; coord_m: Rotation matrix from original data's coordinates (DSL) to
; those specified by the COORD keyword.
; rot_m: Rotation matrix from the the specified coordinates to those
; defined by ROTATION.
; orient_m: Rotation matrix from the coordinates defined by ROTATION to
; the arbitrary coordinates defined by SLICE_ORIENT
; (column matrix of new coord's basis).
; (N/A for 2D interp)
; }
;
;
;CAVEATS: Due to IDL software constraints areas of no data are assigned zeros
; instead of NaNs. This yields a degeneracy between data gaps and
; measured zeros for both plots and numerical data.
; (Such degeneracy would always exist in any logarithmic plot made
; with the contour procedure)
;
;
;CREATED BY: A. Flores
; Based on work by Bryan Kerr and Arjun Raj
;
;EXAMPLES: see the crib file: thm_crib_part_slice2d.pro
;-
pro thm_part_slice2d, datArray, datArray2, datArray3, datArray4, $
; Time options
timewin=timewin, slice_time=slice_time, $
center_time=center_time, $
; Range limits
erange=erange, $
thetarange=thetarange, zdirrange=zdirrange, $
; Orientations
coord=coord, rotation=rotation, $
slice_orient=slice_orient, $
; Type options
type=type, two_d_interp=two_d_interp, $
three_d_interp=three_d_interp, $
; Other options
units=units, resolution=resolution, $
subtract=subtract,smooth=smooth,onecount=onecount, $
regrid=regrid, slice_width=slice_width, $
; Output
part_slice=part_slice, xgrid=xgrid, ygrid=ygrid, $
slice_info=slice_info, $
; Testing
vis=vis, $
; Other
_extra = _extra, fail=fail
compile_opt idl2
if size(datArray,/type) ne 8 then begin
fail = 'Invalid data structure. Input must be valid array of '+ $
'ESA or SST distributions with required tags. '+ $
'See thm_part_dist_array.'
dprint, fail
return
endif
if ~keyword_set(slice_time) then begin
fail = 'Please specifiy a time at which to compute the slice.'
dprint, fail
return
endif
if ~keyword_set(timewin) then begin
fail = 'Please specifiy a time window for the slice."
dprint, fail
return
endif
d_set = [keyword_set(datarray), keyword_set(datarray2), $
keyword_set(datarray3), keyword_set(datarray4)]
if total(d_set) gt 1 and keyword_set(type) then begin
fail = 'Multiple distributions are only supported for the '+ $
'geometric method.'
dprint, fail
return
endif
fail = ''
; Defaults
;---------
if ~keyword_set(coord) then coord='DSL'
if ~keyword_set(rotation) then rotation='xy'
if ~keyword_set(units) then units = 'df'
if ~keyword_set(slice_orient) then slice_orient = [0,0,1.,0]
if rotation eq 'xyz' then rotation = 'xy'
if rotation eq 'perp_xyz' then rotation = 'perp_xy'
probe = keyword_set(datArray[0].spacecraft) ? datArray[0].spacecraft : $
strmid(datArray[0].project_name, 0, /reverse_offset)
if ~stregex(probe, '[abcde]', /bool, /fold_case) then probe = ''
;Type specific:
if keyword_set(two_d_interp) then type = 2
if keyword_set(three_d_interp) then type = 3
if size(type,/type) eq 0 then type = 0
if type gt 3 then type = 0
; 2d interp
if type eq 2 then begin
if size(smooth,/type) eq 0 then smooth = 3
if size(resolution,/type) eq 0 then resolution = 51
if ~keyword_set(thetarange) and ~keyword_set(zdirrange) then begin
thetarange = [-20.,20]
endif
endif
; 2d nn
if type eq 1 then begin
if size(resolution,/type) eq 0 then resolution = 250
if keyword_set(regrid) && n_elements(regrid) eq 3 then begin
regrid = float(regrid)
endif else $
regrid = [64.,32,64]
endif
; geometric
if ~keyword_set(type) then begin
;nothing for now
endif
;resolution
if ~keyword_set(resolution) then begin
resolution = keyword_set(type) ? 150:500
endif
dprint, dlevel=2, 'Processing slice at ' + time_string(slice_time, format=5) +'...'
; Account for multiple data types
; -------------------------------
ds = [ptr_new(datArray), ptr_new(datArray2), $
ptr_new(datArray3), ptr_new(datArray4)]
ds = ds[where(d_set)]
; Get time info
;--------------
; get the boundaries of the time window
if keyword_set(center_time) then begin
tmin = slice_time - timewin/2
tmax = slice_time + timewin/2
endif else begin
tmin = slice_time
tmax = slice_time + timewin
endelse
; Cut dist. arrays to time limits
for i=0, n_elements(ds)-1 do begin
; get indexes of dat structures in reqested time window
times = (*ds[i]).time + ((*ds[i]).end_time - (*ds[i]).time)/2
times_ind = where(times ge tmin AND times le tmax, ndat)
if ndat gt 0 then begin
*ds[i] = (*ds[i])[times_ind]
endif else begin
dprint, dlevel=4, 'No particle distributions exist in data array '+ $
strtrim(i,2)+' that are within the requested time window.'
ds[i] = ptr_new() ;no need to free pointer
continue
endelse
endfor
; Remove arrays that were out of range
intrange = where(ptr_valid(ds), cin)
if cin gt 0 then begin
ds = ds[intrange]
endif else begin
fail = 'There are no distributions within the given time window ('+ $
time_string(slice_time)+ ' + '+strtrim(timewin,2)+' sec).'
dprint, fail
return
endelse
; Produce a concatenated list of cartesian-based
; velocity bins by looping over data types
;-------------------------------
for i=0, n_elements(ds)-1 do begin
;Get cartesian velocities
thm_part_slice2d_getXYZ, *ds[i], units, type, $
onecount=onecount, regrid=regrid, erange=erange, $
orange=org, vrange=vrange, $
velocities=v0, $ ;array of velocities
datapoints=dpoints, $ ;arraw of data
dp=dp0, dt=dt0, dv=dv0, $ ;dphi, dtheta, dr
fail=fail
if keyword_set(fail) then return
;Concatenate velocities from different distributions
if keyword_set(velocities) then begin
velocities = [ velocities, temporary(v0)]
datapoints = [datapoints, temporary(dpoints)]
endif else begin
velocities = temporary(v0)
datapoints = temporary(dpoints)
endelse
;concatenate extra variables for geometric method
if keyword_set(dp0) then begin
if keyword_set(dp) then begin
dp = [ dp, temporary(dp0)]
dt = [ dt, temporary(dt0)]
dv = [ dv, temporary(dv0)]
endif else begin
dp = dp0
dt = dt0
dv = dv0
endelse
endif
orange = keyword_set(orange) ? minmax([orange,org]):org
endfor
;pull last struct from the array to use later
last = (*ds[0])[n_elements(*ds[0])-1]
; For testing only: shows 3-D distribution of nonzero data
if keyword_set(vis) then begin
rr = where(datapoints ne 0)
rg = [-max(velocities),max(velocities)]
iplot, velocities[rr,0], velocities[rr,1], velocities[rr,2], linestyle=6, sym_index=1, $
xrange = rg, yrange=rg, zrange=rg, $
rgb_table=13, vert_colors=bytscl(alog10(datapoints[rr]))
endif
; Apply coordinate transforms and rotations
;------------------------------------------
; Get mag data
bfield = thm_part_slice2d_getmag(last, *ds[0], magdata=magdata, rotation=rotation, fail=fail)
if n_elements(bfield) eq 1 then return
; Get bulk velocity (in m/s)
vbulk = thm_part_slice2d_getvel(last, *ds[0], veldata=veldata, rotation=rotation, fail=fail)
if n_elements(vbulk) eq 1 then return
; Subtract bulk velocity vector?
if keyword_set(subtract) and type ne 0 then begin
thm_part_slice2d_subtract, velocities, vbulk, veldata=veldata
endif
; Rotate to GSM or GSE if requested
thm_part_slice2d_cotrans, velocities, coord, bfield, vbulk, type, probe=probe, $
twin=[tmin,tmax], ct=ct, fail=fail
if keyword_set(fail) then return
; Rotate data to desired coordinates
thm_part_slice2d_rotate, velocities, rotation, bfield, vbulk, type, rot=rot, fail=fail
if keyword_set(fail) then return
; Convert velocities to km/s
factor = 1000.
velocities = temporary(velocities)/factor
vbulk = vbulk/factor
if keyword_set(dv) then dv = dv/factor
if keyword_set(vrange) then vrange = vrange/factor
; rr = where(datapoints ne 0)
; rg = [-max(velocities),max(velocities)]
; iplot, velocities[rr,0], velocities[rr,1], velocities[rr,2], linestyle=6, sym_index=1, $
; xrange = rg, yrange=rg, zrange=rg, $
; rgb_table=13, vert_colors=bytscl(alog10(datapoints[rr]))
; Misc.
;------
; Apply theta or z-axis limits (for old method, see documentation)
thm_part_slice2d_2dlimits, velocities, datapoints, thetarange=thetarange, zdirrange=zdirrange, fail=fail
if keyword_set(fail) then return
; Sort transformed velocity grid
sorted = sort(velocities[*,0])
velocities = velocities[sorted,*]
if keyword_set(dp) then begin
dp = dp[sorted]
dt = dt[sorted]
dv = dv[sorted]
endif
datapoints = datapoints[sorted]
; Create slice:
;
; TYPE=0 Geometric Method
; TYPE=1 2D Nearest Neighbor Interpolation
; TYPE=2 2D Linear Interpolation
; TYPE=3 3D Linear Interpolation
;------------------------------------------
; Linear 2D interpolation from thm_esa_slice2D rewritten
if type eq 2 then begin
dprint, dlevel=4, 'Using 2d linear interpolation'
thm_part_slice2d_2di, datapoints, velocities, resolution, $
part_slice=part_slice, $
xgrid=xgrid, ygrid=ygrid, $
fail=fail
; 3D Interpolated slices
endif else if type eq 3 then begin
dprint, dlevel=4, 'Using 3d linear interpolation'
thm_part_slice2d_3di, datapoints, velocities, resolution, $
slice_orient, orange, $
part_slice=part_slice, $
xgrid=xgrid, ygrid=ygrid, $
fail=fail
; 2D Nearest Neighbor (for testing only)
endif else if type eq 1 then begin
dprint, dlevel=4, 'Using 2d nearest neighbor'
thm_part_slice2d_nn, datapoints, velocities, resolution, slice_orient, $
part_slice=part_slice, $
xgrid=xgrid, ygrid=ygrid, $
slice_width=slice_width, $
shift=keyword_set(subtract) ? vbulk:0, $
fail=fail
; Geometric Method
endif else begin
dprint, dlevel=4, 'Using geometric method'
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
endelse
if keyword_set(fail) then begin
dprint, fail
return
endif
; Apply smoothing
if keyword_set(smooth) then begin
smooth = round(smooth)
if smooth ge 2 then begin
part_slice = smooth(part_slice, smooth)
negidx = where(part_slice lt 0, nneg)
if nneg gt 0 then begin
; dprint, 'Negative values created from smoothing. Removing...'
part_slice[negidx] = 0
endif
endif
endif
; Pass out slice information for plotting
;----------------------------------------
;Misc
tr = [min((*ds[0]).time), last.end_time]
for i=0, n_elements(ds)-1 do begin
;dist type name
dname = keyword_set(dname) ? [dname,((*ds[i]).data_name)[0]]:((*ds[i]).data_name)[0]
;time range of data
tr = [ tr[0] < min((*ds[i]).time), $
tr[1] > max((*ds[i]).end_time) ]
endfor
dname = strjoin(dname,'/')
;Ensure user defined rotations are accounted for and passed out
if type ne 2 then begin
thm_part_slice2d_orientslice, slice_orient, matrix=matrix, fail=fail
if keyword_set(vbulk) then vbulk = matrix ## vbulk
endif
;Velocity range
if keyword_set(slice_orient) then begin
vrange = sqrt( vrange^2 - slice_orient[3]^2 > 0.)
endif
slice_info = {probe: probe, $
dist: dname, $
coord: coord, $
rot: rotation, $
units: units, $
twin: timewin, $
type: type, $
range: orange, $
trange: tr, $ ;time range of slice data
ndists: n_elements((*ds[0])), $ ;ndat ;# distributions used
vrange: vrange , $
shift: keyword_set(subtract) ? vbulk:0, $
bulk: keyword_set(vbulk) ? vbulk:0, $
coord_m: keyword_set(ct) ? ct:-1, $ ;coordinate transform from DSL
rot_m: rot, $ ;rotation matrix
orient_m: keyword_set(matrix) ? matrix:-1 $ ;slice plane orientation matrix
}
return
end