;+ ;Procedure: ; thm_part_slice2d ; ;Purpose: ; Wrapper for using spd_slice2d to create particle distribution contours. ; This essentially performs all THEMIS-specific operations that occured ; in the original thm_part_slice2d and should be backward-compatible. ; ; Returns a 2-D slice of THEMIS ESA/SST particle distributions. ; This procedure works in conjunction with thm_part_dist_array.pro and ; thm_part_slice2d_plot.pro. ; ; There are three methods for generating slices: ; ; Geometric: ; Each point on the plot is given the value of the bin it instersects. ; This allows bin boundaries to be drawn at high resolutions. ; ; 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. ; ; 3D Interpolation: ; The entire 3-dimensional distribution is linearly interpolated onto a ; regular 3d grid and a slice is extracted from the volume. (default) ; ; ;Calling Sequence: ; thm_part_slice2d, datArr, [datArr2, [datArr3, [datArr4]]], $ ; timewin = timewin, slice_time=slice_time, $ ; part_slice=part_slice ; ; ;Arguments: ; DATARR[#]: An array of pointers to 3D data structures. ; See thm_part_dist_array.pro for more. ; ; ;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. ; TIMEWIN: Length in seconds over which to compute the slice. ; CENTER_TIME: Flag that, when set, centers the time window around the time ; specified by SLICE_TIME keyword. ; UNITS: A string specifying the units to be used. ; ('counts', 'DF' (default), 'rate', 'crate', 'flux', 'eflux') ; ENERGY: Flag to plot data against energy (in eV) instead of velocity. ; LOG: Flag to apply logarithmic scaling to the radial mesure (i.e. energy/velocity). ; (on by default if /ENERGY is set) ; ; GEOMETRIC: Flag to use geometric method (described above) ; TWO_D_INTERP: Flag to use 2D interpolation method (described above) ; THREE_D_INTERP: Flag to use 3D interpolation method (described above, the default) ; ; COORD: A string designating the coordinate system in which the slice will be ; oriented. Options are 'DSL', 'GSM', 'GSE' and the following magnetic ; field aligned coordinates (field parallel to z axis). ; ; 'xgse': The x axis is the projection of the GSE x-axis ; 'ygsm': The y axis is the projection of the GSM y-axis ; 'zdsl': The y axis is the projection of the DSL z-axis ; 'RGeo': The x is the projection of radial spacecraft position vector (GEI) ; 'mRGeo': The x axis is the projection of the negative radial spacecraft position vector (GEI) ; 'phiGeo': The y axis is the projection of the azimuthal spacecraft position vector (GEI), positive eastward ; 'mphiGeo': The y axis is the projection of the azimuthal spacecraft position vector (GEI), positive westward ; 'phiSM': The y axis is the projection of the azimuthal spacecraft position vector in Solar Magnetic coords ; 'mphiSM': The y axis is the projection of the negative azimuthal spacecraft position vector in Solar Magnetic coords ; ; ; ROTATION: The rotation keyword specifies the orientation of the slice plane ; within the given coordinates (BV and BE will be invariant between ; coordinate systems). ; ; '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 coordinate's x axis and y is along the coordinate's y axis ; 'xz': The x axis is along the coordinate's x axis and y is along the coordinate's z axis ; 'yz': The x axis is along the coordinate's y axis and y is along the coordinate's z axis ; 'xvel': The x axis is along the coordinate'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 coordinate's x & y axes are projected onto the plane normal to the B field ; 'perp_xz': The coordinate's x & z axes are projected onto the plane normal to the B field ; 'perp_yz': The coordinate'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 ; COORD and ROTATION. Both keywords take 3-vectors as input. ; ; 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 plane perpendicular to DSL z-axis using [3,2,0] as plane's x-axis: ; (this is the same as only using SLICE_X=[3,2,1]) ; COORD='dsl' (default), ROTATION='xyz' (default), SLICE_X=[3,2,1] ; ; Slice plane perp. to GSE x-axis, bulk velocity used to define plane's x-axis: ; COORD='gse', ROTATION='xvel', SLICE_NORM=[1,0,0], SLICE_X=[0,1,0] ; ; Slice plane along the B field and radial position vectors, B field used as slice's x-axis: ; COORD='rgeo', SLICE_NORM=[0,1,0], SLICE_X=[0,0,1] ; 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] ; ; ; VEL_DATA: Name of tplot variable containing the bulk velocity data. ; This will be used for slice plane alignment and subtraction. ; If not set the bulk velocity will be automatically calculated ; from the distribution (when needed). ; MAG_DATA: Name of tplot variable containing magnetic field data. ; This will be used for slice plane alignment. ; ERANGE: Two element array specifying the energy range to be used. ; COUNT_THRESHOLD: Mask bins that fall below this number of counts after averaging. ; (e.g. COUNT_THRESHOLD=1 masks bins with counts below 1) ; 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 the smoothing window in # ; of points. Even entries will be incremented, 0 and 1 are ignored. ; Smoothing is performed with a gaussian convolution. ; ; 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 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 x-axis; ; positive in the right handed direction. This will ; average over all data within that range. ; e.g. [-25,25] will average data within 25 degrees ; of the slice plane about it's x-axis ; ; MSG_OBJ: Object reference to GUI message bar. If included useful ; console messages will also be output to GUI. ; ; ;Output: ; PART_SLICE: Structure to be passed to thm_part_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 ; n_samples: number of distributions averaged to create slice ; ; mass: partile 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=geo, 2=2D interp, 3=3D interp); ; energy: flag that x and y are energy instead of velocity ; 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 bfiend 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 containting 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: ; A. Flores Based on work by Bryan Kerr and Arjun Raj, and Xuzhi Zhou ; ; ;EXAMPLES: ; See the crib file: thm_crib_part_slice2d.pro ; ; ;$LastChangedBy: nikos $ ;$LastChangedDate: 2020-06-09 16:15:53 -0700 (Tue, 09 Jun 2020) $ ;$LastChangedRevision: 28772 $ ;$URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/spdsoft/tags/spedas_4_0/projects/themis/spacecraft/particles/slices/thm_part_slice2d.pro $ ;- pro thm_part_slice2d, ptrArray, ptrArray2, ptrArray3, ptrArray4, $ ; Time options timewin=timewin, slice_time=slice_time_in, $ center_time=center_time, $ ; Range limits erange=erange, $ thetarange=thetarange, zdirrange=zdirrange, $ ; Orientations coord=coord_in, rotation=rotation, $ slice_x=slice_x, slice_norm=slice_z, $ displacement=displacement_in, $ ; Support Data mag_data=mag_data, vel_data=vel_data, $ ; Type options type=type, $ geometric=geometric, $ two_d_interp=two_d_interp, $ three_d_interp=three_d_interp, $ ; Other options units=units, resolution=resolution, $ average_angle=average_angle, smooth=smooth, $ count_threshold=count_threshold, $ subtract_counts=subtract_counts, $ subtract_bulk=subtract_bulk, $ log=log, energy=energy, $ ; Output part_slice=slice, $ ; Dummy vars to conserve backwards compatibility ; (should allow oldscripts to be run withouth issue) xgrid=xgrid, ygrid=ygrid, slice_info=slice_info, $ onecount=onecount, $ ; Other fix_counts=fix_counts, $ msg_obj=msg_obj, fail=fail, $ _extra = _extra compile_opt idl2 fail = '' if size(ptrArray,/type) ne 10 then begin fail = 'Invalid data structure. Input must be valid array of pointers to arrays of '+ $ 'ESA or SST distributions with required tags. '+ $ 'See thm_part_dist_array.' dprint, dlevel=1, fail return endif if undefined(slice_time_in) then begin fail = 'Please specifiy a time at which to compute the slice.' dprint, dlevel=1, fail return endif if undefined(timewin) then begin fail = 'Please specifiy a time window for the slice." dprint, dlevel=1, fail return endif valid_coords = ['dsl','gse','gsm','xgse','ygsm','zdsl','rgeo', $ 'mrgeo','phigeo','mphigeo','phism','mphism'] if ~undefined(coord_in) then begin coord = strlowcase(coord_in) if ~in_set(coord,valid_coords) then begin fail = 'Invalid coordinates requested. See thm_crib_part_slice2d for examples.' dprint, dlevel=1, fail return endif endif ; Defaults ;------------------------------------------------------------ slice_time = time_double(slice_time_in[0]) if undefined(coord) then coord='dsl' if undefined(units) then units = 'df' ;Interpolation type: if undefined(type) then begin 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 endif else begin if type[0] eq 0 then geometric = 1 if type[0] eq 2 then two_d_interp = 1 if type[0] eq 3 then three_d_interp = 1 type = type[0] endelse if undefined(type) then type = 3 ;3D Interpolation is the dafault probe = keyword_set((*ptrArray[0])[0].spacecraft) ? (*ptrArray[0])[0].spacecraft : $ strmid((*ptrArray[0])[0].project_name, 0, /reverse_offset) if ~stregex(probe, '[abcde]', /bool, /fold_case) then probe = '' ;preserve backward compatability if ~keyword_set(count_threshold) && keyword_set(onecount) then count_threshold = onecount ; Get data & info ;------------------------------------------------------------ ;aggregate pointers here to make the unit conversion easier ds = ptrArray if ptr_valid(ptrArray2) then ds = [ds,ptrArray2] if ptr_valid(ptrArray3) then ds = [ds,ptrArray3] if ptr_valid(ptrArray4) then ds = [ds,ptrArray4] ; get the boundaries of the time window to retrieve support data if keyword_set(center_time) then begin trange = [slice_time - timewin/2, $ slice_time + timewin/2 ] endif else begin trange = [slice_time, $ slice_time + timewin ] endelse ; Get any support data that's stored in the distribution structures ; (backward compatibility) ;------------------------------------------------------------ ; Get bulk velocity calculated from dist struct if none is specified if undefined(vel_data) then begin vel_data = thm_part_slice2d_getvel(ds, trange=trange, energy=energy, fail=fail) endif ; Get sun vector if present in structures sun_data = thm_part_slice2d_getsun(ds, trange=trange, fail=fail) if n_elements(sun_data) ne 3 then undefine, sun_data ; Get THEMIS specific rotation matrices ;------------------------------------------------------------ if in_set(coord,['dsl','gse','gsm']) then begin thm_part_slice2d_cotrans, probe=probe, coord=coord, trange=trange, fail=fail, $ matrix=custom_rotation endif else begin thm_part_slice2d_fac, probe=probe, coord=coord, trange=trange, fail=fail, $ mag_data=mag_data, matrix=custom_rotation endelse if keyword_set(fail) then return ; Process the data and generate the slice ;------------------------------------------------------------ ; perform unit conversion, remove sst contamination, apply eclipse corrections, etc thm_part_process, ds, processed, trange=trange, units=units, _extra=_extra slice = spd_slice2d(processed, $ ; Time options time=slice_time, $ window=timewin, $ center_time=center_time, $ ; Orientations custom_rotation=custom_rotation, $ rotation=rotation, $ slice_norm=slice_z, $ slice_x=slice_x, $ 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, $ ; Other msg_obj=msg_obj, $ fail=fail, $ _extra = _extra) if keyword_set(fail) then begin fail = 'Error producing slice: ' + fail dprint, dlevel=0, fail return endif ;set coordinates for plot labels slice.coord=coord ; Apply N-count threshold/subtraction if requested ; -create copy of the data set to N counts ; -produce identical slice from the synthetic data ; -use new slice as mask ;------------------------------------------------------------ if keyword_set(count_threshold) or keyword_set(subtract_counts) then begin fix_counts = keyword_set(subtract_counts) ? subtract_counts:count_threshold thm_part_copy, ds, fixed thm_part_set_counts, fixed, fix_counts, /set_units ;no background or contamination removal should be performed here! thm_part_process, fixed, fixed_processed, units=units, esa_bgnd_remove=0 mask = spd_slice2d(fixed_processed, $ ; Time options time=slice_time, $ window=timewin, $ center_time=center_time, $ ; Orientations custom_rotation=custom_rotation, $ rotation=rotation, $ slice_norm=slice_z, $ slice_x=slice_x, $ 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, $ ; Other msg_obj=msg_obj, $ fail=fail, $ _extra = _extra) if keyword_set(fail) then begin fail = 'Error calculating count threshold: ' + fail dprint, dlevel=0, fail return endif if keyword_set(count_threshold) then begin btidx = where(slice.data lt mask.data,nbt) if nbt gt 0 then begin slice.data[btidx] = 0 endif endif else if keyword_set(subtract_counts) then begin slice.data = (slice.data - mask.data) > 0 endif endif end