;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; plot_3d_vect.pro ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Author: wd (Wolfgang.Dobler@kis.uni-freiburg.de)
;;; Date: 18-Aug-2000
;;; Version: 1.37
;;; CVS: $Revision: 1.8 $
;;;
;;; Description:
;;; Given either a 3d array on a 2d grid (e.g. B(0:64,0:64,15,0:3)),
;;; or three 1d arrays on a 2d grid (e.g. Bx(0:64,0:64,15),
;;; By(..), Bz(..)), plot the vertical component (the third cone) as
;;; coloured background and overlay the horizontal components with
;;; wdvelovect.
;;;
;;; Usage:
;;; PLOT_3D_VECT, VEC, [X, Y,]
;;; SAMPLE=sample,
;;; PERMUTE=permute, AUTOXYTITLE=autoxy, $
;;; PS_RGB=ps_rgb, PS_COL_MIN=ps_col_min, $
;;; TV_RGB=tv_rgb, $
;;; PS_ENH=ps_enh, $
;;; BITS=bits, $
;;; ZRANGE=zrange, $
;;; /KEEP_CT, $
;;; /TRUE_COLOR, $
;;; /PS_FULL_RANGE
;;;
;;; or
;;; PLOT_3D_VECT, VEC_x, VEC_y, VEC_z, [X, Y,] $
;;; .........
;;;
;;; Arguments:
;;; VEC --- An array of size (NX, NY, 3) containing the vector
;;; field to be plotted
;;; VEC_X/Y/Z --- Three arrays of size (NX,NY) containing the vector
;;; field to be plotted
;;; X, Y --- Arrays of size (NX) and (NY) containing the
;;; (equidistant) coordinates of the grid axis
;;;
;;; Keywords:
;;; SAMPLE --- A scalar or 1-d array containing a factor by
;;; which the vec_x and vec_y should be downsampled.
;;; PERMUTE --- A 3-d array, specifying the order in which the
;;; data array components are to be used, e.g.
;;; plot_3d_vect, uu[*,0,*], PERMUTE=[0,2,1], x, z
;;; will plot an x-z section
;;; AUTOXYTITLE --- Automatically generate XTITLE and YTITLE based
;;; on PERMUTE (and assuming the canonical order x-y-z)
;;; NEGATE3 --- Multiply 3rd component (the colour-coded one) with -1
;;; PS_RGB --- A (n_cols,3) array defining the red, green, and
;;; blue values of the PostScript colortable
;;; PS_COL_MIN --- Index of the first colour (usually the darkest
;;; one) to be used in PostScript.This is to make the
;;; overlayed arrows visible.
;;; TV_RGB --- Either: a (n_cols,3) array defining the red, green,
;;; and blue values of the computer display colour
;;; table;
;;; or: a number, indexing one of several predefined
;;; colour tables
;;; PS_ENH --- Factor by which the number of grid points should
;;; be enhanced (by interpolating data) for the
;;; vertical component
;;; BITS --- Colour depth for PostSCript device (default = 8)
;;; ZRANGE --- interval for z component (mapped to colours
;;; [0, ncolours])
;;; KEEP_CT --- Flag for keeping the colour table instead of
;;; creating a new one.
;;; TRUE_COLOR --- Flag for explicitly choosing TrueColor (if set to
;;; non-zero value) or PseudoColor (if set to zero).
;;; If this keyword is not set, plot_3d_vect
;;; determines the visual on its own.
;;; PS_FULL_RANGE --- Flag telling plot_3d_vect to use the full
;;; available colour range for PostScript. This
;;; will in general shift the colour for value zero
;;; XRANGE,YRANGE --- Set x- and y-range
;;; DEBUG --- Write diagnostic output
;;;
;;; Any other keywords are handed on to WDVELOVECT.
;;;
PRO plot_3d_vect, arg1, arg2, arg3, $
arg4, arg5, $
SAMPLE=sample, $
PERMUTE=perm, AUTOXYTITLE=autoxy, NEGATE3=negate3, $
PS_RGB=ps_rgb, PS_COL_MIN=ps_col_min, $
TV_RGB=tv_rgb, $
PS_ENH=ps_enh, $
BITS=bits, $
ZRANGE=zrange, $
KEEP_CT=keep_ct, $
TRUE_COLOR=true_col, $
PS_FULL_RANGE=ps_full_range, $
XRANGE=xrange, YRANGE=yrange, $
DEBUG=debug, HELP=help, _EXTRA=extra
COMMON display1, visual, called, tv_rgb_list ; Remember these variables
if (keyword_set(help)) then extract_help, 'plot_3d_vect'
; ON_ERROR, 2 ;Return to caller if an error occurs
pmulti0 = !p.multi
default, ps_col_min, 0
default, bits, 8
default, perm, [0,1,2]
default, negate3, 0
default, debug, 0
xytit = '!8' + ['x','y','z'] + '!X'
if (debug) then $
print, FORMAT='(A, 10(I3))', $
'PLOT_3D_VECT: Initial !p.multi = ', !p.multi
if (n_elements(xrange) gt 1) then begin
xbounds = xrange
endif else begin
xbounds = !x.range
endelse
if (n_elements(yrange) gt 1) then begin
ybounds = yrange
endif else begin
ybounds = !y.range
endelse
;; Make sure these ranges make sense
if ((xbounds[1]-xbounds[0]) eq 0) then xbounds = [-1e30,1e30]
if ((ybounds[1]-ybounds[0]) eq 0) then ybounds = [-1e30,1e30]
IF (N_ELEMENTS(called) LT 1) THEN BEGIN ; Do this just at the first call
if (!p.multi[0] eq 0) then begin
erase ; Sometimes necessary to get n_cols right
endif
ENDIF
n_cols = !D.TABLE_SIZE ; Number of colours
IF (NOT KEYWORD_SET(ps_rgb)) THEN BEGIN
ps_rgb = intarr(n_cols,3)
ps_rgb(*,0) = ps_col_min + $
findgen(n_cols)*(n_cols-ps_col_min-1.)/(n_cols-1.)
ps_rgb(*,1) = ps_rgb(*,0)
ps_rgb(*,2) = ps_rgb(*,1)
ENDIF
;; Set up a few colour tables
IF (N_ELEMENTS(called) LT 1) THEN BEGIN ; Do this just at the first call
tv_rgb_list = intarr(n_cols,3,2)
phi = (findgen(n_cols-2)+1)/(n_cols-1)*!PI/2
;; First colour table
tv_rgb_list(*,0,0) = byte([0, 255*sin(phi), 255])
tv_rgb_list(*,1,0) = byte([0, 255/sqrt(2.)*sin(2*phi)^2, 255])
tv_rgb_list(*,2,0) = byte([0, 255*cos(phi), 255])
;; Second colour table
tv_rgb_list(*,0,1) = [0, bindgen(n_cols-2)+1, 255]
tv_rgb_list(*,1,1) = [0, bindgen(n_cols/2-1)+1, $
n_cols-bindgen((n_cols-1)/2)-1-n_cols/2, 255]
tv_rgb_list(*,2,1) = [0, n_cols-bindgen(n_cols-2)-1, 255]
ENDIF
IF (N_ELEMENTS(tv_rgb) EQ 0) THEN BEGIN
tv_rgb = reform(tv_rgb_list(*,*,0))
ENDIF ELSE BEGIN
IF ((SIZE(tv_rgb))(0) EQ 0) THEN BEGIN ; If tv_rgb is scalar then..
IF (tv_rgb GE (size(tv_rgb_list))(3)) THEN BEGIN
MESSAGE, "Not enough colour tables in list."
ENDIF
tv_rgb = reform(tv_rgb_list(*,*,tv_rgb)) ; ..use it as index for the list
ENDIF
ENDELSE
IF (KEYWORD_SET(ps_full_range)) THEN absolut = 0 ELSE absolut = 1
;; Get the name of the visual from IDL
if (!d.name eq 'X') then begin
if ((n_elements(called) lt 1) or (n_elements(visual) le 0)) then begin
;; Do this just at the first call or if the first call was for a
;; non-X device
device, GET_VISUAL_NAME=visual
endif
if (visual eq 'TrueColor') then true_color = 1
endif
IF (N_ELEMENTS(true_col) NE 0) THEN BEGIN ; Overrides any previous choices
IF (true_col NE 0) THEN true_color = 1 ELSE true_color = 0
ENDIF
IF (N_ELEMENTS(true_color) EQ 0) THEN true_color = 0
IF (NOT KEYWORD_SET(keep_ct)) THEN keep_ct = 0
IF (keep_ct) THEN true_color = 0 ; Helps a lot
;; Process the data
vec_x = reform(arg1)
IF ((size(vec_x))[0] EQ 3) THEN BEGIN
;; If the first argument has dimensions (NX, NY, 3)
IF (N_ELEMENTS(arg3) NE 0) THEN BEGIN
;; three arguments: data_matrix, x, y
x = arg2
y = arg3
ENDIF
vec = vec_x
vec_x = vec(*,*,0)
vec_y = vec(*,*,1)
vec_z = vec(*,*,2)
ENDIF ELSE BEGIN ; First ARG1 has dimensions (NX, NY)
vec_x = arg1
vec_y = arg2
vec_z = arg3
IF (N_ELEMENTS(arg5) NE 0) THEN BEGIN
;; five arguments: data_x, data_y, data_z, x, y
x = arg4
y = arg5
ENDIF
ENDELSE
vec_all = [[[reform(vec_x)]], [[reform(vec_y)]], [[reform(vec_z)]]]
; vec_x = reform(vec_x)
; vec_y = reform(vec_y)
; vec_z = reform(vec_z)
;; Apply permutation
;; 1. Determine parity of permutation.
;; Even permutation (parity=+1) --> right-handed system
;; Odd permutation (parity=-1) --> left-handed system
parity = sign(perm[1]-perm[0])+sign(perm[2]-perm[1])+sign(perm[0]-perm[2])
if (negate3) then parity = -parity
vec_x = vec_all[*,*,perm[0]]
vec_y = vec_all[*,*,perm[1]]
vec_z = vec_all[*,*,perm[2]]*parity
v_size = size(vec_x)
nx = v_size(1)
ny = v_size(2)
nx_s = nx
ny_s = ny
IF (N_ELEMENTS(x) EQ 0) THEN x = findgen(nx)
IF (N_ELEMENTS(y) EQ 0) THEN y = findgen(ny)
x_s = congrid(x, nx_s)
y_s = congrid(y, ny_s)
;; Honour xrange, yrange:
goodx = where((x_s ge xbounds[0]) and (x_s le xbounds[1]))
if (goodx[0] lt 0) then message, 'No valid points within xrange'
x_s = x_s[goodx]
nx_s = n_elements(x_s)
vec_x = vec_x[goodx,*]
vec_y = vec_y[goodx,*]
vec_z = vec_z[goodx,*]
;
goody = where((y_s ge ybounds[0]) and (y_s le ybounds[1]))
if (goody[0] lt 0) then message, 'No valid points within yrange'
y_s = y_s[goody]
ny_s = n_elements(y_s)
vec_x = vec_x[*,goody]
vec_y = vec_y[*,goody]
vec_z = vec_z[*,goody]
dx_s = x_s(1)-x_s(0)
dy_s = y_s(1)-y_s(0)
IF (NOT KEYWORD_SET(ps_enh)) THEN ps_enh = 1
IF (!D.NAME EQ 'PS') THEN DEVICE, /COLOR, BITS=bits
called = 1
;;;;;;;
;; Establish coordinates and size of plot
xr = minmax(x_s)+dx_s*[-0.5,0.5]
yr = minmax(y_s)+dy_s*[-0.5,0.5]
idxr = (xr-x[0])/(x[1]-x[0]) ; index range to be plotted
idyr = (yr-y[0])/(y[1]-y[0]) ; index range to be plotted
;;
default, xrange, !x.range
default, yrange, !y.range
;; Do the sampling if needed
if (n_elements(sample) ne 0) then begin
samp_x = sample[0] > 1
if (n_elements(sample) eq 1) then begin
samp_y = sample[0] > 1
endif else begin
samp_y = sample[1] > 1
endelse
endif else begin
samp_x = 1
samp_y = 1
endelse
ixvel = indgen(nx_s/samp_x)*samp_x
iyvel = indgen(ny_s/samp_y)*samp_y
;; Awkward, but IDL does not support start:stop:stride notation, nor
;; does it allow for the array subscripts a la vec[ixvel,iyvel]
vec_x = vec_x[ixvel,*]
vec_x = vec_x[*,iyvel]
vec_y = vec_y[ixvel,*]
vec_y = vec_y[*,iyvel]
wdvelovect, vec_x, vec_y, x_s[ixvel], y_s[iyvel], $
XRANGE=xrange, YRANGE=yrange, /NOPLOT, _EXTRA=extra
sx = (xr[1]-xr[0])*!D.X_VSIZE*!X.S[1] ; Size ...
sy = (yr[1]-yr[0])*!D.Y_VSIZE*!Y.S[1] ; Size ...
xpos = (!X.S[1]*xr[0] + !X.S[0]) * !D.X_VSIZE + 1 ; ..and position of image
ypos = (!Y.S[1]*yr[0] + !Y.S[0]) * !D.Y_VSIZE + 1 ; ..and position of image
IF (!D.NAME EQ 'PS') THEN BEGIN
; TVLCT, ps_rgb(*,0), ps_rgb(*,1), ps_rgb(*,2)
if (n_elements(zrange) gt 0) then message, $
"Warning: keyword ZRANGE not yet implemented for PS output", $
/INFORMATIONAL
wdtvscl, $
interpolate(vec_z, $
(findgen(nx_s*ps_enh)+0.5)/ps_enh-0.5, $
(findgen(ny_s*ps_enh)+0.5)/ps_enh-0.5, $
/grid), $
xpos, ypos, $
XSIZE=sx, YSIZE=sy, $
BOTTOM=PS_COL_MIN, $
ABSOLUT=absolut
ENDIF ELSE BEGIN
IF ((NOT true_color) AND (NOT keep_ct)) THEN $
TVLCT, tv_rgb(*,0), tv_rgb(*,1), tv_rgb(*,2)
if (n_elements(zrange) eq 0) then begin
;; need a copy, so PLOT_#D_VECT, ..,ZRANGE=undef
;; will not set undef to a value if it was previously undefined
zrange1=[-1,1]*MAX([-MIN(vec_z), MAX(vec_z)])
endif else begin
zrange1 = zrange
endelse
x_enh = sx/(idxr[1]-idxr[0])
y_enh = sy/(idyr[1]-idyr[0])
nx1 = floor(nx_s*x_enh)
ny1 = floor(ny_s*y_enh)
image = 1+(interpolate(vec_z, $
(findgen(nx1)+0.5)/x_enh-0.5, $
(findgen(ny1)+0.5)/y_enh-0.5, $
/grid) $
-zrange1[0])/(zrange1[1]-zrange1[0])*(n_cols-3.)
IF (true_color) THEN BEGIN
true_image = lonarr(nx1,ny1,3)
FOR j = 0,2 DO BEGIN
true_image[*,*,j] = reform(tv_rgb[image,j], nx1, ny1)
ENDFOR
TV, true_image, xpos, ypos, true=3
ENDIF ELSE BEGIN
TV, image, xpos, ypos
ENDELSE
ENDELSE
pmulti1 = !p.multi
!p.multi = pmulti0
if (keyword_set(autoxy)) then begin
wdvelovect, vec_x , vec_y, x_s[ixvel], y_s[iyvel], /NOERASE, $
XRANGE=xrange, YRANGE=yrange, $
XTITLE=xytit[perm[0]], YTITLE=xytit[perm[1]], $
_EXTRA=extra
endif else begin
wdvelovect, vec_x , vec_y, x_s[ixvel], y_s[iyvel], /NOERASE, $
XRANGE=xrange, YRANGE=yrange, _EXTRA=extra
endelse
!p.multi = pmulti1
if (debug) then begin
print, FORMAT='(A, 10(I3))', $
'PLOT_3D_VECT: Final !p.multi = ', !p.multi
print, '------------------------------------------------'
endif
END
; End of file plot_3d_vect.pro