;+ ; ; Name: ; spd_slice2d_geo.pro ; ; Purpose: ; Helper function for spd_slice2d.pro ; ; Produces slices showing each bin's boundaries by assigning ; each bin's value to all points on the slice plane that ; fall within that bin's boundaries. ; ; This is essentially a hack to allow plotting of bin boundaries ; with contour instead of creating a new plotting routine. ; ; Input: ; data: N element array of data values ; rad: N element array of radial values ; phi: N element array of phi values ; theta: N element array of theta values ; ; dp: N element array of phi ranges ; dt: N element array of theta ranges ; dr: N element array of velocity ranges ; ; resolution: Single value (R) giving the number of points in each ; dimension of the slice ; average_angle: 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. ; e.g. [-25,25] will average data within 25 degrees ; of the slice plane about it's x-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 x-axis. ; e.g. [-25,25] will sum data within 25 degrees ; of the slice plane about it's x-axis ; ; custom_matrix: Rotation matrix from native -> user specified coordinates ; (applied first) ; rotation_matrix: Rotation matrix from given coordinates to built in ; rotated coordinates, e.g. GSM -> BV, perp_xy ; (applied second) ; slice_matrix: Rotation matrix from specified coords/rotation into ; the slice plane's final coordinates as defined by ; the user specified normal and x projection. ; (applied last) ; shift: Vector by which the slice should be shifted (e.g. bulk velocity subtraction). ; The slice plane will be shifted by the z value and the x & y values will be ; subtracted from the corresponding axes' grids (it should already be in the ; slice plane's coordinates). ; ; Output: ; xgrid: R element array of x-axis values for the slice ; ygrid: R element array of y-axis values for the slice ; slice: RxR element array containing the slice data ; ; Other Keywords: ; msg_obj: dprint display object reference ; msg_prefix: String prefix to be printed with progress messages. ; ; Notes: ; -This routine will slow as the number of bins (N) increases. ; Averaging will significantly lengthen the required time. ; ; $LastChangedBy: egrimes $ ; $LastChangedDate: 2018-04-12 13:51:19 -0700 (Thu, 12 Apr 2018) $ ; $LastChangedRevision: 25039 $ ; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/spdsoft/tags/spedas_5_0/general/science/spd_slice2d/core/spd_slice2d_geo.pro $ ;- pro spd_slice2d_geo, data=data, resolution=resolution, $ rad=r, phi=phi, theta=theta, $ dr=dr, dp=dp, dt=dt, $ custom_matrix=ct, rotation_matrix=rot, orient_matrix=mt, $ shift=shift, $ average_angle=average_angle, $ sum_angle=sum_angle, $ ; Data Output slice=slice, xgrid=xgrid, ygrid=ygrid, $ ; Info Output fail=fail, msg_obj=msg_obj, msg_prefix=msg_prefix, $ ; Other _extra=_extra compile_opt idl2, hidden ;for progress messages previous_time = systime(/sec) n = float(resolution) rd = 180./!dpi tolerance = 5d-7 ;Initialize slice and coordinates slice = dblarr(n,n) ;Create grid of coordinates for entire slice plane vrange = max(abs(r))*[-1,1] xgrid = interpol(vrange + [-1,1]*max(dr), n) ygrid = interpol(vrange + [-1,1]*max(dr), n) z = undefined(shift) ? 0. : shift[2] u = [ [reform(xgrid # replicate(1.,n), n^2)], $ ;x [reform(replicate(1.,n) # ygrid, n^2)], $ ;y [reform(replicate(z,[n,n]), n^2)] ] ;z ;Rotate slice coordinates to desired location. ;The "ct" and "rot" matrices transform INTO the slice plane's coordinates. ;To determine which bins intersect the slice plane the transformation ;is reversed and the applied to the slice itself to transform it ;into the data's native coordinates. m = mt if keyword_set(ct) then begin if keyword_set(rot) then m = (ct # rot) # m $ else m = ct # m endif else begin if keyword_set(rot) then m = rot # m endelse u = u # transpose(m) ;Get necessary rotations for averaging/summing. ;The coordinates of the slice plane will be rotated through the ;specified angle range. A separate search will be done at each ;new plane. The number of planes is determined from the desired ;resolution and the width of the average. if keyword_set(average_angle) or keyword_set(sum_angle) then begin if keyword_set(sum_angle) then alpha = minmax(sum_angle) if keyword_set(average_angle) then alpha = minmax(average_angle) ;number of additional slices to average/sum over na = (2 * sqrt(n) * (alpha[1]-alpha[0])/90.) > 2 ;copy slice's x-vector xv = m[*,0] ## replicate(1d,na) ;interpolate across the angle range a = ([dindgen(na-1)/(na-1),1] * (alpha[1]-alpha[0])) + alpha[0] ;construct quaternion array to get rotation matrices qs = qcompose(xv,a/rd, /free) ;quaternions to rotate about x by a ms = qtom(qs) ;get matrices if n_elements(ms) eq 1 then begin fail = 'Error: Cannot construct rotation matrices for angular averaging.' dprint, dlevel=0, fail return endif ; Testing vectorization... ; ut = u ; for i=0, na-1 do begin ; ut = [ [temporary(ut)], [qrotv2(qs[i,*],u)] ] ; endfor ; ; u = temporary(ut) ; ; ix = lindgen(na)*3 ; iy = ix+1 ; iz = ix+2 ; ; ;Convert transformed slice coordinates to spherical ; pcoords = rd * atan(u[*,iy],u[*,ix]) ;phi ; tcoords = rd * atan(u[*,iz], sqrt(u[*,ix]^2 + u[*,iy]^2)) ;theta ; rcoords = sqrt(u[*,ix]^2 + u[*,iy]^2 + temporary(u[*,iz]^2));r endif else na = 0 ;simplify first loop below ;Loop over slice planes (if averaging over angle) ;A vectorized method was tested and took longer. weight = bytarr(size(slice,/dim)) np = n_elements(data) for j=-1, na-1 do begin ut = j ge 0 ? reform(ms[j,*,*]) ## u:u ;Convert transformed slice coordinates to spherical pcoords = rd * atan(ut[*,1],ut[*,0]) ;phi tcoords = rd * atan(ut[*,2], sqrt(total(ut[*,0:1]^2,2)) ) ;theta rcoords = sqrt(ut[*,0]^2 + ut[*,1]^2 + ut[*,2]^2) ;r ;Loop over bins to determine what region each bin covers on the slice plane. for i=0, np-1 do begin if data[i] eq 0 then continue ; Theta-- tlim = [ theta[i]-0.5*dt[i], theta[i]+0.5*dt[i] ] ;account for rounding errors ;this is particularly important is slice plane is at zero elevation tr = where( abs(tlim - round(tlim)) lt tolerance, ntr) if ntr gt 0 then tlim[tr] = round(tlim[tr]) t = (tcoords gt tlim[0]) and (tcoords le tlim[1]) if total(t) lt 1 then continue ; Phi-- plim = [ phi[i]-0.5*dp[i], phi[i]+0.5*dp[i] ] ; keep limits within [-180,180] over = where(plim gt 180, no) under = where(plim lt -180, nu) if no gt 0 then plim[over] += -360 if nu gt 0 then plim[under] += 360 ;account for rounding errors pr = where( abs(plim - round(plim)) lt tolerance, npr) if npr gt 0 then plim[pr] = round(plim[pr]) ;determine which region ( p0->p1 or p1->p0) the bin spans if plim[0] gt plim[1] then begin p = (pcoords gt plim[0]) or (pcoords le plim[1]) endif else begin p = (pcoords gt plim[0]) and (pcoords le plim[1]) endelse if total(p) lt 1 then continue ; R (velocity/energy)-- s = (rcoords ge r[i]-0.5*dr[i]) and (rcoords lt r[i]+0.5*dr[i]) ; Combine criteria and assign values-- bidx = where( p and t and s, nb) if nb gt 0 then begin weight[bidx]++ slice[bidx] += data[i] endif ; Output progress messages every 6 seconds-- if systime(/sec) - previous_time gt 6 then begin msg = strtrim(long( 100.*((j+1)*np + i)/((na>1)*np) ),2) + '% complete' if keyword_set(msg_prefix) then msg = msg_prefix + msg dprint, dlevel=2, msg, display_object=msg_obj previous_time = systime(/sec) endif endfor endfor ;Average areas where bins overlapped adj = where(weight eq 0, na) if na gt 0 then weight[adj] = 1b if undefined(sum_angle) then slice = slice / weight ;Align the x & y axis values if there was a shift if ~undefined(shift) then begin xgrid -= shift[0] ygrid -= shift[1] endif return end