[FLASH-BUGS] FLASH / IDL 3D improvements
Robi Banerjee
banerjee at physics.mcmaster.ca
Thu Feb 13 10:34:18 CST 2003
Hey,
we did some improvements concerning the visualization of
3D flash data with your IDL routines:
1. plotting of a 1D slice
2. showing velocity vectors
3. manual setting of the resolution for a 2D slice
4. speed up of the merge_amr routine
We enclosed the new/modified IDL files; maybe you are interested to
integrate these routines to the next flash release?
Robi
====================================
Robi Banerjee
Department of Physics and Astronomy
ABB-320, McMaster University
Hamilton, ON L8S 4M1
CANADA
e-mail: banerjee at physics.mcmaster.ca
phone : (905) 525-9140 x 23189
fax : (905) 546-1252
-------------- next part --------------
;=============================================================================
; xflash -- general FLASH data plotting interface
;
; plot 1, 2, or 3-d data with velocity vectors, contours, and
; particle data overlayed.
;
; type xflash at the idl prompt to run
;
; modify xflash_defaults.pro to add additional problems
;
; xflash is an interface to the collection of idl routines which
; work with flash data. xflash collects the plot options and passes
; them on to xplot_amr, which does the actual plotting.
;
;
; modified:
; improved to handle 3D data 2003/02 robi banerjee
;=============================================================================
; event handler
pro xflash_event, ev
; common all the controls into the widget
common variables, varnames, native_num
common save, filename, pblm_ctr_lim, orientation, ndim
common save_color_info, current_color_id, current_color_name
common save_default_info, current_default_id, current_default_name
common save_contour_info, contours
common save_velocity_info, velocity
common save_histogram_info, nbins, max_scale
common save_number_info, current_number_id, current_number_name
common save_ptr, tree_ptr, data_ptr, params_ptr, data2d, $
xmerge, ymerge, zmerge
common save_path, current_path
common save_filename, filename_old, filename_last_read, old_start_sfx, $
old_end_sfx
widget_control, ev.top, get_uvalue=info
widget_control, ev.id, get_uvalue=uval
; take action pased on which widget was changed
case uval of
;-----------------------------------------------------------------------------
; file actions
;-----------------------------------------------------------------------------
'fopen': begin
filename = dialog_pickfile(FILTER='*hdf*', /MUST_EXIST, $
PATH=current_path)
; if no file was selected then skip every thing else
if (filename NE "") then begin
; update the current path for next time
pathEnd = strpos(filename, '/', /REVERSE_SEARCH)
current_path = strcompress(strmid(filename, 0, pathEnd))
baseNmEnd = strpos(filename, '_', /REVERSE_SEARCH)
filenameBase = strcompress(strmid(filename, 0, baseNmEnd+1))
oldBaseNmEnd = strpos(filename_old, '_', /REVERSE_SEARCH)
oldFilenameBase = strcompress(strmid(filename_old, 0, oldBaseNmEnd+1))
same_base = 0
if (oldFilenameBase EQ filenameBase) then same_base = 1
; make sure the file appears to be valid
itype = determine_file_type(filename)
if (itype EQ -1) then begin
result = dialog_message(['ERROR: file does not appear to be', $
'a valid FLASH data file.'], $
title='ERROR', $
dialog_parent=info.mainBase)
endif else begin
widget_control, info.finfo, sensitive=1
pathEnd = strpos(filename, '/', /reverse_search)
; ---- update the prototype file label ----------------------------------------
widget_control, info.fprototype, $
set_value='Prototype File: '+ $
strcompress(strmid(filename,pathEnd+1, $
strlen(filename)-pathEnd))
; ---- get the dimensionality -------------------------------------------------
ndim = get_dimensionality(filename)
; ---- setup some pointers to pass the data -----------------------------------
; if (n_elements(tree_ptr) GT 0) then ptr_free, tree_ptr
; if (n_elements(params_ptr) GT 0) then ptr_free, params_ptr
; if (n_elements(data_ptr) GT 0) then ptr_free, data_ptr
tree = {lrefine:0l, $
nodeType:0l, $
gid:lonarr(2*ndim+1+2^ndim), $
coord:fltarr(ndim), $
size:fltarr(ndim), $
bndBox:fltarr(2,ndim)}
tree_ptr = ptr_new(tree)
params = {totBlocks:0, $
corners:0, $
ndim:0, $
nvar:0, $
nxb:0, $
nyb:0, $
nzb:0, $
ntopx:1, $
ntopy:1, $
ntopz:1, $
time:0.0, $
dt:0.0}
params_ptr = ptr_new(params)
data_ptr = ptr_new(0.0)
; ---- get the suffix ---------------------------------------------------------
baseNmEnd = strpos(filename, '_', /REVERSE_SEARCH)
suffix = strcompress(strmid(filename,baseNmEnd+1,4))
; ---- sensitize the different plot options -----------------------------------
widget_control, info.fstartSuffix, sensitive=1
widget_control, info.fstartSuffix, SET_VALUE=suffix
widget_control, info.fendSuffix, sensitive=1
widget_control, info.fstep, sensitive=1
widget_control, info.out_buttons, sensitive=1
widget_control, info.hsize, sensitive=1
widget_control, info.vsize, sensitive=1
widget_control, info.opt_buttons, sensitive=1
widget_control, info.zoomLabel, sensitive=1
widget_control, info.varMin, sensitive=1
widget_control, info.varMax, sensitive=1
widget_control, info.autoData, sensitive=1
; ---- populate the variable drop box -----------------------------------------
; get the names of the unknowns stored in the file
varnames = get_var_list(filename)
if n_elements(varnames) EQ 1 then begin
if (varnames[0] EQ -1) then begin
; files from FLASH1 will not have the unknown names stored
result = dialog_message(['Warning: file does not have' , $
'variable name information.', $
'Unable to proceed . . .'], $
title = 'Warning', $
dialog_parent = info.mainBase)
endif
endif
if (same_base NE 1) then begin
; ---- check for particles ----------------------------------------------------
numParticles = get_particle_number(filename)
; ---- determine which derived variables are possible -------------------------
native_num = (size(varnames))[1]
; total velocity
totVelDefined = 0
case ndim of
1: begin
if (var_index('velx') NE -1) then begin
varnames = add_var(temporary(varnames), 'tot_vel')
totVelDefined = 1
endif
end
2: begin
if (var_index('velx') NE -1 AND $
var_index('vely') NE -1) then begin
varnames = add_var(temporary(varnames), 'tot_vel')
totVelDefined = 1
endif
end
3: begin
if (var_index('velx') NE -1 AND $
var_index('vely') NE -1 AND $
var_index('velz') NE -1) then begin
varnames = add_var(temporary(varnames), 'tot_vel')
totVelDefined = 1
endif
end
endcase
; sound speed
if var_index('pres') NE -1 AND $
var_index('dens') NE -1 AND $
var_index('gamc') NE -1 then $
varnames = add_var(temporary(varnames), 'snd_spd')
; mach number
if (var_index('pres') NE -1 AND $
var_index('dens') NE -1 AND $
var_index('gamc') NE -1 AND $
totVelDefined EQ 1) then $
varnames = add_var(temporary(varnames), 'mach')
; internal energy and ratio of kinetic to internal
if (var_index('ener') NE -1 AND $
totVelDefined EQ 1) then begin
varnames = add_var(temporary(varnames), 'int_ener')
varnames = add_var(temporary(varnames), 'ekin/eint')
endif
widget_control, info.vars, /destroy
; make sure to store the new instance of the droplist back to the info
; structure so it is available to others
info.vars = widget_droplist(info.varBase, $
title = 'Variables: ', $
uvalue = 'var', value = varnames)
; droplist for resolution
widget_control, info.resolution, /destroy
info.resolution = widget_droplist(info.varBase, $
title = 'Resolution: ', $
uvalue = 'res', $
value = ['0','1','2','3'])
; ---- load the generic defaults ----------------------------------------------
xflash_defaults, 'Generic', $
CTR_LIM=ctr_lim, ISWP=iswp, $
CONTOURS=contours_def, VELOCITY=velocity_def, $
NBINS = nbins, HIST_SCALE = max_scale
; save the contours
contours = contours_def
velocity = velocity_def
; we need to save ctr_lim, so other events can access it
pblm_ctr_lim = ctr_lim
orientation = iswp
; ---- activate the options based on the # of dims ----------------------------
widget_control, info.zoomLabel, sensitive = 1
widget_control, info.xmin, sensitive = 1
widget_control, info.xmax, sensitive = 1
case ndim of
1: begin
widget_control, info.ymin, sensitive = 0
widget_control, info.ymax, sensitive = 0
widget_control, info.zmin, sensitive = 0
widget_control, info.zmax, sensitive = 0
widget_control, info.slicePlane, sensitive = 0
widget_control, info.contourOpts, sensitive = 0
widget_control, info.velocityOpts, sensitive = 0
widget_control, info.particleOpts, sensitive = 0
widget_control, info.histogramOpts, sensitive = 1
end
2: begin
widget_control, info.ymin, sensitive = 1
widget_control, info.ymax, sensitive = 1
widget_control, info.zmin, sensitive = 0
widget_control, info.zmax, sensitive = 0
widget_control, info.slicePlane, sensitive = 0
widget_control, info.contourOpts, sensitive = 1
widget_control, info.velocityOpts, sensitive = 1
widget_control, info.histogramOpts, sensitive = 1
numParticles = get_particle_number(filename)
if (numParticles GT 0) then widget_control, $
info.particleOpts, sensitive = 1
end
3: begin
widget_control, info.ymin, sensitive = 1
widget_control, info.ymax, sensitive = 1
; by default, we are in the x-y plane, so only one z value is needed
widget_control, info.zmin, sensitive = 1
widget_control, info.zmax, sensitive = 0
widget_control, info.slicePlane, sensitive = 1
widget_control, info.contourOpts, sensitive = 0
widget_control, info.velocityOpts, sensitive = 1
widget_control, info.particleOpts, sensitive = 0
widget_control, info.histogramOpts, sensitive = 1
end
endcase
endif
; ---- activate the plot button -----------------------------------------------
widget_control, info.plot, sensitive = 1
widget_control, info.histogram, sensitive = 1
; ---- store the new value of info.vars ---------------------------------------
widget_control, ev.top, set_uvalue = info, /no_copy
endelse
filename_old = filename
endif else begin
filename = filename_old
endelse
end
'finfo': begin
xfile_info, filename
end
;-----------------------------------------------------------------------------
; default options
;-----------------------------------------------------------------------------
'def_item': begin
widget_control, ev.id, Get_Value=defaultProblemName
; to fake the appearance of a checkmark next to the selected menu
; item, the first two characters of the problem name are spaces, or
; an x followed by a space. Strip these off when calling
; xflash_defaults
; load the defaults for the selected problem
xflash_defaults, strmid(defaultProblemName,2), $
CTR_LIM=ctr_lim, DEFAULT_NAMES=problems, ISWP=iswp, $
CONTOURS=contours_def, VELOCITY=velocity_def
contours = contours_def
velocity = velocity_def
; remove the x off the old problem name and put it on the new one
widget_control, current_default_id, GET_VALUE=oldName
widget_control, current_default_id, SET_VALUE=' ' + strmid(oldName,2)
widget_control, ev.id, SET_VALUE='x ' + strmid(defaultProblemName,2)
; store the index
current_default_id = ev.id
current_default_name = strmid(defaultProblemName,2)
; store the data ranges globally
pblm_ctr_lim = ctr_lim
orientation = iswp
; get the variable number
ivar = widget_info(info.vars, /droplist_select)
; update the data ranges
widget_control, info.varMin, set_value = pblm_ctr_lim[ivar,0]
widget_control, info.varMax, set_value = pblm_ctr_lim[ivar,1]
; reset the zoom options
widget_control, info.xmin, set_value = -1
widget_control, info.xmax, set_value = -1
widget_control, info.ymin, set_value = -1
widget_control, info.ymax, set_value = -1
widget_control, info.zmin, set_value = -1
widget_control, info.zmax, set_value = -1
end
;-----------------------------------------------------------------------------
; color options
;-----------------------------------------------------------------------------
'colorItem': begin
widget_control, ev.id, GET_VALUE=colorName
; I could not figure out how to have IDL put a checkmark next to the
; color name that is currently selected off the menu, so we fake this
; below
; remove the 'x' off the old color
widget_control, current_color_id, GET_VALUE=oldName
widget_control, current_color_id, SET_VALUE=' ' + strmid(oldName,2)
; put an 'x' on the new one to indicate that it is selected
widget_control, ev.id, SET_VALUE = 'x ' + strmid(colorName,2)
; store the new id and name so we can use them next time
current_color_id = ev.id
current_color_name = strmid(colorName,2)
end
;-----------------------------------------------------------------------------
; xy count options
;-----------------------------------------------------------------------------
'xycountItem': begin
widget_control, ev.id, GET_VALUE=numberName
widget_control, current_number_id, GET_VALUE=oldName
widget_control, current_number_id, SET_VALUE=' ' + strmid(oldName,2)
widget_control, ev.id, SET_VALUE = 'x ' + strmid(numberName,2)
current_number_id = ev.id
current_number_name = strmid(numberName,2)
end
;-----------------------------------------------------------------------------
; variable changed
;-----------------------------------------------------------------------------
'var': begin
; get the new variable number
ivar = widget_info(info.vars, /droplist_select)
; put the new contour limits in
widget_control, info.varMin, set_value = pblm_ctr_lim[ivar,0]
widget_control, info.varMax, set_value = pblm_ctr_lim[ivar,1]
end
;-----------------------------------------------------------------------------
; output changed
;-----------------------------------------------------------------------------
'output': begin
widget_control, info.out_buttons, get_value = ipost
if ipost EQ 1 then begin
widget_control, info.hsize, sensitive = 0
widget_control, info.vsize, sensitive = 0
endif else begin
widget_control, info.hsize, sensitive = 1
widget_control, info.vsize, sensitive = 1
endelse
end
;-----------------------------------------------------------------------------
; automatic data range limits selected
;-----------------------------------------------------------------------------
'auto': begin
; get the value of the auto button
widget_control, info.autoData, get_value = iauto
if iauto[0] EQ 1 then begin
widget_control, info.varMin, sensitive=0
widget_control, info.varMax, sensitive=0
endif else begin
widget_control, info.varMin, sensitive=1
widget_control, info.varMax, sensitive=1
endelse
end
;-----------------------------------------------------------------------------
; slice plane changed (3d)
;-----------------------------------------------------------------------------
'slice': begin
widget_control, info.slicePlane, get_value = islice_dir
case islice_dir of
0: begin
widget_control, info.xmax, sensitive=1
widget_control, info.ymax, sensitive=1
widget_control, info.zmax, sensitive=0
end
1: begin
widget_control, info.xmax, sensitive=1
widget_control, info.ymax, sensitive=0
widget_control, info.zmax, sensitive=1
end
2: begin
widget_control, info.xmax, sensitive=0
widget_control, info.ymax, sensitive=1
widget_control, info.zmax, sensitive=1
end
endcase
end
;-----------------------------------------------------------------------------
; contour options pressed
;-----------------------------------------------------------------------------
'ctropt': begin
xcontour, VARIABLES=varnames
end
;-----------------------------------------------------------------------------
; velocity options pressed
;-----------------------------------------------------------------------------
'velopt': begin
xvelocity
end
;-----------------------------------------------------------------------------
; histogram options pressed
;-----------------------------------------------------------------------------
'histopt': begin
xhist
end
;-----------------------------------------------------------------------------
; plot pressed
;-----------------------------------------------------------------------------
'plot': begin
; file
widget_control, info.fstartSuffix, get_value = startSuffix
startSuffix = fix(startSuffix[0])
widget_control, info.fendSuffix, get_value = endSuffix
endSuffix = fix(endSuffix[0])
if (endSuffix LT startSuffix) then endSuffix = startSuffix
widget_control, info.fstep, get_value = step
step = fix(step[0])
if (step EQ 0) then step = 1
fileInfo = {prototype:filename, $
startSuffix:startSuffix, $
endSuffix:endSuffix, $
step:step}
; figure out if we need to read the data once again. If only one
; file is specified by the fileInfo structure (i.e. startSuffix =
; endSuffix), and it is the same as the last file read. then we
; don't need to read the data again.
baseNmEnd = strpos(fileInfo.prototype, '_', /REVERSE_SEARCH)
filenameBase = strcompress(strmid(fileInfo.prototype, 0, baseNmEnd+1))
filename_read = filenameBase + $
string(fileInfo.startSuffix, format = '(i4.4)')
print, '^^^^ old filename = ', filename_last_read
print, '^^^^ new filename = ', filename_read
; if we are doing a loop over several files, or if we did a loop last
; time, reread the file this time through.
if (fileInfo.startSuffix EQ fileInfo.endSuffix AND $
old_start_sfx EQ old_end_sfx) then begin
; create the filename that we are examining
if (filename_read EQ filename_last_read) then begin
iread = 0
print, '^^^^ dont need to read again'
endif else begin
iread = 1
print, '^^^^ need to read again'
endelse
endif else begin
iread = 1
print, '^^^^ suffixes dont match -- need to read again'
endelse
filename_last_read = filename_read
old_start_sfx = fileInfo.startSuffix
old_end_sfx = fileInfo.endSuffix
; output
widget_control, info.out_buttons, get_value = outputType
widget_control, info.vsize, get_value = vpixels
widget_control, info.hsize, get_value = hpixels
; create a structure to hold the output information
output = {type:outputType, $ ; flag for output medium
hsize:float(hpixels[0]), $ ; # of pixels for graphic in x
vsize:float(vpixels[0])} ; # of pixels for graphic in y
; variable
ivar = widget_info(info.vars, /droplist_select)
; plotting options
widget_control, info.opt_buttons, get_value = opts
; data range
widget_control, info.varMin, get_value = dataMin
widget_control, info.varMax, get_value = dataMax
widget_control, info.autoData, get_value = auto
pblm_ctr_lim[ivar,0] = dataMin
pblm_ctr_lim[ivar,1] = dataMax
variable = {name:varnames(ivar), $
min:dataMin, $
max:dataMax, $
auto:auto[0]}
; create a structure to hold the options
options = {blocks:opts[3], $ ; draw block boundaries
log:opts[0], $ ; take the log of data
annotate:opts[4], $ ; plot time and credits
colorbar:opts[5], $ ; plot colorbar
max:opts[2], $ ; plot max of variable
abs:opts[1], $ ; take abs value of var
tick:opts[6]} ; display tick marks
; zoom options
widget_control, info.slicePlane, get_value = plane
widget_control, info.xmin, get_value = xmin
widget_control, info.xmax, get_value = xmax
widget_control, info.ymin, get_value = ymin
widget_control, info.ymax, get_value = ymax
widget_control, info.zmin, get_value = zmin
widget_control, info.zmax, get_value = zmax
; create a structure to hold the zoom information
zoom = {plane:plane, $
xmin:xmin, $
xmax:xmax, $
ymin:ymin, $
ymax:ymax, $
zmin:zmin, $
zmax:zmax}
particle = {enabled:0, $
plot_vel:1, $
sym_size: 1.0, $
show_tag: 1, $
typical_velocity:10.}
problemInfo = {orientation:orientation, $
name:current_default_name}
if (ndim EQ 1) then begin
xplot1d_amr_new, FILE_INFO = fileInfo, $
VARIABLE_INFO = variable, $
OPTIONS = options, $
OUTPUT = output, $
ZOOM = zoom, $
COLORMAP = current_color_name, $
WIDGET_INFORMATION = info, $
PROBLEM_INFO = problemInfo, $
KNOWN_VARIABLES = varnames, $
TREE_PTR = tree_ptr, $
PARAMS_PTR = params_ptr, $
DATA_PTR = data_ptr, $
PLOT_COUNT = current_number_name
if (current_number_name EQ '1') then begin
widget_control, info.query, SENSITIVE=1
endif else begin
widget_control, info.query, SENSITIVE=0
endelse
widget_control, info.slice1d, SENSITIVE=0
endif else if (ndim EQ 2) then begin
xplot_amr_new, FILE_INFO = fileInfo, $
READ = iread, $
VARIABLE_INFO = variable, $
OPTIONS = options, $
OUTPUT = output, $
VELOCITY = velocity, $
ZOOM = zoom, $
COLORMAP = current_color_name, $
WIDGET_INFORMATION = info, $
CONTOUR_OPT = contours, $
PARTICLE_OPT = particle, $
PROBLEM_INFO = problemInfo, $
KNOWN_VARIABLES = varnames, $
TREE_PTR = tree_ptr, $
PARAMS_PTR = params_ptr, $
DATA_PTR = data_ptr, $
PLOT_COUNT = current_number_name
if (current_number_name EQ '1') then begin
widget_control, info.query, SENSITIVE=1
widget_control, info.slice1d, SENSITIVE=1
endif else begin
widget_control, info.query, SENSITIVE=0
widget_control, info.slice1d, SENSITIVE=0
endelse
endif else if (ndim EQ 3) then begin
xplot3d_amr, FILE_INFO = fileInfo, $
READ = iread, $
VARIABLE_INFO = variable, $
OPTIONS = options, $
OUTPUT = output, $
VELOCITY = velocity, $
CONTOUR_OPT = contours, $
PROBLEM_INFO = problemInfo, $
SLICE_DIR = plane, $
ZOOM = zoom, $
TREE_PTR = tree_ptr, $
PARAMS_PTR = params_ptr, $
DATA_PTR = data_ptr, $
SHOW_DATA = data2d, $
XMERGE = xmerge, $
YMERGE = ymerge, $
ZMERGE = zmerge, $
COLORMAP = current_color_name, $
WIDGET_INFORMATION = info
if (current_number_name EQ '1') then begin
widget_control, info.query, SENSITIVE=0
widget_control, info.slice1d, SENSITIVE=1
endif else begin
widget_control, info.query, SENSITIVE=0
widget_control, info.slice1d, SENSITIVE=0
endelse
endif
end
;-----------------------------------------------------------------------------
; histogram pressed
;-----------------------------------------------------------------------------
'histogram': begin
; file
widget_control, info.fstartSuffix, get_value = startSuffix
startSuffix = fix(startSuffix[0])
widget_control, info.fendSuffix, get_value = endSuffix
endSuffix = fix(endSuffix[0])
if (endSuffix LT startSuffix) then endSuffix = startSuffix
widget_control, info.fstep, get_value = step
step = fix(step[0])
if (step EQ 0) then step = 1
fileInfo = {prototype:filename, $
startSuffix:startSuffix, $
endSuffix:endSuffix, $
step:step}
; output
widget_control, info.out_buttons, get_value = outputType
widget_control, info.vsize, get_value = vpixels
widget_control, info.hsize, get_value = hpixels
; create a structure to hold the output information
output = {type:outputType, $ ; flag for output medium
hsize:float(hpixels[0]), $ ; # of pixels for graphic in x
vsize:float(vpixels[0])} ; # of pixels for graphic in y
; variable
ivar = widget_info(info.vars, /droplist_select)
; data range
widget_control, info.varMin, get_value = dataMin
widget_control, info.varMax, get_value = dataMax
widget_control, info.autoData, get_value = auto
pblm_ctr_lim[ivar,0] = dataMin
pblm_ctr_lim[ivar,1] = dataMax
variable = {name:varnames(ivar), $
min:dataMin[0], $
max:dataMax[0], $
auto:auto[0]}
widget_control, info.opt_buttons, get_value = opts
options = {blocks:opts[3], $ ; draw block boundaries
log:opts[0], $ ; take the log of data
annotate:opts[4], $ ; plot time and credits
colorbar:opts[5], $ ; plot colorbar
max:opts[2], $ ; plot max of variable
abs:opts[1], $ ; take abs value of var
tick:opts[6]} ; display tick marks
hist_driver, FILE_INFO = fileInfo, $
VARIABLE_INFO = variable, $
OPTIONS = options, $
NBINS = nbins, $
HIST_SCALE = max_scale, $
OUTPUT = output, $
KNOWN_VARIABLES = varnames
end
;-----------------------------------------------------------------------------
; query pressed
;-----------------------------------------------------------------------------
'query': begin
widget_control, info.status, $
set_value = 'click in the domain to get info'
case orientation of
0: cursor, x_query, y_query
1: cursor, y_query, x_query
2: cursor, x_query, y_query
endcase
query, X_CURS = x_query, Y_CURS = y_query, WIDGET_INFORMATION = info, $
TREE_PTR = tree_ptr, DATA_PTR = data_ptr, PARAMS_PTR = params_ptr, $
ORIENTATION = orientation, VAR_NAMES = varnames
widget_control, info.status, $
set_value = 'status: awaiting orders'
end
;-----------------------------------------------------------------------------
; slice pressed
;-----------------------------------------------------------------------------
'1dslice': begin
widget_control, info.status, $
set_value = 'click for 1-d slice: left (vertical); right (horizontal)'
print, 'going to get the cursor position'
widget_control, info.slicePlane, get_value = plane
iout = 0
case orientation of
0: begin
cursor, x_query, y_query
button = !MOUSE.BUTTON
if (button EQ 1) then begin
dir = 1
endif else if (button EQ 4) then begin
dir = 0
endif else begin
widget_control, info.status, $
set_value = 'Error: invalid direction'
iout = 1
endelse
end
1: begin
print, 'in 1'
cursor, y_query, x_query, /DOWN, /DATA
print, 'x,y = ', x_query, y_query
button = !MOUSE.BUTTON
if (button EQ 1) then begin
dir = 0
endif else if (button EQ 4) then begin
dir = 1
endif else begin
widget_control, info.status, $
set_value = 'Error: invalid direction'
iout = 1
endelse
end
2: begin
cursor, x_query, y_query
button = !MOUSE.BUTTON
if (button EQ 1) then begin
dir = 1
endif else if (button EQ 4) then begin
dir = 0
endif else begin
widget_control, info.status, $
set_value = 'Error: invalid direction'
iout = 1
endelse
end
endcase
; make sure we are in the domain
; get the limits of the current plot domain
plt_x_min = !x.crange[0]
plt_x_max = !x.crange[1]
plt_y_min = !y.crange[0]
plt_y_max = !y.crange[1]
if ((orientation EQ 0 AND ((x_query LT plt_x_min) OR $
(x_query GT plt_x_max) OR $
(y_query LT plt_y_min) OR $
(y_query GT plt_y_max))) OR $
(orientation EQ 1 AND ((y_query LT plt_x_min) OR $
(y_query GT plt_x_max) OR $
(x_query LT plt_y_min) OR $
(x_query GT plt_y_max))) OR $
(orientation EQ 2 AND ((x_query LT plt_x_min) OR $
(x_query GT plt_x_max) OR $
(y_query GT plt_y_min) OR $
(y_query LT plt_y_max)))) $
then begin
widget_control, info.status, $
set_value = 'status: outside of domain'
iout = 1
endif
if (iout EQ 0) then begin
ivar = widget_info(info.vars, /droplist_select)
print, 'plot a line through ', x_query, y_query
print, 'about to extract_line', dir, x_query, y_query
slice1d = extract3d_line(data2d, $
POINT = [x_query,y_query], $
XMERGE = xmerge, $
YMERGE = ymerge, $
ZMERGE = zmerge, $
SLICE_DIR = plane, $
DIRECTION = dir, $
COORDS = coords1d)
; for i = 0, (size(coords1d))[1]-1, 10 do begin
; print, coords1d[i], slice1d[i]
; endfor
; get the plotting options -- log is the 0th one
widget_control, info.opt_buttons, get_value = opts
log = opts[0]
print, 'current window = ', !D.WINDOW
; save the setting for the master plot window (0)
p0 = !P & x0 = !X & y0 = !Y
window, 1
print, 'new current window = ', !D.WINDOW
; make some plot limits based on this slice
if (log EQ 0) then begin
data1dMin = 0.8*min(slice1d)
data1dMax = 1.2*max(slice1d)
endif else begin
data1dMin = alog10(0.8*min(slice1d))
data1dMax = alog10(1.2*max(slice1d))
endelse
; TO DO: consider orientation !!!!!!!!!!
case plane of
0: begin ; x-y plane
sx = x_query
sy = y_query
sz = zmerge[0]
case dir of
0: sdir = "x"
1: sdir = "y"
endcase
end
1: begin ; x-z plane
sx = x_query
sy = ymerge[0]
sz = y_query
case dir of
0: sdir = "x"
1: sdir = "z"
endcase
end
2: begin ; y-z plane
sx = xmerge[0]
sy = x_query
sz = y_query
case dir of
0: sdir = "x"
1: sdir = "z"
endcase
end
endcase
title = sdir + " slice through (x,y,z) = (" $
+ strcompress(string(sx, FORMAT='(g12.5)')) + ", " $
+ strcompress(string(sy, FORMAT='(g12.5)')) + ", " $
+ strcompress(string(sz, FORMAT='(g12.5)')) + ")"
if (log EQ 0) then begin
plot, coords1d, slice1d, $
YRANGE = [data1dMin,data1dMax], $
YSTYLE = 1, XSTYLE = 1, $
XTITLE = sdir + ' [cm]', $
YTITLE = varnames(ivar), $
TITLE = title
endif else begin
plot, coords1d, alog10(slice1d), $
YRANGE = [data1dMin,data1dMax], $
YSTYLE = 1, XSTYLE = 1, $
XTITLE = sdir + ' [cm]', $
YTITLE = "log("+varnames(ivar)+")", $
TITLE = title
endelse
; switch back to the master plot window
wset, 0
!P = p0 & !X = x0 & !Y = y0
widget_control, info.status, $
set_value = 'status: awaiting orders'
endif
end
;-----------------------------------------------------------------------------
; exit pressed
;-----------------------------------------------------------------------------
'fexit': widget_control, ev.top, /destroy
else:
endcase
end
;==============================================================================
; create the widget that handles 1, 2, and 3-d FLASH data.
;
; Once a prototype file is selected, the main information is read from
; this file to determine the variables stored, existence of particles,
; dimensionality, ... This information will be used to make different
; parts of the main widget sensitive
;
;==============================================================================
pro xflash3d
; common all the controls into the widget
common variables, varnames, native_num
common save_color_info, current_color_id, current_color_name
common save_default_info, current_default_id, current_default_name
common save_number_info, current_number_id, current_number_name
common save_path, current_path
common save_filename, filename_old, filename_last_read, old_start_sfx, $
old_end_sfx
print, 'initializing xflash'
version = !VERSION.RELEASE
print, '... IDL version = ', version
if (version LE 5.3) then begin
print, '... GIF graphics will be used'
endif else begin
print, '... PNG graphics will be used'
endelse
; get the xflash directory from the xflash_dir environmental variable
xflash_dir = getenv('XFLASH_DIR')
; make sure the path ends with a `/'
if strlen(xflash_dir) NE 0 then begin
if strmid(xflash_dir,strlen(xflash_dir)-1,1) NE '/' then $
xflash_dir = xflash_dir + '/'
endif
if (strlen(xflash_dir) NE 0) then begin
print, '... xflash directory is', xflash_dir
endif else begin
print, '... xflash directory is not found. Please set the'
print, ' XFLASH_DIR environment variable'
endelse
spawn, 'pwd', current_path
current_path = current_path[0]
filename_old = ' '
filename_last_read = ' '
old_start_sfx = -99
old_end_sfx = -99
; create the main widget base
mainBase = widget_base(/column, title = 'xflash', MBAR=bar)
; divide this base into vertical sections
fileBase = widget_base(mainBase, /column, /frame)
outputBase = widget_base(mainBase, /column, /frame)
varBase = widget_base(mainBase, /row, /frame)
optBase = widget_base(mainBase, /column, /frame)
rangeBase = widget_base(mainBase, /row, /frame)
zoomBase = widget_base(mainBase, /column, /frame)
advancedBase = widget_base(mainBase, /row, /frame)
endBase = widget_base(mainBase, /row)
statusBase = widget_base(mainBase, /row, /frame)
;------------------------------------------------------------------------------
; load the problem specific information and get the problem names
;------------------------------------------------------------------------------
xflash_defaults, 'Generic', $
NUM_DEFAULTS=num_defaults, DEFAULT_NAMES=default_names, $
CONTOURS=contours, NUM_CONTOURS=num_contours, $
VELOCITY=velocity, $
PARTICLE=particle, $
CTR_LIM=ctr_lim
;------------------------------------------------------------------------------
; create the menus
;------------------------------------------------------------------------------
; ---- file menu --------------------------------------------------------------
fileMenu = widget_button(bar, VALUE='File', /MENU)
fopen = widget_button(fileMenu, VALUE='Open prototype...', UVALUE='fopen')
finfo = widget_button(fileMenu, VALUE='Information', UVALUE='finfo', $
SENSITIVE=0)
fexit = widget_button(fileMenu, VALUE='Exit', UVALUE='fexit')
; ---- default menu -----------------------------------------------------------
defaults_menu = widget_button(bar, VALUE='Defaults', /MENU)
def_item = lonarr(num_defaults)
; put a 'x' next to the currently selected problem name. Pad all
; problems with 2 characters
for i = 0, num_defaults-1 do begin
if (i EQ 0) then begin
def_item[i] = widget_button(defaults_menu, $
VALUE='x '+ default_names[i], $
uvalue='def_item')
endif else begin
def_item[i] = widget_button(defaults_menu, $
VALUE=' '+ default_names[i], $
uvalue='def_item')
endelse
endfor
; store the value of the widget that is the current problem
current_default_id = def_item[0]
current_default_name = default_names[0]
; ---- color menu -------------------------------------------------------------
colorMenu = widget_button(bar, VALUE='Colormap', /MENU)
; get the available colors
dummy = color_index('Grayscale', GET_NAMES=colors)
numcolors = (size(colors))[1]
colorItem = lonarr(numcolors)
; the names on the menu will have two extra characters to indicate the
; current selection
for i = 0, numcolors-1 do begin
if (i EQ 0) then begin
colorItem[i] = widget_button(colorMenu, $
VALUE='x ' + colors[i], $
uvalue='colorItem')
endif else begin
colorItem[i] = widget_button(colorMenu, $
VALUE=' ' + colors[i], $
uvalue='colorItem')
endelse
endfor
; store the value of the widget that is the current selected color
current_color_id = colorItem[0]
current_color_name = colors[0]
; ---- number menu ------------------------------------------------------------
numberMenu = widget_button(bar, VALUE='X/Y plot count')
xycounts = ['1', '1x2', '2x1', '2x2', '2x3', '3x2', '3x3']
numcounts = (size(xycounts))[1]
xycountItem = lonarr(numcounts)
for i = 0, numcounts-1 do begin
if (i EQ 0) then begin
xycountItem[i] = widget_button(numberMenu, $
VALUE = 'x ' + xycounts[i], $
uvalue='xycountItem')
endif else begin
xycountItem[i] = widget_button(numberMenu, $
VALUE = ' ' + xycounts[i], $
uvalue='xycountItem')
endelse
endfor
current_number_id = xycountItem[0]
current_number_name = xycounts[0]
;------------------------------------------------------------------------------
; file selection options
;------------------------------------------------------------------------------
fprototype = widget_label(fileBase, /ALIGN_LEFT, $
value='Prototype File: not yet defined', $
/DYNAMIC_RESIZE)
suf_base = widget_base(fileBase, /row)
fstartSuffix = cw_field(suf_base, title = 'suffix: ', value = '1', $
xsize = 4, uvalue = 'stsfx')
fendSuffix = cw_field(suf_base, title = 'to', value = ' ', $
xsize = 4, uvalue = 'endsfx')
fstep = cw_field(suf_base, title = ' step', value = '1', $
xsize = 4, uvalue = 'step')
widget_control, fstartSuffix, sensitive=0
widget_control, fendSuffix, sensitive=0
widget_control, fstep, sensitive=0
;------------------------------------------------------------------------------
; output options
;------------------------------------------------------------------------------
if (!VERSION.RELEASE LE 5.3) then begin
output = ['screen', 'postscript', 'gif']
endif else begin
output = ['screen', 'postscript', 'png']
endelse
out_buttons = cw_bgroup(outputBase, output, column = 3, $
label_left = 'Output: ', $
/exclusive, set_value = 0, uvalue = 'output')
; create a sub-base for the basename and suffix range to plot
size_base = widget_base(outputBase, /row)
hsize = cw_field(size_base, title = 'Plot size, horizontal: ', $
value = 800, xsize = 5, uvalue = 'hsize')
vsize = cw_field(size_base, title = ' vertical: ', value = 600, $
xsize = 5, uvalue = 'vsize')
widget_control, out_buttons, sensitive=0
widget_control, hsize, sensitive=0
widget_control, vsize, sensitive=0
;------------------------------------------------------------------------------
; variables
;------------------------------------------------------------------------------
varnames = ['----']
vars = widget_droplist(varBase, title = 'Variables: ', uvalue = 'var', $
value = varnames)
widget_control, vars, sensitive = 0
;------------------------------------------------------------------------------
; resolution
;------------------------------------------------------------------------------
resolution = widget_droplist(varBase, title = 'Resolution: ', uvalue = 'res', $
value = ['0'])
widget_control, resolution, sensitive = 0
;------------------------------------------------------------------------------
; data range
;------------------------------------------------------------------------------
varMin = cw_field(rangeBase, title = 'Data range: ', $
xsize = 12, uvalue = 'varMin', value = ctr_lim[0,0])
varMax = cw_field(rangeBase, title = 'to: ', xsize = 12, uvalue = 'varMax', $
value = (ctr_lim[0,1]))
widget_control, varMin, sensitive=0
widget_control, varMax, sensitive=0
dataOpt = ['auto']
autoData = cw_bgroup(rangeBase, dataOpt, $
/nonexclusive, $
set_value = [0], uvalue = 'auto')
widget_control, autoData, sensitive = 0
;------------------------------------------------------------------------------
; general plot options
;------------------------------------------------------------------------------
options = ['log ', $
'abs. value', $
'max', $
'show blocks', $
'annotate', $
'colorbar', $
'show ticks']
opt_buttons = cw_bgroup(optBase, options, column = 4, $
label_left = 'Options: ', $
/nonexclusive, set_value = [1, 0, 0, 0, 1, 1, 1], $
uvalue = 'opt')
widget_control, opt_buttons, sensitive = 0
;------------------------------------------------------------------------------
; zoom
;------------------------------------------------------------------------------
slice_options = ['x-y', 'x-z', 'y-z']
slicePlane = cw_bgroup(zoomBase, slice_options, column = 3, $
label_left = 'Slice Plane: ', $
/exclusive, set_value = 0, uvalue = 'slice')
zoomLabel = widget_label(zoomBase, /align_left, value = 'Zoom:' + $
' (set = -1 for default)')
xbase = widget_base(zoomBase, /row)
ybase = widget_base(zoomBase, /row)
zbase = widget_base(zoomBase, /row)
xmin = cw_field(xbase, title = ' xrange: ', xsize = 9, $
uvalue = 'xmin', value=-1)
xmax = cw_field(xbase, title = 'to ', xsize = 9, $
uvalue = 'xmax', value=-1)
ymin = cw_field(ybase, title = ' yrange: ', xsize = 9, $
uvalue = 'ymin', value=-1)
ymax = cw_field(ybase, title = 'to ', xsize = 9, $
uvalue = 'ymax', value=-1)
zmin = cw_field(zbase, title = ' zrange: ', xsize = 9, $
uvalue = 'zmin', value=-1)
zmax = cw_field(zbase, title = 'to ', xsize = 9, $
uvalue = 'zmax', value=-1)
widget_control, xmin, sensitive = 0
widget_control, xmax, sensitive = 0
widget_control, ymin, sensitive = 0
widget_control, ymax, sensitive = 0
widget_control, zmin, sensitive = 0
widget_control, zmax, sensitive = 0
widget_control, slicePlane, sensitive = 0
widget_control, zoomLabel, sensitive = 0
;------------------------------------------------------------------------------
; advanced options
;------------------------------------------------------------------------------
contourOpts = widget_button(advancedBase, value = 'Contour Options', $
uvalue = 'ctropt')
widget_control, contourOpts, sensitive = 0
velocityOpts = widget_button(advancedBase, value = 'Velocity Options', $
uvalue = 'velopt')
widget_control, velocityOpts, sensitive = 0
particleOpts = widget_button(advancedBase, value = 'Particle Options', $
uvalue = 'partopt')
widget_control, particleOpts, sensitive = 0
histogramOpts = widget_button(advancedBase, value = 'Histogram Options', $
uvalue = 'histopt')
widget_control, histogramOpts, sensitive = 0
;------------------------------------------------------------------------------
; plot buttons and status
;------------------------------------------------------------------------------
plot = widget_button(endBase, value = 'Plot', uvalue = 'plot')
histogram = widget_button(endBase, value = 'Histogram', uvalue = 'histogram')
query = widget_button(endBase, value = 'Query', uvalue = 'query')
slice1d = widget_button(endBase, value = '1-d Slice', uvalue = '1dslice')
widget_control, query, sensitive = 0
widget_control, slice1d, sensitive = 0
widget_control, plot, sensitive = 0
widget_control, histogram, sensitive = 0
status = widget_label(statusBase, value = 'status: ', $
uvalue = 'status', /align_left, /dynamic_resize)
;------------------------------------------------------------------------------
; draw the widget
;------------------------------------------------------------------------------
widget_control, mainBase, /realize
; setup a structure to hold the widget information
info = {mainBase:mainBase, $ ; main base
fileBase:fileBase, $ ; file base widget
fileMenu:fileMenu, $
fopen:fopen, $
finfo:finfo, $
fexit:fexit, $
fprototype:fprototype, $
fstartSuffix:fstartSuffix, $ ; field widget for starting suffix
fendSuffix:fendSuffix, $ ; field widget for ending suffix
fstep:fstep, $ ; field widget for step
colorMenu:colorMenu, $
colorItem:colorItem, $
numberMenu:numberMenu, $
xycountItem:xycountItem, $
outputBase:outputBase, $ ; output base widget
out_buttons:out_buttons, $ ; button group for output
hsize:hsize, $ ; field widget for horizontal image size
vsize:vsize, $ ; field widget for vertical image size
varBase:varBase, $ ; variable base widget
vars:vars, $ ; droplist widget for variable selection
resolution:resolution, $ ; resolution for merge function
optBase:optBase, $ ; option base widget
opt_buttons:opt_buttons, $ ; button group for options
rangeBase:rangeBase, $ ; contour base widget
varMin:varMin, $ ; field widget for min contour level
varMax:varMax, $ ; field widget for max contour level
autoData:autoData, $
advancedBase:advancedBase, $
contourOpts:contourOpts, $ ; contour options button
velocityOpts:velocityOpts, $ ; velocity options button
particleOpts:particleOpts, $ ; particle options button
histogramOpts:histogramOpts, $
zoomBase:zoomBase, $ ; zoom base widget
slicePlane:slicePlane, $
zoomLabel:zoomLabel, $
xmin:xmin, $ ; field widget for min x coord
ymin:ymin, $ ; field widget for min y coord
xmax:xmax, $ ; field widget for max x coord
ymax:ymax, $ ; field widget for max y coord
zmin:zmin, $
zmax:zmax, $
endBase:endBase, $ ; end base widget
plot:plot, $ ; button for plotting
histogram:histogram, $
query:query, $ ; button for querying
slice1d:slice1d, $ ; button for slicing
statusBase:statusBase, $
status:status} ; status bar
; register the info structure with the widget base
widget_control, mainBase, set_uvalue=info, /no_copy
xmanager, 'xflash', mainBase
end
-------------- next part --------------
;-----------------------------------------------------------------------------
; xplot3d_amr.pro -- MZ 3-11-00
;
; plot output from flash -- for 3d data slices
;
; this routine is typically driven by xflash3d
;
; general flow:
; -- read in data
; -- convert data from block structure to 2d slice
; -- scale data to colormap and display
; -- plot block boundaries if desired
; -- plot velocity vectors, if desired
;
;-----------------------------------------------------------------------------
pro xplot3d_amr, FILE_INFO = fileInfo, $
READ = iread, $
VARIABLE_INFO = variable, $
OPTIONS = options, $
VELOCITY = velocity, $
CONTOUR_OPT = contours, $
CORNERS = corners, $
OUTPUT = output, $
PROBLEM_INFO = problem, $
SLICE_DIR = islice_dir, $
ZOOM = zoom, $
TREE_PTR = tree_ptr, $
PARAMS_PTR = params_ptr, $
DATA_PTR = data_ptr, $
SHOW_DATA = show_data, $
XMERGE = xmerge, $
YMERGE = ymerge, $
ZMERGE = zmerge, $
COLORMAP = colormap, $
WIDGET_INFORMATION = info
common system, xflash_dir
; create a common to store the max variable quantity. During the loop
; over the files, it is stored here
common maxvar, max_var
common save_state, params, tree, unk
; load the colortable -- xflash uses special colortables, the first 12
; colors are reserved and constant across the different tables
xflash_dir = get_xflash_path()
; get the colortable
iclrmap = color_index(colormap, MIN_VALUE=colorMin, MAX_VALUE=colorMax)
loadct, iclrmap, FILE = xflash_dir + 'flash_colors.tbl', /SILENT
; Set the long version of the variable name, for the plot title
titleName = variable.name
; override the short name for the standard cases
case variable.name of
'dens': titleName = 'Density (g/cm!E3!N)'
'temp': titleName = 'Temperature (K)'
'gamc': titleName = 'Gamma C'
'game': titleName = 'Gamma E'
'velx': titleName = 'X Velocity (cm/s)'
'vely': titleName = 'Y Velocity (cm/s)'
'velz': titleName = 'Z Velocity (cm/s)'
'pres': titleName = 'Pressure (erg/cm!E3!N)'
'ener': titleName = 'Energy (ergs/g)'
'enuc': titleName = 'Nuclear Energy Generation Rate (ergs/g/s)'
'tot_vel': titleName = 'Total Velocity (cm/s)'
'int_ener': titleName = 'Specific Internal Energy (erg/g)'
'ekin/eint': titleName = 'Kinetic Energy / Internal Energy'
'snd_spd': titleName = 'Sound Speed (cm/s)'
'mach': titleName = 'Mach Number'
else:
endcase
if options.log then titleName = 'Log10 ' + titleName
if options.max then titleName = 'Max ' + titleName
set_plot, 'X'
!p.background = color('white')
!p.color = color('black')
!p.multi = [0,1,1]
if (output.type EQ 1) then begin
!p.charsize = 1.0
endif else begin
if output.hsize GT 1200 then begin
!p.charsize = 3.00
!p.charthick = 2.00
endif else if output.hsize GT 600 then begin
!p.charsize = 1.25
!p.charthick = 1.00
endif else begin
!p.charsize = 0.90
!p.charthick = 1.00
endelse
endelse
; FLASH files end with a 4-digit file number that is incremented as
; files are produced. To get the base, strip the end off of the
; filename. For now, ignore the possiblity of a .gz ending.
fileBase = strmid(fileInfo.prototype,0,strlen(fileInfo.prototype)-4)
pathEnd = strpos(fileBase, '/', /REVERSE_SEARCH)
fileBaseClean = strcompress(strmid(fileBase, pathEnd+1, $
strlen(fileBase)-pathEnd))
;==============================================================================
; loop over the files
;==============================================================================
for ifile = fileInfo.startSuffix, fileInfo.endSuffix, fileInfo.step do begin
; ---- determine the filename -------------------------------------------------
filename = fileBase + string(ifile, format = '(i4.4)')
; later on, we'll figure out whether the file has been read already,
; or if we need to unzip it. For now read the file
itype = determine_file_type(filename)
if (itype EQ -1) then begin
result = dialog_message('Error: file does not exist', $
/error, title = 'ERROR!!')
return
endif
if (n_elements(info) NE 0) then widget_control, info.status, $
set_value = 'status: reading from file'
if iread EQ 1 then begin
case itype of
1: read_amr, filename, TREE=tree, DATA=unk, PARAMETERS=params
2: read_amr_hdf5, filename, TREE=tree, DATA=unk, PARAMETERS=params
endcase
tree_ptr = tree
data_ptr = unk
params_ptr = params
if (n_elements(info) NE 0) then $
widget_control, info.status, $
set_value = 'status: done reading'
endif else begin
print, '^^^^ skipping the read'
endelse
; dump a few statistics on the screen
print, ' number of blocks: ', params.totBlocks
lrefine_max = max(tree[*].lrefine)
top_level = (size(where(tree[*].lrefine EQ lrefine_max)))[1]
print, ' number of top level zones: ', top_level*params.nxb*params.nyb*params.nzb
uniform = params.ntopx*params.nxb*2^(lrefine_max-1)* $
params.ntopy*params.nyb*2^(lrefine_max-1)* $
params.ntopz*params.nzb*2^(lrefine_max-1)
print, ' fraction of uniform grid: ', $
float(top_level*params.nxb*params.nyb*params.nzb)/float(uniform)
if var_index('velx') GE 0 then begin
print, 'min/max(velx)', min(unk(var_index('velx'),*,*,*,*),max=m), m
endif
if var_index('vely') GE 0 then begin
print, 'min/max(vely)', min(unk(var_index('vely'),*,*,*,*),max=m), m
endif
if var_index('velz') GE 0 then begin
print, 'min/max(velz)', min(unk(var_index('velz'),*,*,*,*),max=m), m
endif
if params.time LT 1.e-9 then begin
timeout = 'time = ' + $
string(format = '(f7.3)',params.time*1e12) + ' ps'
endif else if params.time LT 1.e-6 then begin
timeout = 'time = ' + $
string(format = '(f7.3)',params.time*1e9) + ' ns'
endif else if params.time LT 1.e-3 then begin
timeout = 'time = ' + $
string(format = '(f7.3)',params.time*1e6) + ' !4l!3s'
endif else begin
timeout = 'time = ' + $
string(format = '(f7.3)',params.time) + ' s'
endelse
; ---- create the output filename ---------------------------------------------
outfile = fileBaseClean
if options.blocks then outfile = outfile + 'blk_'
outfile = outfile + strcompress(variable.name, /REMOVE_ALL) + $
string(ifile, format = '(i4.4)')
case output.type of
1: outfile = outfile + '.ps'
2: begin
; newer version of IDL no longer support GIFs
if (!VERSION.RELEASE LE 5.3) then begin
outfile = outfile + '.gif'
endif else begin
outfile = outfile + '.png'
endelse
end
else:
endcase
;-----------------------------------------------------------------------------
; load the variable to be plotted into temporary storage
;-----------------------------------------------------------------------------
temp_arr = reform(unk[var_index(variable.name),*,*,*,*])
if options.abs EQ 1 then temp_arr = abs(temp_arr)
print, variable.name + ' extrema: ', min(temp_arr), max(temp_arr)
; if we are automagically scaling, reset the data range
if variable.auto EQ 1 then begin
variable.min = min(temp_arr)
variable.max = max(temp_arr)
endif
;-----------------------------------------------------------------------------
; find the size of the plotting area
;-----------------------------------------------------------------------------
; by default, set the plot extrema to the domain limits. If this was
; overridden by the zoom boxes, then make sure that the value given
; for the limit is within the computational domain.
xmin_tot = min(tree[*].bndBox[0,0])
xmax_tot = max(tree[*].bndBox[1,0])
ymin_tot = min(tree[*].bndBox[0,1])
ymax_tot = max(tree[*].bndBox[1,1])
zmin_tot = min(tree[*].bndBox[0,2])
zmax_tot = max(tree[*].bndBox[1,2])
; right now, a value of -1.0 means use the default. This should be
; changed at some point to allow any valid negative coordinate.
; mimimum coordinates
if (zoom.xmin NE -1.0d0) then begin
xmin = (zoom.xmin > xmin_tot) < max(tree[*].bndBox[0,0])
endif else begin
xmin = xmin_tot
endelse
if (zoom.ymin NE -1.0d0) then begin
ymin = (zoom.ymin > ymin_tot) < max(tree[*].bndBox[0,1])
endif else begin
ymin = ymin_tot
endelse
if (zoom.zmin NE -1.0d0) then begin
zmin = (zoom.zmin > zmin_tot) < max(tree[*].bndBox[0,2])
endif else begin
zmin = zmin_tot
endelse
; maximum coordinates
if (zoom.xmax NE -1.0d0) then begin
xmax = (zoom.xmax < xmax_tot) > xmin*1.00001
endif else begin
xmax = xmax_tot
endelse
if (islice_dir EQ 2) then xmax = xmin
if (zoom.ymax NE -1.0d0) then begin
ymax = (zoom.ymax < ymax_tot) > ymin*1.00001
endif else begin
ymax = ymax_tot
endelse
if (islice_dir EQ 1) then ymax = ymin
if (zoom.zmax NE -1.0d0) then begin
zmax = (zoom.zmax < zmax_tot) > zmin*1.00001
endif else begin
zmax = zmax_tot
endelse
if (islice_dir EQ 0) then zmax = zmin
; compute the aspect ratio, and define the plot size
; determine the window size and
case islice_dir of
0: begin ; x-y plane
dvertical = ymax - ymin
dhorizontal = xmax - xmin
vfrac = dvertical/(ymax_tot - ymin_tot)
hfrac = dhorizontal/(xmax_tot - xmin_tot)
end
1: begin ; x-z plane
dvertical = zmax - zmin
dhorizontal = xmax - xmin
vfrac = dvertical/(zmax_tot - zmin_tot)
hfrac = dhorizontal/(xmax_tot - xmin_tot)
end
2: begin ; y-z plane
dvertical = zmax - zmin
dhorizontal = ymax - ymin
vfrac = dvertical/(zmax_tot - zmin_tot)
hfrac = dhorizontal/(ymax_tot - ymin_tot)
end
endcase
case problem.orientation of
0: dataAspectRatio = dvertical/dhorizontal
1: dataAspectRatio = dhorizontal/dvertical
2: dataAspectRatio = dvertical/dhorizontal
endcase
; ---- initialize the device --------------------------------------------------
if (n_elements(info) NE 0) then $
widget_control, info.status, set_value = 'status: initialize device'
if output.type EQ 1 then begin
current_device = !d.name
set_plot, 'PS'
if dataAspectRatio GE 1. then begin
iorient = 0
endif else begin
iorient = 1
endelse
case iorient of
; portrait orientation
0: begin
xsize = 7.5
ysize = 9.5
device, FILE = outfile, XSIZE = xsize, YSIZE = ysize, $
XOFF = 0.5, YOFF = 0.75, /INCH, /COLOR, $
BITS_PER_PIXEL = 8
end
; landscape orientation
1: begin
xsize = 10.
ysize = 6.
device, FILE = outfile, XSIZE = xsize, YSIZE = ysize, $
XOFF = 1.25, YOFF = 10.5, /INCH, /COLOR, $
BITS_PER_PIXEL = 8, /LANDSCAPE
end
endcase
deviceAspect = ysize/xsize
endif else begin
deviceAspect = output.vsize/output.hsize
endelse
aspect = dataAspectRatio/deviceAspect
; ---- determine the bounds for the plot --------------------------------------
; set the normal coordinates of the portion of the display/page you
; wish to use -- leave a little margin so the plots don't run to the
; edge of the page
page_nx1 = 0.05
page_nx2 = 0.95
page_ny1 = 0.05
page_ny2 = 0.90
; if we are considering multiple plots / page, we'd set some
; boundaries here -- just take the full page for now
dpagex = page_nx2 - page_nx1
dpagey = page_ny2 - page_ny1
nx1 = page_nx1
nx2 = page_nx2
ny1 = page_ny1
ny2 = page_ny2
; aspect compares the aspect ratio of the data to that of the device,
; aspect = (dy_data/dy_device)/(dx_data/dx_device). If aspect > 1,
; it means that the y-coordinate is the one that is going to set the
; overall scaling. If aspect < 1, then the x-coordinate will set the
; scaling. Consider each case separately below.
if (aspect GE 1.) then begin
; the y size of the data sets the scaling of the plot
; set the initial values of the y min and max normal coordinates. We
; leave some room for the axis labels
py1 = ny1
py2 = ny2
; compute the x size, using the aspect ratio
dpy = py2 - py1
dpx = dpy/aspect < .70* dpagex
; recompute dpy, in case dpx was adjusted
dpy = aspect*dpx
; compute the plot coordinates
px1 = nx1
px2 = px1 + dpx
py_center = 0.5*(ny1 + ny2)
py1 = py_center - .5*dpy
py2 = py_center + .5*dpy
; set the plot and legend bounding box -- the legend will be vertical
plot_pos = [px1, py1, px2, py2]
; the legend goes in the remaining space in the x direction -- figure
; out how much room is left there
dx_legend = dpagex + page_nx1 - px2
lwidth = dx_legend/4 < 0.25*(px2 - px1)
lcenter = 0.5* (dpagex + page_nx1 + px2)
legend_pos = [lcenter-0.5*lwidth, py1, lcenter+0.5*lwidth, py2]
endif else begin
; the x size of the data sets the scaling of the plot
; set the initial x min and max normal coordiantes
px1 = nx1
px2 = nx2
dpx = px2 - px1
dpy = aspect*dpx < .7* dpagey
; recompute dpx, in case dpy was adjusted
dpx = dpy/aspect
; recompute the plot coordinates
px_center = 0.5*(nx1 + nx2)
px1 = px_center - .5*dpx
px2 = px_center + .5*dpx
py2 = ny2
py1 = py2 - dpy
; set the plot and legend bounding box -- the legend will be horizontal
plot_pos = [px1, py1, px2, py2]
; the legend goes in the remaining space in the lower y direction --
; figure out how much room is left there
dy_legend = py1 - (page_ny2 - dpagey)
lheight = dy_legend/4 < 0.25*(py2 - py1)
lcenter = 0.5*(py1 + (page_ny2 - dpagey))
legend_pos = [px1, lcenter-0.5*lheight, px2, lcenter+0.5*lheight]
endelse
; create the plot window if necessary
if (ifile EQ fileInfo.startSuffix) then begin
if (output.type NE 1) then $
window, XSIZE = output.hsize, YSIZE = output.vsize, $
TITLE = 'AMR plot'
endif
;-----------------------------------------------------------------------------
; create the axes -- just get the data coordinates defined
;-----------------------------------------------------------------------------
; create the title -- if we are printing it
if (options.annotate EQ 1) then begin
title = titleName
endif else begin
title = ' '
endelse
case islice_dir of
0: begin ; x-y
case problem.orientation of
0: begin
plot, [xmin, xmax], [ymin, ymax], POS = plot_pos, $
XSTYLE = 5, YSTYLE = 5, TITLE = title
oplot, [xmin, xmax], [ymax, ymin]
end
1: begin
plot, [ymin, ymax], [xmin, xmax], POS = plot_pos, $
XSTYLE = 5, YSTYLE = 5, TITLE = title
oplot, [ymin, ymax], [xmax, xmin]
end
2: begin
plot, [xmin, xmax], [ymin, ymax], POS = plot_pos, $
XSTYLE = 5, YSTYLE = 5, YRANGE=[ymax,ymin], TITLE = title
oplot, [xmin, xmax], [ymax, ymin]
end
endcase
end
1: begin ; x-z
case problem.orientation of
0: begin
plot, [xmin, xmax], [zmin, zmax], POS = plot_pos, $
XSTYLE = 5, YSTYLE = 5, TITLE = title
oplot, [xmin, xmax], [zmax, zmin]
end
1: begin
plot, [zmin, zmax], [xmin, xmax], POS = plot_pos, $
XSTYLE = 5, YSTYLE = 5, TITLE = title
oplot, [zmin, zmax], [xmax, xmin]
end
2: begin
plot, [xmin, xmax], [zmin, zmax], POS = plot_pos, $
XSTYLE = 5, YSTYLE = 5, TITLE = title
oplot, [xmin, xmax], [zmax, zmin]
end
endcase
end
2: begin ; y-z
case problem.orientation of
0: begin
plot, [ymin, ymax], [zmin, zmax], POS = plot_pos, $
XSTYLE = 5, YSTYLE = 5, TITLE = title
oplot, [ymin, ymax], [zmax, zmin]
end
1: begin
plot, [zmin, zmax], [ymin, ymax], POS = plot_pos, $
XSTYLE = 5, YSTYLE = 5, TITLE = title
oplot, [zmin, zmax], [ymax, ymin]
end
2: begin
plot, [ymin, ymax], [zmin, zmax], POS = plot_pos, $
XSTYLE = 5, YSTYLE = 5, TITLE = title
oplot, [ymin, ymax], [zmax, zmin]
end
endcase
end
endcase
;------------------------------------------------------------------------------
; determine whether to sub-sample the data
;------------------------------------------------------------------------------
; assume 1024 pixels on screen and 4096 on paper
sample = 0
max_refine = max(tree[*].lrefine)
; compute the number of pixels in a uniform grid of our sub-domain
case islice_dir of
0: max_pixels = params.ntopx*params.nxb*2.^(max_refine-1)*hfrac > $
params.ntopy*params.nyb*2.^(max_refine-1)*vfrac
1: max_pixels = params.ntopx*params.nxb*2.^(max_refine-1)*hfrac > $
params.ntopz*params.nzb*2.^(max_refine-1)*vfrac
2: max_pixels = params.ntopy*params.nyb*2.^(max_refine-1)*hfrac > $
params.ntopz*params.nzb*2.^(max_refine-1)*vfrac
endcase
if output.type EQ 1 then begin
ideal = 4096
if max_pixels GT ideal then $
sample = fix(alog(max_pixels/ideal) / alog(2.)) > 0
endif else begin
ideal = 1024
if max_pixels GT ideal then $
sample = fix(alog(max_pixels/ideal) / alog(2.)) > 0
endelse
; make sure we have atleast one point in every block
if (2^sample GT (params.nxb < params.nyb < params.nzb)) then begin
sample = fix(alog(params.nxb < params.nyb < params.nzb)/alog(2.))
endif
;------------------------------------------------------------------------------
; scale the data
;------------------------------------------------------------------------------
if (n_elements(info) NE 0) then $
widget_control, info.status, SET_VALUE = 'status: scaling data'
; first put it on a uniform grid
help, temp_arr
print, 'ranges:'
print, xmin, xmax
print, ymin, ymax
print, zmin, zmax
; if we are max'ing the variable, treat it specially
if options.max EQ 0 then begin
; set 'sample' to user resolution if less than max. resolution
resolution = widget_info(info.resolution, /droplist_select)
if resolution gt sample then sample = resolution
print, 'resolution = ', sample
; merge it -- corners are handled by merge_amr automagically
nsdata = merge_amr(temp_arr, TREE=tree, PARAMETERS=params, $
XMERGE=x, YMERGE=y, ZMERGE=z, $
SAMPLE = sample, $
XRANGE = [xmin,xmax], $
YRANGE = [ymin,ymax], $
ZRANGE = [zmin,zmax])
nsdata = reform(temporary(nsdata))
show_data = nsdata
xmerge = x
ymerge = y
zmerge = z
if options.log EQ 1 then nsdata = alog10(temporary(nsdata))
case options.log of
0: nsdata = scale_color(temporary(nsdata), $
VARMAX = variable.max, $
VARMIN = variable.min, $
COLORMAP_MIN = colorMin, $
COLORMAP_MAX = colorMax)
1: nsdata = scale_color(temporary(nsdata), $
VARMAX = variable.max, $
VARMIN = variable.min, /log, $
COLORMAP_MIN = colorMin, $
COLORMAP_MAX = colorMax)
endcase
endif else begin
; merge it manually
stempvar = merge_amr(temp_arr, TREE=tree, PARAMETERS=params, $
XMERGE = x, YMERGE = y, ZMERGE = z, $
SAMPLE = sample, $
XRANGE = [xmin,xmax], $
YRANGE = [ymin,ymax], $
ZRANGE = [zmin,zmax])
stempvar = reform(temporary(stempvar))
if options.log EQ 1 then stempvar = alog10(temporary(stempvar))
; take the maximum
if ifile EQ fileInfo.startSuffix then begin
max_var = stempvar
print, 'max = ', max(max_var), min(max_var)
endif else begin
max_var = stempvar > max_var
endelse
undefine, stempvar
; scale it
case options.log of
0: nsdata = scale_color(max_var, $
VARMAX = variable.max, $
VARMIN = variable.min, $
COLORMAP_MIN = colorMin, $
COLORMAP_MAX = colorMax)
1: nsdata = scale_color(max_var, $
VARMAX = variable.max, $
VARMIN = variable.min, /log, $
COLORMAP_MIN = colorMin, $
COLORMAP_MAX = colorMax)
endcase
endelse
if problem.orientation EQ 1 then nsdata = transpose(temporary(nsdata))
;------------------------------------------------------------------------------
; plot the data
;------------------------------------------------------------------------------
if (n_elements(info) NE 0) then $
widget_control, info.status, $
set_value = 'status: plotting data'
; convert the limits of the domain into normal coordinates
case islice_dir of
0: begin
case problem.orientation of
0: begin
lower = convert_coord([xmin, ymin], /DATA, /TO_NORMAL)
upper = convert_coord([xmax, ymax], /DATA, /TO_NORMAL)
end
1: begin
lower = convert_coord([ymin, xmin], /DATA, /TO_NORMAL)
upper = convert_coord([ymax, xmax], /DATA, /TO_NORMAL)
end
2: begin
lower = convert_coord([xmin, ymax], /DATA, /TO_NORMAL)
upper = convert_coord([xmax, ymin], /DATA, /TO_NORMAL)
nsdata = rotate(temporary(nsdata),7)
end
endcase
end
1: begin
case problem.orientation of
0: begin
lower = convert_coord([xmin, ymin], /DATA, /TO_NORMAL)
upper = convert_coord([xmax, ymax], /DATA, /TO_NORMAL)
end
1: begin
lower = convert_coord([zmin, xmin], /DATA, /TO_NORMAL)
upper = convert_coord([zmax, xmax], /DATA, /TO_NORMAL)
end
2: begin
lower = convert_coord([xmin, ymax], /DATA, /TO_NORMAL)
upper = convert_coord([xmax, ymin], /DATA, /TO_NORMAL)
nsdata = rotate(temporary(nsdata),7)
end
endcase
end
2: begin
case problem.orientation of
0: begin
lower = convert_coord([ymin, zmin], /DATA, /TO_NORMAL)
upper = convert_coord([ymax, zmax], /DATA, /TO_NORMAL)
end
1: begin
lower = convert_coord([zmin, ymin], /DATA, /TO_NORMAL)
upper = convert_coord([zmax, ymax], /DATA, /TO_NORMAL)
end
2: begin
lower = convert_coord([ymin, zmax], /DATA, /TO_NORMAL)
upper = convert_coord([ymax, zmin], /DATA, /TO_NORMAL)
nsdata = rotate(temporary(nsdata),7)
end
endcase
end
endcase
sz = size(nsdata)
help, nsdata
tvimage, nsdata, pos = [lower[0], lower[1], upper[0], upper[1]], /overplot
;------------------------------------------------------------------------------
; plot the velocity vectors
;------------------------------------------------------------------------------
if velocity.enabled then begin
if (n_elements(info) NE 0) then $
widget_control, info.status, $
set_value = 'status: plotting velocity vectors'
;
; NOTE, the shown plane is allways named 'xy'-plane and
; the velocity vectors in the plane are always called vx and vy !!
;
case islice_dir of
0: begin ; x-y
if var_index('velx') GE 0 AND var_index('vely') GE 0 then begin
velx = reform(unk(var_index('velx'),*,*,*,*))
vely = reform(unk(var_index('vely'),*,*,*,*))
vx = merge_amr(velx, TREE=tree, PARAMETERS=params, $
XMERGE=x, YMERGE=y, $
SAMPLE = sample, $
XRANGE = [xmin,xmax], $
YRANGE = [ymin,ymax], $
ZRANGE = [zmin,zmax])
vx = reform(temporary(vx))
vy = merge_amr(vely, TREE=tree, PARAMETERS=params, $
SAMPLE = sample, $
XRANGE = [xmin,xmax], $
YRANGE = [ymin,ymax], $
ZRANGE = [zmin,zmax])
vy = reform(temporary(vy))
endif
end
1: begin ; x-z
if var_index('velx') GE 0 AND var_index('velz') GE 0 then begin
velx = reform(unk(var_index('velx'),*,*,*,*))
vely = reform(unk(var_index('velz'),*,*,*,*))
vx = merge_amr(velx, TREE=tree, PARAMETERS=params, $
XMERGE=x, ZMERGE=y, $
SAMPLE = sample, $
XRANGE = [xmin,xmax], $
YRANGE = [ymin,ymax], $
ZRANGE = [zmin,zmax])
vx = reform(temporary(vx))
vy = merge_amr(vely, TREE=tree, PARAMETERS=params, $
SAMPLE = sample, $
XRANGE = [xmin,xmax], $
YRANGE = [ymin,ymax], $
ZRANGE = [zmin,zmax])
vy = reform(temporary(vy))
endif
end
2: begin ; y-z
if var_index('vely') GE 0 AND var_index('velz') GE 0 then begin
velx = reform(unk(var_index('vely'),*,*,*,*))
vely = reform(unk(var_index('velz'),*,*,*,*))
vx = merge_amr(velx, TREE=tree, PARAMETERS=params, $
YMERGE=x, ZMERGE=y, $
SAMPLE = sample, $
XRANGE = [xmin,xmax], $
YRANGE = [ymin,ymax], $
ZRANGE = [zmin,zmax])
vx = reform(temporary(vx))
vy = merge_amr(vely, TREE=tree, PARAMETERS=params, $
SAMPLE = sample, $
XRANGE = [xmin,xmax], $
YRANGE = [ymin,ymax], $
ZRANGE = [zmin,zmax])
vy = reform(temporary(vy))
endif
end
endcase
if (problem.orientation EQ 1) then begin
tmp = transpose(vy)
vy = transpose(vx)
vx = tmp
undefine, tmp
ty = y
y = x
x = ty
undefine, ty
endif
xpts = (size(x))[1]
ypts = (size(y))[1]
xt = x # replicate(1,ypts)
yt = replicate(1,xpts) # y
help, vx, vy, xt, yt
partvelvec, vx, vy, xt, yt, /over, color = color('verygray'), $
minmag = velocity.min_velocity, maxmag = velocity.max_velocity, $
typvel = velocity.typical_velocity, $
xskip = velocity.xskip, yskip = velocity.yskip, legend = [.70,.03]
undefine, velx
undefine, vely
endif
;-----------------
; end of plot velocity vectors
;-----------------
if options.blocks then draw_blocks, color('ltblue'), TREE = tree, PARAMETERS = params, ORIENTATION = problem.orientation
;------------------
; restore the axes
;------------------
case islice_dir of
0: begin
case problem.orientation of
0: begin
if (problem.name EQ 'xray_cyl') then begin
case options.tick of
0: begin ; no tick marks
end
1: begin ; display tick marks
axis, xmin, ymin, xaxis = 0, XSTYLE = 1, $
xtitle = 'r (cm)'
axis, xmin, ymin, yaxis = 0, YSTYLE = 1, $
ytitle = 'z (cm)'
axis, xmin, ymax, xaxis = 1, XSTYLE = 1, $
xtickformat = 'nolabel'
axis, xmax, ymin, yaxis = 1, YSTYLE = 1, $
ytickformat = 'nolabel'
end
endcase
endif else begin
case options.tick of
0: begin ; no tick marks
end
1: begin ; display tick marks
axis, xmin, ymin, xaxis = 0, XSTYLE = 1, $
xtitle = 'x (cm)'
axis, xmin, ymin, yaxis = 0, YSTYLE = 1, $
ytitle = 'y (cm)'
axis, xmin, ymax, xaxis = 1, XSTYLE = 1, $
xtickformat = 'nolabel'
axis, xmax, ymin, yaxis = 1, YSTYLE = 1, $
ytickformat = 'nolabel'
end
endcase
endelse
end
1: begin
case options.tick of
0: begin ; no tick marks
end
1: begin ; display tick marks
axis, ymin, xmin, xaxis = 0, XSTYLE = 1, xtitle = 'y (cm)'
axis, ymin, xmin, yaxis = 0, YSTYLE = 1, ytitle = 'x (cm)'
axis, ymin, xmax, xaxis = 1, XSTYLE = 1, xtickformat = 'nolabel'
axis, ymax, xmin, yaxis = 1, YSTYLE = 1, ytickformat = 'nolabel
end
endcase
end
2: begin
case options.tick of
0: begin ; no tick marks
end
1: begin ; display tick marks
axis, xmin, ymax, xaxis = 0, XSTYLE = 1, xtitle = 'x (cm)'
axis, xmin, ymax, yaxis = 0, YSTYLE = 1, ytitle = 'y (cm)'
axis, xmin, ymin, xaxis = 1, XSTYLE = 1, $
xtickformat = 'nolabel'
axis, xmax, ymax, yaxis = 1, YSTYLE = 1, $
ytickformat = 'nolabel'
end
endcase
end
endcase
end
1: begin
case problem.orientation of
0: begin
case options.tick of
0: begin ; no tick marks
end
1: begin ; display tick marks
axis, xmin, zmin, xaxis = 0, XSTYLE = 1, xtitle = 'x (cm)'
axis, xmin, zmin, yaxis = 0, YSTYLE = 1, ytitle = 'z (cm)'
axis, xmin, zmax, xaxis = 1, XSTYLE = 1, xtickformat = 'nolabel'
axis, xmax, zmin, yaxis = 1, YSTYLE = 1, ytickformat = 'nolabel'
end
endcase
end
1: begin
case options.tick of
0: begin ; no tick marks
end
1: begin ; display tick marks
axis, zmin, xmin, xaxis = 0, XSTYLE = 1, xtitle = 'z (cm)'
axis, zmin, xmin, yaxis = 0, YSTYLE = 1, ytitle = 'x (cm)'
axis, zmin, xmax, xaxis = 1, XSTYLE = 1, xtickformat = 'nolabel'
axis, zmax, xmin, yaxis = 1, YSTYLE = 1, ytickformat = 'nolabel'
end
endcase
end
2: begin
case options.tick of
0: begin ; no tick marks
end
1: begin ; display tick marks
axis, xmin, zmin, xaxis = 0, XSTYLE = 1, xtitle = 'x (cm)'
axis, xmin, zmin, yaxis = 0, YSTYLE = 1, ytitle = 'z (cm)'
axis, xmin, zmax, xaxis = 1, XSTYLE = 1, xtickformat = 'nolabel'
axis, xmax, zmin, yaxis = 1, YSTYLE = 1, ytickformat = 'nolabel'
end
endcase
end
endcase
end
2: begin
case problem.orientation of
0: begin
case options.tick of
0: begin ; no tick marks
end
1: begin ; display tick marks
axis, ymin, zmin, xaxis = 0, XSTYLE = 1, xtitle = 'y (cm)'
axis, ymin, zmin, yaxis = 0, YSTYLE = 1, ytitle = 'z (cm)'
axis, ymin, zmax, xaxis = 1, XSTYLE = 1, xtickformat = 'nolabel'
axis, ymax, zmin, yaxis = 1, YSTYLE = 1, ytickformat = 'nolabel'
end
endcase
end
1: begin
case options.tick of
0: begin ; no tick marks
end
1: begin ; display tick marks
axis, zmin, ymin, xaxis = 0, XSTYLE = 1, xtitle = 'z (cm)'
axis, zmin, ymin, yaxis = 0, YSTYLE = 1, ytitle = 'y (cm)'
axis, zmin, ymax, xaxis = 1, XSTYLE = 1, xtickformat = 'nolabel'
axis, zmax, ymin, yaxis = 1, YSTYLE = 1, ytickformat = 'nolabel'
end
endcase
end
2: begin
case options.tick of
0: begin ; no tick marks
end
1: begin ; display tick marks
axis, ymin, zmin, xaxis = 0, XSTYLE = 1, xtitle = 'y (cm)'
axis, ymin, zmin, yaxis = 0, YSTYLE = 1, ytitle = 'z (cm)'
axis, ymin, zmax, xaxis = 1, XSTYLE = 1, xtickformat = 'nolabel'
axis, ymax, zmin, yaxis = 1, YSTYLE = 1, ytickformat = 'nolabel'
end
endcase
end
endcase
end
endcase
; save the plot window status so we can play with it later
p_hydro = !p & x_hydro = !x & y_hydro = !y
;-----------------------
; add a colorbar legend if desired
;-----------------------
if options.colorbar EQ 1 then begin
if aspect LT 1 then begin
case options.log of
0: colorbar2, float(variable.min), float(variable.max), legend_pos, $
COLORMIN = colorMin, COLORMAX = colorMax, CHARSIZE = 1.5
1: colorbar2, alog10(variable.min), alog10(variable.max), $
legend_pos, COLORMIN = colorMin, COLORMAX = colorMax, $
CHARSIZE = 1.5
endcase
endif else begin
case options.log of
0: vcolorbar, float(variable.min), float(variable.max), legend_pos, $
COLORMIN = colorMin, COLORMAX = colorMax, CHARSIZE = 1.5
1: vcolorbar, alog10(variable.min), alog10(variable.max), $
legend_pos, COLORMIN = colorMin, COLORMAX = colorMax, $
CHARSIZE = 1.5
endcase
endelse
endif
; restore the plot to the hydro plot so we can play with it
!p = p_hydro & !x = x_hydro & !y = y_hydro
; --------------------------------------------------------
; finally add some information to the bottom of the plot
; --------------------------------------------------------
gridtxt = 'number of blocks = ' + string(format = '(i7)', params.totBlocks)
gridtxt2 = 'AMR levels = ' + string(format = '(i5)', max(tree[*].lrefine))
; if desired, print the time, filename, and credits
if options.annotate EQ 1 then begin
if output.type EQ 1 then begin
xyouts, .05, .05, timeout, color = color('black'), /normal
xyouts, .05, .02, gridtxt + ', ' + gridtxt2, $
color = color('black'), /normal
xyouts, .05, -.05, filename, color = color('dkgray'), $
/normal, charsize = .8
endif else begin
xyouts, .05, .07, timeout, color = color('black'), /normal
xyouts, .05, .04, gridtxt, color = color('black'), /normal
xyouts, .05, .01, gridtxt2, color = color('black'), /normal
endelse
endif
; empty the buffer
empty
case output.type of
1: begin
device, /close
set_plot, current_device
end
2: color_bitmap, outfile
else:
endcase
if (n_elements(info) NE 0) then $
widget_control, info.status, $
set_value = 'status: awaiting orders . . .'
endfor
end
-------------- next part --------------
;
; returns a 1d slice through 2d data 2003/02 robi banerjee
;
; TO DO: consider problem.orientation
function extract3d_line, data2d, $
POINT = point, $
XMERGE = x, $
YMERGE = y, $
ZMERGE = z, $
SLICE_DIR = islice_dir, $
DIRECTION = dir, $
PROBLEM_INFO = problem, $
COORDS = coords1d
help, data2d
; default return values
coords1d = 0.0
data1d = 0.0
nx = n_elements(x)
ny = n_elements(y)
nz = n_elements(z)
;
; assume equidistant spacing
;
dxh = abs((x[nx-1]-x[0])/float(nx-1))/2.
dyh = abs((y[ny-1]-y[0])/float(ny-1))/2.
dzh = abs((z[nz-1]-z[0])/float(nz-1))/2.
case islice_dir of
0: begin ; x-y plane
case dir of ; along x direction
0: begin
idx = where( y-dyh le point[1] AND y+dyh gt point[1], rc)
if rc gt 0 then begin
print, 'draw x-line at (y,z) = ',y[idx[0]],z[0]
data1d = data2d(*,idx[0])
coords1d = x
endif
end
1: begin ; along y direction
idx = where( x-dxh le point[0] AND x+dxh gt point[0], rc)
if rc gt 0 then begin
print, 'draw y-line at (x,z) = ',x[idx[0]],z[0]
data1d = data2d(idx[0],*)
coords1d = y
endif
end
endcase
end
1: begin ; x-z plane
case dir of ; along x direction
0: begin
idx = where( z-dzh le point[1] AND z+dzh gt point[1], rc)
if rc gt 0 then begin
print, 'draw x-line at (y,z) = ',y[0],z[idx[0]]
data1d = data2d(*,idx[0])
coords1d = x
endif
end
1: begin ; along z direction
idx = where( x-dxh le point[0] AND x+dxh gt point[0], rc)
if rc gt 0 then begin
print, 'draw z-line at (x,y) = ',x[idx[0]],y[0]
data1d = data2d(idx[0],*)
coords1d = z
endif
end
endcase
end
2: begin ; y-z plane
case dir of ; along y direction
0: begin
idx = where( z-dzh le point[1] AND z+dzh gt point[1], rc)
if rc gt 0 then begin
print, 'draw y-line at (x,z) = ',x[0],z[idx[0]]
data1d = data2d(*,idx[0])
coords1d = y
endif
end
1: begin ; along z direction
idx = where( y-dyh le point[0] AND y+dyh gt point[0], rc)
if rc gt 0 then begin
print, 'draw z-line at (x,y) = ',x[0],y[idx[0]]
data1d = data2d(idx[0],*)
coords1d = z
endif
end
endcase
end
endcase
help, data1d, coords1d
return, data1d
end
-------------- next part --------------
;------------------------------------------------------------------------------
; merge_amr.pro
;
; take a hydro variable in AMR block format and merge it into a
; single, uniformly gridded array at the desired AMR resolution.
;
; this version is for 2-d or 3-d data. The dimensionality of the
; dataset is determined from the PARAMETERS structure that is passed
; in along with the dataset.
;
; corner data is supported automatically, if params.corners is set to
; 1. The data will be interpolated to the cell centers, before
; scaling to a uniform grid.
;
; arguments:
;
; tvar -- the variable to be scaled, tvar[maxblocks,nxb,nyb
;
; TREE -- the tree information structure, as returned from the
; read_amr routine.
;
; PARAMETERS -- the parameter structure, as returned from the
; read_amr routine.
;
; optional arguments
;
; XRANGE -- the minimum and maximum x coord of the desired uniform
; domain
;
; YRANGE -- the minimum and maximum y coord of the desired uniform
; domain
;
; ZRANGE -- the minimum and maximum z coord of the desired uniform
; domain (3-D only).
;
; Note, for 3-d data, if the range in any 1 coordinate direction is
; set equal to each other, then a 2-d plane will be produced,
; passing through that coord point.
;
; Ex: XRANGE=[0.,1.], YRANGE=[0.,1.], ZRANGE=[0.5,0.5]
;
; will merge the data onto a 2-d plane in the x-y plane passing
; through z=0.5.
;
;
; XMERGE \ return the x, y, and z (3-d only) positions of the cells
; YMERGE > in the new master scaled array
; ZMERGE /
;
; sample -- # of levels down from finest resolution to sample
; into the new array. Blocks with finer resolution
; will be averaged. Note: sample = 0 is the same as
; omitting the sample keyword.
;
; DOUBLE -- use double precision -- this requires that the AMR
; data be double precision, usually the result of
; passing /DOUBLE to the read_amr routine.
;
;
; This routine can be used to select a subset of the domain for
; uniform gridding. Only this region is allocated as a uniform box.
; If the minimum and maximum range in a direction are equal, then the
; data is sampled into a plane passing through that coordinate value.
;
; Description of variables:
;
; scaling is the factor that a block is expanded (or shruck) by to put
; it at the same resolution as the scaled domain. Coarsely refined
; blocks need to be expanded greatly when mapping into the uniform domain
;
; {x,y,z}range_min_index is the offset, in units of uniform zones from
; the computational zone boundary, and the boundary of the domain we
; are scaling.
;
; {x,y,z}start are the location into the scaled block to copy. In
; most cases, this will be 0, as the beginning of the block normally
; falls into the scaled domain. If the domain intersects a block,
; then this will be 0 < {x,y,z}start < scaled*n{x,y,z}b - 1
;
; {x,y,z}index is the offset from the block lower boundary to the
; lower boundary of the computation domain
;
;
; +--------------------------------------------------------+
; | |
; | |
; | ................................ |
; | . . |
; | . . |
; | . . |
; | . . |
; | . . |
; | . . |
; | . . |
; | . xxxxxxxx . |
; | . x x . |
; | . x x . |
; | . x x . |
; | . x x . |
; | . xxxxxxxx . |
; | . . |
; | . . |
; | ................................ |
; | |
; | |
; +--------------------------------------------------------+
;
; |<- A ->|
; |<-B ->|
;
;
; | marks the computation domain boundary
; . marks the region we are scaling to a uniform mesh
; x marks the current block
;
; A is xrange_min_index
; B is xindex
;
;
; The case is more complicated when the block goes out of the domain
; we are scaling
;
;
; +--------------------------------------------------------+
; | |
; | |
; | ................................ |
; | . . |
; | . . |
; | . . |
; | . . |
; | . . |
; | . . |
; | . . |
; | xxxxxxxxxxxxxxxxxx . |
; | x . x . |
; | x . x . |
; | x . x . |
; | x . x . |
; | xxxxxxxxxxxxxxxxxx . |
; | . . |
; | . . |
; | ................................ |
; | |
; | |
; +--------------------------------------------------------+
;
; |<- A ->|
; ->| B|<-
;
; now, xindex (B) is negative
; and xstart is -B -- we start in the center of the block when copying
; data into the uniform grid
;
;
;------------------------------------------------------------------------------
function merge_amr, temp_arr, $
TREE = tree, PARAMETERS = params, $
DOUBLE = double, $
XRANGE = txrange, YRANGE = tyrange, ZRANGE = tzrange, $
XMERGE = xout, YMERGE = yout, ZMERGE = zout, $
SAMPLE = sample
; find the minimum and maximum refinement levels of the good data
index_good = where(tree[*].nodetype EQ 1)
lmax = max(tree[index_good].lrefine)
lmin = min(tree[index_good].lrefine)
if n_elements(sample) EQ 0 then begin
lwant = lmax
endif else begin
lwant = lmax - sample
if lwant LE 0 then begin
print, 'Warning cannot sample below coarsest resolution'
lwant = 1
endif
endelse
if n_elements(double) EQ 0 then double = 0
; if no domain ranges are specified, then set them to the extreme of
; the domain
if n_elements(txrange) NE 2 then begin
xrange = fltarr(2)
xrange[0] = min(tree[*].bndBox[0,0]) ; set to minimum of x coord
xrange[1] = max(tree[*].bndBox[1,0]) ; set to maximum of x coord
endif else begin
xrange = float(txrange)
endelse
if n_elements(tyrange) NE 2 then begin
yrange = fltarr(2)
yrange[0] = min(tree[*].bndBox[0,1]) ; set to minimum of y coord
yrange[1] = max(tree[*].bndBox[1,1]) ; set to maximum of y coord
endif else begin
yrange = float(tyrange)
endelse
if (params.ndim EQ 3) then begin
if n_elements(tzrange) NE 2 then begin
zrange = fltarr(2)
zrange[0] = min(tree[*].bndBox[0,2]) ; set to minimum of z coord
zrange[1] = max(tree[*].bndBox[1,2]) ; set to maximum of z coord
endif else begin
zrange = float(tzrange)
endelse
endif
if params.corners EQ 1 then begin
print, 'corner data detected, interpolating to cell centers'
nxb_centered = params.nxb - 1
nyb_centered = params.nyb - 1
if (params.ndim EQ 3) then nzb_centered = params.nzb - 1
endif
;------------------------------------------------------------------------------
; compute the number of zones in each direction for the user's domain
; limits, assuming uniform grid
;------------------------------------------------------------------------------
case params.corners of
0: begin ; no corners
dx_fine = (max(tree[*].bndBox[1,0]) - $
min(tree[*].bndBox[0,0]))/ $
(params.ntopx*params.nxb*2^(lwant-1))
if xrange[1] EQ xrange[0] then xrange[1] = xrange[0] + dx_fine
dx_merge = xrange[1] - xrange[0]
dy_fine = (max(tree[*].bndBox[1,1]) - $
min(tree[*].bndBox[0,1]))/ $
(params.ntopy*params.nyb*2^(lwant-1))
if yrange[1] EQ yrange[0] then yrange[1] = yrange[0] + dy_fine
dy_merge = yrange[1] - yrange[0]
if (params.ndim EQ 3) then begin
dz_fine = (max(tree[*].bndBox[1,2]) - $
min(tree[*].bndBox[0,2]))/ $
(params.ntopz*params.nzb*2^(lwant-1))
if zrange[1] EQ zrange[0] then zrange[1] = zrange[0] + dz_fine
dz_merge = zrange[1] - zrange[0]
endif
end
1: begin ; corners
dx_fine = (max(tree[*].bndBox[1,0]) - $
min(tree[*].bndBox[0,0]))/ $
(params.ntopx*nxb_centered*2^(lwant-1))
if xrange[1] EQ xrange[0] then xrange[1] = xrange[0] + dx_fine
dx_merge = xrange[1] - xrange[0]
dy_fine = (max(tree[*].bndBox[1,1]) - $
min(tree[*].bndBox[0,1]))/ $
(params.ntopy*nyb_centered*2^(lwant-1))
if yrange[1] EQ yrange[0] then yrange[1] = yrange[0] + dy_fine
dy_merge = yrange[1] - yrange[0]
if (params.ndim EQ 3) then begin
dz_fine = (max(tree[*].bndBox[1,2]) - $
min(tree[*].bndBox[0,2]))/ $
(params.ntopz*nzb_centered*2^(lwant-1))
if zrange[1] EQ zrange[0] then zrange[1] = zrange[0] + dz_fine
dz_merge = zrange[1] - zrange[0]
endif
end
endcase
nx = fix(dx_merge/dx_fine) > 1
ny = fix(dy_merge/dy_fine) > 1
if params.ndim EQ 3 then nz = fix(dz_merge/dz_fine) > 1
;------------------------------------------------------------------------------
; declare the storage for the uniformly gridded result
;------------------------------------------------------------------------------
case params.ndim of
2: begin
if (double) then begin
temp_merge = dblarr(nx,ny)
endif else begin
temp_merge = fltarr(nx,ny)
endelse
; initialize temp_merge to negative numbers, if there are any holes in
; the domain (i.e. non-rectangular), then these will remain negative
; and scale will color these black
temp_merge[*,*] = -1.e30
end
3: begin
if (double) then begin
temp_merge = dblarr(nx,ny,nz)
endif else begin
temp_merge = fltarr(nx,ny,nz)
endelse
; initialize temp_merge to negative numbers
temp_merge[*,*] = -1.e30
end
endcase
;------------------------------------------------------------------------------
; compute the x, y, and z coordinates of the data for the entire domain
;------------------------------------------------------------------------------
xmax = max(tree[*].bndBox[1,0])
xmin = min(tree[*].bndBox[0,0])
ymax = max(tree[*].bndBox[1,1])
ymin = min(tree[*].bndBox[0,1])
if (params.ndim EQ 3) then begin
zmax = max(tree[*].bndBox[1,2])
zmin = min(tree[*].bndBox[0,2])
endif
case params.corners of
0: begin ; no corners
nx_domain = params.ntopx*params.nxb*2^(lwant-1)
ny_domain = params.ntopy*params.nyb*2^(lwant-1)
if params.ndim EQ 3 then nz_domain = $
params.ntopz*params.nzb*2^(lwant-1)
end
1: begin ; corners
nx_domain = params.ntopx*nxb_centered*2^(lwant-1)
ny_domain = params.ntopy*nyb_centered*2^(lwant-1)
if params.ndim EQ 3 then nz_domain = $
params.ntopz*nzb_centered*2^(lwant-1)
end
endcase
dx = (xmax - xmin)/nx_domain
x = (findgen(nx_domain) + .5)*dx + xmin
dy = (ymax - ymin)/ny_domain
y = (findgen(ny_domain) + .5)*dy + ymin
if (params.ndim EQ 3) then begin
dz = (zmax - zmin)/nz_domain
z = (findgen(nz_domain) + .5)*dz + zmin
endif
;------------------------------------------------------------------------------
; find out where the range extreme would fall in a uniform grid of the
; total domain -- compute these using nx and ny above to avoid errors
; from rounding
;------------------------------------------------------------------------------
xrange_min_index = (where(xrange[0] LE x))[0]
xrange_max_index = xrange_min_index + nx - 1
yrange_min_index = (where(yrange[0] LE y))[0]
yrange_max_index = yrange_min_index + ny - 1
if (params.ndim EQ 3) then begin
zrange_min_index = (where(zrange[0] LE z))[0]
zrange_max_index = zrange_min_index + nz - 1
endif
;------------------------------------------------------------------------------
; loop over the blocks and store them at the resolution of the
; finest block into the master array -- difficulties arise when
; only a part of a block falls into the user's domain subset
;------------------------------------------------------------------------------
i = (lonarr(1))[0]
for i = 1l, params.totBlocks do begin
if tree[i-1].nodetype EQ 1 then begin
; compute the scale factor to bring it to the resolution of the finest block
scaling = 2.^(lwant - tree[i-1].lrefine)
; find out where it should live in the master array
xindex = (where(tree[i-1].bndBox[0,0] LE x))[0] - xrange_min_index
yindex = (where(tree[i-1].bndBox[0,1] LE y))[0] - yrange_min_index
if (params.ndim EQ 3) then $
zindex = (where(tree[i-1].bndBox[0,2] LE z))[0] - zrange_min_index
; set the limits for the local block array
case params.corners of
0: begin ; no corners
xstart = 0
xspan = scaling*params.nxb-1
ystart = 0
yspan = scaling*params.nyb-1
if (params.ndim EQ 3) then begin
zstart = 0
zspan = scaling*params.nzb-1
endif
end
1: begin ; corners
xstart = 0
xspan = scaling*nxb_centered-1
ystart = 0
yspan = scaling*nyb_centered-1
if (params.ndim EQ 3) then begin
zstart = 0
zspan = scaling*nzb_centered-1
endif
end
endcase
; first check if we go over the maximum
if (xindex + xspan GT (xrange_max_index-xrange_min_index)) then $
xspan = (xrange_max_index-xrange_min_index) - xindex
; now check if we go below the minimum
if (xindex LT 0 AND xindex + xspan GE 0) then begin
xstart = -xindex
xspan = xindex+xspan
xindex = 0
endif
if (yindex + yspan GT (yrange_max_index-yrange_min_index)) then $
yspan = (yrange_max_index-yrange_min_index) - yindex
if (yindex LT 0 AND yindex + yspan GE 0) then begin
ystart = -yindex
yspan = yindex+yspan
yindex = 0
endif
if (params.ndim EQ 3) then begin
if (zindex + zspan GT (zrange_max_index-zrange_min_index)) then $
zspan = (zrange_max_index-zrange_min_index) - zindex
if (zindex LT 0 AND zindex + zspan GE 0) then begin
zstart = -zindex
zspan = zindex+zspan
zindex = 0
endif
endif
case params.ndim of
2: begin
if (xindex + xspan GE 0 AND $
yindex + yspan GE 0 AND $
xindex LE (xrange_max_index-xrange_min_index) AND $
yindex LE (yrange_max_index-yrange_min_index)) then begin
; store it
sub_array = congrid(reform(temp_arr[i-1,*,*]), $
scaling*params.nxb, $
scaling*params.nyb)
temp_merge[xindex:xindex+xspan, $
yindex:yindex+yspan] = $
sub_array[xstart:xstart+xspan, $
ystart:ystart+yspan]
endif
end
3: begin
if (xindex + xspan GE 0 AND $
yindex + yspan GE 0 AND $
zindex + zspan GE 0 AND $
xindex LE (xrange_max_index-xrange_min_index) AND $
yindex LE (yrange_max_index-yrange_min_index) AND $
zindex LE (zrange_max_index-zrange_min_index)) then begin
; store it
case params.corners of
0: begin ; no corners
sub_array = rebin(reform(temp_arr[i-1,*,*,*]), $
scaling*params.nxb, $
scaling*params.nyb, $
scaling*params.nzb, /SAMPLE)
end
1: begin ; corners
; create a temporary array to hold the interpolated block's data
case double of
0: var_cc = fltarr(nxb_centered, $
nyb_centered, $
nzb_centered)
1: var_cc = dblarr(nxb_centerd, $
nyb_centerd, $
nzb_centered)
endcase
; get the temp array by averaging to cell centers
for kz = 0, nzb_centered-1 do begin
for jz = 0, nyb_centered-1 do begin
for iz = 0, nxb_centered-1 do begin
var_cc[iz,jz,kz] = .125* $
(temp_arr[i-1,iz, jz, kz ] + $
temp_arr[i-1,iz+1,jz, kz ] + $
temp_arr[i-1,iz, jz+1,kz ] + $
temp_arr[i-1,iz+1,jz+1,kz ] + $
temp_arr[i-1,iz, jz, kz+1] + $
temp_arr[i-1,iz+1,jz, kz+1] + $
temp_arr[i-1,iz, jz+1,kz+1] + $
temp_arr[i-1,iz+1,jz+1,kz+1])
endfor
endfor
endfor
sub_array = rebin(reform(var_cc), $
scaling*nxb_centered, $
scaling*nyb_centered, $
scaling*nzb_centered, /SAMPLE)
end
endcase
temp_merge[xindex:xindex+xspan, $
yindex:yindex+yspan, $
zindex:zindex+zspan] = $
sub_array[xstart:xstart+xspan, $
ystart:ystart+yspan, $
zstart:zstart+zspan]
endif
end
endcase
endif ; nodetype
endfor
help, temp_merge
xout = x[xrange_min_index:xrange_max_index]
yout = y[yrange_min_index:yrange_max_index]
if params.ndim EQ 3 then zout = z[zrange_min_index:zrange_max_index]
return, temp_merge
end
More information about the flash-bugs
mailing list