; 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