diff --git a/Makefile b/Makefile
index fc9a1f1..36b1cbc 100644
--- a/Makefile
+++ b/Makefile
@@ -11,7 +11,7 @@ PYTHON = anaconda3
 
 # git directories we should have here
 
-GIT_DIRS = athena athenak
+GIT_DIRS = athena athenak pythena
 
 # URLs that we'll need
 
@@ -117,6 +117,13 @@ build_nemo:    nemo
 ## build_python:    build your private anaconda3
 build_python:  python
 
+## build_pythena:   testing a new pyThena, rebuild from scratch
+build_pythena:   
+	rm -rf pyThena
+	git clone $(URL6)
+	(cd pyThena; make build_athenak; make run1; make test1)
+
+
 # a few sample runs
 
 ##
@@ -237,9 +244,4 @@ agui_a:
 agui_k:
 	git clone $(URL5c) agui_k
 
-# testing a pyThena from scratch:
 
-pythena:
-	rm -rf pyThena
-	git clone $(URL6)
-	(cd pyThena; make build_athenak; make run1; make test1)
diff --git a/TODO.md b/TODO.md
index e427bf9..41084b3 100644
--- a/TODO.md
+++ b/TODO.md
@@ -2,40 +2,77 @@
 
 - python vs. python3 (just a worry if some have no python, but do have python3)
 
-- better detect what's available (tab, vtk, bin)
-  - use the -d flag in plot1.py for directory, it should auto-detect *.tab or tab/*.tab  [done]
-  - use vtk reader to grab vtk files (1d and 2d)  [no more vtk, we use bin]
-  - use bin reader (can do AMR, where vtk cannot) - check between
+- pythana (as package name) is already in use by PyPi:   https://pypi.org/project/pythena/
+  so either we change the name, or live with duplicate.   Plenty of those for NEMO already.
+
+  Also: if we allow 'pip install' we need to build a proper wheel, such that athenak gets
+  pulled and compiled.
+
 
 - for output (run) directory, use basename of the problem set, and then run1, run2, ..... inside of those?
   - the startup tool could look inside the problem run directories
-  - this closes the loop in running and reviewing old runs  
+  - this closes the loop in running and reviewing old runs
+  - could clean up after running?
 
 - GUI peculiarities
-  - when browsing, a cancel does not repopulate with the old entry
+  - when browsing, a cancel does not repopulate with the old entry [fixed]
   - killing parent will not kill the children (could be an advanced option,
     because it's nice killing a parent keeps the children alive.
+  - resizing and sizing of sub-areas in the GUI (avoid negative space)
+  - the delay(ms) in the animator doesn't seem to work well, it's erratic?
+    For example I had it taking 5.8sec to do run2, but that was then
+    independent of the delay setting , which is now 1..20 ms.
+  - in the animation, the picture at T=0 when nothing is loaded yet looks
+  - it's easy to accidentally hut the RUN button twice.... athena will run twice. Do
+    we need to make a lock file?
+
+- GUI time slider:  it's hard to find a specific time, it would be good if the scrollbar
+  and/or click to the left or right of the slider would advance/decrease the slider by
+  one unit. 
+
+- Important need for our GUI tags.   We need a new TAGNAME, that clones a keyword value, e.g.
+       nx=1   # help for a    #>   SLIDER 0:10:1
+       ny=1   # help for b    #>   LINK nx
+  or in athena there are blocks, so there it would need to be something like
+       ny=1   # help for b    #>   LINK mesh/nx
+  this would cause mesh/ny to get the value of mesh/nx.
+
+- GUI ideas
+  - when in animation and changing variables, perhaps stick to the time so one can flip at that time
+    between e.g. x1v-dens and x1v-velx
+  - when in a zoomed mode, remember the new viewing frame and keep that as the animation continues
+    or flipping between line and dots, or flipping between variables?
 
-- option to integrate with a running athena
+- option to integrate with a running athena [advanced]
 
-- plot1d:  bring up multiple plots.  issue will be what scaling to use
+- plot1d:  bring up multiple graphs in one plot.  issue will then how to solve the scaling issue
+
+- plot1d:  allow a fixed scaling [the -f flag does that based on first plot, but not always good enough]
 
 - plot2d:  allow option of contour with color, or just contour or just color. this would
-  (like in plot1d) give us option to combine two variables
+  (like in plot1d) give us option to combine two variables and compare them.
 
+- ensure if with id= we can handle different *.tab files. Probably need a new --id flag in plot1d.py for this
+  In the end, extracting them from the bin file might be faster.
+  
+- save animation as mp4? 
   
 
 # Older vis software from athenak
 
 Awkward to use, but here it is:
 
-      ./plot_tab.py  -i tab/LinWave.hydro_w.00000.tab  -n 100 -v dens
+      ./plot_tab.py  -i run1/tab/LinWave.hydro_w.00000.tab  -n 100 -v dens
 
 and for 2D:
 
       ./plot_slice.py run2/bin/OrszagTang.mhd_w_bcc.00001.bin dens show --notex
 
-but this only brings up one plot, no animation. our plot2d should do that
+but this only brings up one plot, no animation. our plot2d should do that.
+
+And for hst files
+
+      ./plot_hst.py ./plot_hst.py  -i run2/OrszagTang.mhd.hst  -v dt
 
 
 # Features
diff --git a/athena_read.py b/athena_read.py
new file mode 100644
index 0000000..8afb1da
--- /dev/null
+++ b/athena_read.py
@@ -0,0 +1,466 @@
+# Various functions to read Athena++ output data files
+
+# Python modules
+import re
+import warnings
+from io import open  # Consistent binary I/O from Python 2 and 3
+import struct
+
+import matplotlib.colors as colors
+import matplotlib.patches as patches
+import matplotlib.pyplot as plt
+
+
+# Other Python modules
+import numpy as np
+
+check_nan_flag = False
+
+
+# Check input NumPy array for the presence of any NaN entries
+def check_nan(data):
+    if np.isnan(data).any():
+        raise FloatingPointError("NaN encountered")
+    return
+
+
+# Wrapper to np.loadtxt() for checks used in regression tests
+def error_dat(filename, **kwargs):
+    data = np.loadtxt(filename,
+                      dtype=np.float64,
+                      ndmin=2,  # prevent NumPy from squeezing singleton dim
+                      **kwargs)
+    if check_nan_flag:
+        check_nan(data)
+    return data
+
+
+# Read .tab files and return dict.
+def tab(filename, show_vars=False):
+
+    # Parse header
+    data_dict = {}
+    with open(filename, 'r') as data_file:
+        line = data_file.readline()
+        attributes = re.search(r'time=(\S+)\s+cycle=(\S+)', line)
+        line = data_file.readline()
+        headings = line.split()[1:]
+    headings = headings[1:]
+
+    # Go through lines
+    data_array = []
+    num_lines = 0
+    with open(filename, 'r') as data_file:
+        first_line = True
+        for line in data_file:
+            # Skip comments
+            if line.split()[0][0] == '#':
+                continue
+
+            # Extract cell indices
+            vals = line.split()
+            if first_line:
+                num_entries = len(vals) - 1
+                first_line = False
+
+            # Extract cell values
+            vals = vals[1:]
+            data_array.append([float(val) for val in vals])
+            num_lines += 1
+
+    # Reshape array
+    array_shape = (num_lines, num_entries)
+    array_transpose = (1, 0)
+    data_array = np.transpose(np.reshape(data_array, array_shape),
+                              array_transpose)
+
+
+    
+    # Finalize data
+    for n, heading in enumerate(headings):
+        if check_nan_flag:
+            check_nan(data_array[n, ...])
+        data_dict[heading] = data_array[n, ...]
+    if show_vars:
+        return list(data_dict.keys())
+    data_dict['time'] = float(attributes.group(1))
+    data_dict['cycle'] = int(attributes.group(2))
+    return data_dict
+
+
+# Read .hst files and return dict of 1D arrays.
+# Keyword arguments:
+# raw -- if True, do not prune file to remove stale data
+# from prev runs (default False)
+def hst(filename, raw=False):
+    # Read data
+    with open(filename, 'r') as data_file:
+        # Find header
+        header_found = False
+        multiple_headers = False
+        header_location = None
+        line = data_file.readline()
+        while len(line) > 0:
+            if line == '# Athena++ history data\n':
+                if header_found:
+                    multiple_headers = True
+                else:
+                    header_found = True
+                header_location = data_file.tell()
+            line = data_file.readline()
+        if multiple_headers:
+            warnings.warn('Multiple headers found; using most recent data')
+        if header_location is None:
+            raise RuntimeError('athena_read.hst: Could not find header')
+
+        # Parse header
+        data_file.seek(header_location)
+        header = data_file.readline()
+        data_names = re.findall(r'\[\d+\]=(\S+)', header)
+        if len(data_names) == 0:
+            raise RuntimeError('athena_read.hst: Could not parse header')
+
+        # Prepare dictionary of results
+        data = {}
+        for name in data_names:
+            data[name] = []
+
+        # Read data
+        for line in data_file:
+            for name, val in zip(data_names, line.split()):
+                data[name].append(float(val))
+
+    # Finalize data
+    for key, val in data.items():
+        data[key] = np.array(val)
+    if not raw:
+        if data_names[0] != 'time':
+            raise AthenaError('Cannot remove spurious data because time '
+                              'column could not be identified')
+        branches_removed = False
+        while not branches_removed:
+            branches_removed = True
+            for n in range(1, len(data['time'])):
+                if data['time'][n] <= data['time'][n-1]:
+                    branch_index = np.where((data['time'][:n] >=
+                                             data['time'][n]))[0][0]
+                    for key, val in data.items():
+                        data[key] = np.concatenate((val[:branch_index],
+                                                    val[n:]))
+                    branches_removed = False
+                    break
+        if check_nan_flag:
+            for key, val in data.items():
+                check_nan(val)
+    return data
+
+# Read .bin files and return dict with numpy array of variables and WCS
+# This is a Z-only code ripped from athenak's plot_slice.py
+# It returns not only all numpy arrays, but also a few meta-data
+# named:   'time', 'xlim', 'ylim'
+def bin(filename, show_vars=False, **kwargs):
+    
+    # Read data
+    with open(filename, 'rb') as f:
+
+        # Get file size
+        f.seek(0, 2)
+        file_size = f.tell()
+        f.seek(0, 0)
+
+        # Read header metadata
+        line = f.readline().decode('ascii')
+        if line != 'Athena binary output version=1.1\n':
+            print(line)
+            raise RuntimeError('Unrecognized data file format.')
+        next(f)
+        line = f.readline().decode('ascii')
+        if line[:7] != '  time=':
+            raise RuntimeError('Could not read time.')
+        sim_time = float(line[7:])
+        next(f)
+        line = f.readline().decode('ascii')
+        if line[:19] != '  size of location=':
+            raise RuntimeError('Could not read location size.')
+        location_size = int(line[19:])
+        line = f.readline().decode('ascii')
+        if line[:19] != '  size of variable=':
+            raise RuntimeError('Could not read variable size.')
+        variable_size = int(line[19:])
+        next(f)
+        line = f.readline().decode('ascii')
+        if line[:12] != '  variables:':
+            raise RuntimeError('Could not read variable names.')
+        variable_names_base = line[12:].split()
+        line = f.readline().decode('ascii')
+        if line[:16] != '  header offset=':
+            raise RuntimeError('Could not read header offset.')
+        header_offset = int(line[16:])
+
+        # Process header metadata
+        if location_size not in (4, 8):
+            raise RuntimeError('Only 4- and 8-byte integer types supported for '
+                               'location data.')
+        location_format = 'f' if location_size == 4 else 'd'
+        if variable_size not in (4, 8):
+            raise RuntimeError('Only 4- and 8-byte integer types supported for cell '
+                               'data.')
+        variable_format = 'f' if variable_size == 4 else 'd'
+        num_variables_base = len(variable_names_base)
+        if show_vars:
+            return variable_names_base
+            
+        if True:
+            variable_name = kwargs['variable']
+            if variable_name not in variable_names_base:
+                raise RuntimeError('Variable "{0}" not found; options are {{{1}}}.'
+                                   .format(variable_name,
+                                           ', '.join(variable_names_base)))
+            variable_names = [variable_name]
+            variable_ind = 0
+            while variable_names_base[variable_ind] != variable_name:
+                variable_ind += 1
+            variable_inds = [variable_ind]
+        variable_names_sorted = \
+            [name for _, name in sorted(zip(variable_inds, variable_names))]
+        variable_inds_sorted = \
+            [ind for ind, _ in sorted(zip(variable_inds, variable_names))]
+
+        # @todo loop over variables
+        retval = {}
+
+        # Read input file metadata
+        input_data = {}
+        start_of_data = f.tell() + header_offset
+        while f.tell() < start_of_data:
+            line = f.readline().decode('ascii')
+            if line[0] == '#':
+                continue
+            if line[0] == '<':
+                section_name = line[1:-2]
+                input_data[section_name] = {}
+                continue
+            key, val = line.split('=', 1)
+            input_data[section_name][key.strip()] = val.split('#', 1)[0].strip()
+
+        # Extract number of ghost cells from input file metadata
+        try:
+            num_ghost = int(input_data['mesh']['nghost'])
+        except:  # noqa: E722
+            raise RuntimeError('Unable to find number of ghost cells in input file.')
+
+        # Prepare lists to hold results
+        max_level_calculated = -1
+        block_loc_for_level = []
+        block_ind_for_level = []
+        num_blocks_used = 0
+        extents = []
+        quantities = {}
+        for name in variable_names_sorted:
+            quantities[name] = []
+
+        # Go through blocks
+        first_time = True
+        while f.tell() < file_size:
+
+            # Read grid structure data
+            block_indices = np.array(struct.unpack('@6i', f.read(24))) - num_ghost
+            block_i, block_j, block_k, block_level = struct.unpack('@4i', f.read(16))
+
+            # Process grid structure data
+            if first_time:
+                block_nx = block_indices[1] - block_indices[0] + 1
+                block_ny = block_indices[3] - block_indices[2] + 1
+                block_nz = block_indices[5] - block_indices[4] + 1
+                cells_per_block = block_nz * block_ny * block_nx
+                block_cell_format = '=' + str(cells_per_block) + variable_format
+                variable_data_size = cells_per_block * variable_size
+                if True:
+                    # if kwargs['dimension'] == 'z':
+                    if block_nx == 1:
+                        raise RuntimeError('Data in file has no extent in x-direction.')
+                    if block_ny == 1:
+                        raise RuntimeError('Data in file has no extent in y-direction.')
+                    block_nx1 = block_nx
+                    block_nx2 = block_ny
+                    slice_block_n = block_nz
+                    slice_location_min = float(input_data['mesh']['x3min'])
+                    slice_location_max = float(input_data['mesh']['x3max'])
+                    slice_root_blocks = (int(input_data['mesh']['nx3'])
+                                         // int(input_data['meshblock']['nx3']))
+                slice_normalized_coord = (kwargs['location'] - slice_location_min) \
+                    / (slice_location_max - slice_location_min)
+                first_time = False
+
+            # Determine if block is needed
+            if block_level > max_level_calculated:
+                for level in range(max_level_calculated + 1, block_level + 1):
+                    if kwargs['location'] <= slice_location_min:
+                        block_loc_for_level.append(0)
+                        block_ind_for_level.append(0)
+                    elif kwargs['location'] >= slice_location_max:
+                        block_loc_for_level.append(slice_root_blocks - 1)
+                        block_ind_for_level.append(slice_block_n - 1)
+                    else:
+                        slice_mesh_n = slice_block_n * slice_root_blocks * 2 ** level
+                        mesh_ind = int(slice_normalized_coord * slice_mesh_n)
+                        block_loc_for_level.append(mesh_ind // slice_block_n)
+                        block_ind_for_level.append(mesh_ind - slice_block_n
+                                                   * block_loc_for_level[-1])
+                max_level_calculated = block_level
+            # z
+            if block_k != block_loc_for_level[block_level]:
+                f.seek(6 * location_size + num_variables_base * variable_data_size, 1)
+                continue
+            num_blocks_used += 1
+
+            # Read coordinate data
+            block_lims = struct.unpack('=6' + location_format, f.read(6 * location_size))
+            # z
+            extents.append((block_lims[0], block_lims[1], block_lims[2],
+                            block_lims[3]))
+
+            # Read cell data
+            cell_data_start = f.tell()
+            for ind, name in zip(variable_inds_sorted, variable_names_sorted):
+                if ind == -1:
+                    # z
+                    quantities[name].append(np.full((block_ny, block_nx),
+                                                    block_level))
+                else:
+                    f.seek(cell_data_start + ind * variable_data_size, 0)
+                    cell_data = (np.array(struct.unpack(block_cell_format,
+                                                        f.read(variable_data_size)))
+                                 .reshape(block_nz, block_ny, block_nx))
+                    block_ind = block_ind_for_level[block_level]
+                    # z
+                    quantities[name].append(cell_data[block_ind, :, :])
+            f.seek((num_variables_base - ind - 1) * variable_data_size, 1)
+
+    # Prepare to calculate derived quantity
+    for name in variable_names_sorted:
+        quantities[name] = np.array(quantities[name])
+
+    # Extract quantity without derivation
+    quantity = quantities[variable_name]
+
+    if kwargs['output_file'] == None:
+        if num_blocks_used > 1:
+            raise RuntimeError('too many blocks, mesh and meshblock not the same')
+        quantities['time'] = sim_time
+
+        x1_min = float(input_data['mesh']['x1min'])
+        x1_max = float(input_data['mesh']['x1max'])
+        x2_min = float(input_data['mesh']['x2min'])
+        x2_max = float(input_data['mesh']['x2max'])
+        quantities['xlim'] = (x1_min,x1_max)
+        quantities['ylim'] = (x2_min,x2_max)
+            
+        return quantities
+
+    # Calculate colors
+    if kwargs['vmin'] is None:
+        vmin = np.nanmin(quantity)
+    else:
+        vmin = kwargs['vmin']
+    if kwargs['vmax'] is None:
+        vmax = np.nanmax(quantity)
+    else:
+        vmax = kwargs['vmax']
+
+    # Choose colormap norm
+    if kwargs['norm'] == 'linear':
+        norm = colors.Normalize(vmin, vmax)
+        vmin = None
+        vmax = None
+    elif kwargs['norm'] == 'log':
+        norm = colors.LogNorm(vmin, vmax)
+        vmin = None
+        vmax = None
+    else:
+        norm = kwargs['norm']
+
+    # Prepare figure
+    plt.figure()
+
+    x1_labelpad = 2.0
+    x2_labelpad = 2.0
+    dpi = 300
+    
+    # Plot data
+    for block_num in range(num_blocks_used):
+        d = quantity[block_num]
+        print("block_num:",block_num,d.shape,extents[block_num])
+        
+        plt.imshow(quantity[block_num], cmap=kwargs['cmap'], norm=norm, vmin=vmin,
+                   vmax=vmax, interpolation='none', origin='lower',
+                   extent=extents[block_num])
+    # Make colorbar
+    plt.colorbar()
+
+    # Adjust axes
+    # z
+    x1_min = float(input_data['mesh']['x1min'])
+    x1_max = float(input_data['mesh']['x1max'])
+    x2_min = float(input_data['mesh']['x2min'])
+    x2_max = float(input_data['mesh']['x2max'])
+    print("Mesh:  X: %g %g    Y: %g %g" % (x1_min,x1_max,x2_min,x2_max))
+    if kwargs['x1_min'] is not None:
+        x1_min = kwargs['x1_min']
+    if kwargs['x1_max'] is not None:
+        x1_max = kwargs['x1_max']
+    if kwargs['x2_min'] is not None:
+        x2_min = kwargs['x2_min']
+    if kwargs['x2_max'] is not None:
+        x2_max = kwargs['x2_max']
+    plt.xlim((x1_min, x1_max))
+    plt.ylim((x2_min, x2_max))
+    # z
+    plt.xlabel('$x$', labelpad=x1_labelpad)
+    plt.ylabel('$y$', labelpad=x2_labelpad)
+    
+    # Adjust layout
+    plt.tight_layout()
+
+
+    # Save or display figure
+    if kwargs['output_file'] != 'show':
+        plt.savefig(kwargs['output_file'], dpi=dpi)
+    else:
+        plt.show()
+
+
+# General exception class for these functions
+class AthenaError(RuntimeError):
+    pass
+
+
+# testing the bin function
+if __name__ == "__main__":
+    import sys
+    
+    kwargs = {}
+    kwargs['variable'] = 'dens'
+    kwargs['dimension'] = 'z'
+    kwargs['location'] = 0
+    kwargs['vmin'] = None
+    kwargs['vmax'] = None
+    kwargs['norm'] = 'linear'
+    kwargs['cmap'] = 'viridis'
+    kwargs['x1_min'] = None
+    kwargs['x1_max'] = None
+    kwargs['x2_min'] = None
+    kwargs['x2_max'] = None
+    kwargs['x2_max'] = None
+    kwargs['output_file'] = 'show'
+    # kwargs['output_file'] = None
+
+    if False:
+        print(bin(sys.argv[1],True))
+        d = bin(sys.argv[1],False,**kwargs)    
+        print('data',d)
+    else:
+        print(tab(sys.argv[1],True))
+        d = tab(sys.argv[1],False)
+        print('data',d)
+        
diff --git a/plot1d.py b/plot1d.py
index b9c503a..bcad50a 100755
--- a/plot1d.py
+++ b/plot1d.py
@@ -161,7 +161,7 @@ def mouse_moved(e):
 argparser.add_argument('-f', '--fix', action='store_true', help='fixes the x and y axes of the animation based on the animation\'s first frame')
 args = argparser.parse_args()
 
-# fnames='run1/tab/LinWave*tab'
+# @todo   there can be several ID types of the form   BASENAME.ID.NNNNN.tab
 if args.hst:
     f = glob.glob(args.dir + '/*.hst')
 else:
diff --git a/plot2d.py b/plot2d.py
new file mode 100755
index 0000000..a8cead8
--- /dev/null
+++ b/plot2d.py
@@ -0,0 +1,319 @@
+#! /usr/bin/env python
+#
+#     plot2d:    animate images
+#
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.widgets import RadioButtons, Button, Slider, CheckButtons, TextBox
+from argparse import ArgumentParser
+import athena_read
+import glob
+import os
+import sys
+import matplotlib.style as mplstyle
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+#mplstyle.use(['ggplot', 'fast'])
+
+# TODO list
+# round buttons
+# round check boxes
+# nonlinear (logarithmic?) slider
+
+kwargs = {}
+kwargs['variable'] = 'dens'
+kwargs['dimension'] = 'z'
+kwargs['location'] = 0
+kwargs['vmin'] = None
+kwargs['vmax'] = None
+kwargs['norm'] = 'linear'
+kwargs['cmap'] = 'rainbow'
+kwargs['cmap'] = 'gist_rainbow'
+kwargs['cmap'] = 'viridis'
+kwargs['x1_min'] = None
+kwargs['x1_max'] = None
+kwargs['x2_min'] = None
+kwargs['x2_max'] = None
+kwargs['output_file'] = None
+
+zvar = 'dens'
+
+# https://saturncloud.io/blog/how-to-animate-the-colorbar-in-matplotlib-a-guide/
+
+def animate(i):
+    global xcol, ycol, current_frame, xlim, ylim, kwargs
+    kwargs['variable'] = zvar
+    if False:
+        d = athena_read.bin(f[i],False,**kwargs)
+    else:
+        d = data[i][zvar][0]
+    time = data[i]['time']
+    xlim = data[i]['xlim']
+    ylim = data[i]['ylim']
+    extent=[xlim[0],xlim[1],ylim[0],ylim[1]]
+    #sax.set_xlim((-0.5,0.5))
+    
+
+    ax.clear()
+    #ax.set_xlim(xlim)
+    #ax.set_ylim(ylim)
+    if True:
+        im = ax.imshow(d,cmap=kwargs['cmap'],  # norm=norm, vmin=vmin, vmax=vmax,
+                       interpolation='none', origin='lower', extent=extent)
+                       
+        
+        # @todo    this fig.colorbar() will recursively die
+        # fig.colorbar(im, cax=cax, orientation='vertical')
+    else:
+        # now broken after using extent=
+        im = ax.imshow(d,cmap=kwargs['cmap'],  # norm=norm, vmin=vmin, vmax=vmax,
+                       interpolation='none', origin='lower', extent=extent)
+        x = np.arange(d.shape[0])
+        y = np.arange(d.shape[1])
+        cs = ax.contour(x,y,d,extent=extent)
+                       
+    if not xlim and args.fix: # just need to check one since its either both or none
+        xlim = ax.get_xlim()
+        ylim = ax.get_ylim()
+    ax.set_title(f'{zvar}  Time: {float(time)}', loc='left')
+    ax.set_xlabel('x')
+    ax.set_ylabel('y')
+
+# hard-pauses the animation
+# that is, the only way to unpause is by using the play or restart button
+def hpause(self=None):
+    global hard_paused, is_playing
+    hard_paused = True
+    is_playing = False
+    #bpause.color = '0.5'
+    bplay.label.set_text('$\u25B6$')
+    pause()
+
+def play(self=None):
+    global is_playing
+    if is_playing:
+        hpause()
+    else:
+        resume()
+
+# pauses the animation
+def pause(self=None):
+    global is_playing, fig
+    if is_playing:
+        fig.canvas.stop_event_loop()
+        is_playing = False
+
+# plays the animation (either starts or resumes)
+def resume(self=None):
+    global is_playing, current_frame, ax, length, frame_sliding, hard_paused
+    if not is_playing:
+        bplay.label.set_text('$\u25A0$')
+        is_playing = True
+        hard_paused = False
+        #bpause.color = '0.85'
+        while current_frame < length and is_playing:
+            animate(current_frame)
+            current_frame += 1
+            fig.canvas.draw_idle()
+            fig.canvas.start_event_loop(delay)
+            frame_sliding = False
+            fslider.set_val(current_frame)
+        is_playing = False
+        if loop and current_frame == length:
+            restart()
+
+# loops the animation
+def loopf(self=None):
+    global loop, bloop
+    if loop:
+        bloop.color = '0.85' # this is the default color
+    else:
+        bloop.color = '0.5' # highlighted color
+    loop = not loop
+
+# restarts the animation from the beginning
+def restart(self=None):
+    global current_frame, xlim
+    xlim = None
+    pause()
+    current_frame=0
+    resume()
+
+# select the variable
+def select_v(label):
+    global current_frame, xcol, xlim, zvar
+    update_cols(xcol, label)
+    if not hard_paused:
+        restart()
+    else:
+        xlim = None
+        animate(current_frame)
+        fig.canvas.draw_idle()
+    print('select_v',label)
+    zvar = label
+    kwargs['variable'] = zvar
+    reload_data()
+
+def update_cols(x, y):
+    global xcol, ycol, ixcol, iycol
+    xcol=x
+    ycol=y
+    ixcol = variables.index(x)
+    iycol = variables.index(y)
+
+def update_delay(x):
+    global delay
+    delay = x / 1000
+    #  something odd about setting the delay
+    # print("DELAY:",delay)
+
+def update_fslider(n):
+    global frame_sliding, current_frame
+    current_frame = n
+    animate(current_frame - 1)
+    fig.canvas.draw_idle()
+
+'''def fix_axes(v):
+    global xlim, ylim
+    if v == 'Fix X':
+        xlim = ax.get_xlim() if not xlim else None
+    else:
+        ylim = ax.get_ylim() if not ylim else None'''
+
+def mouse_moved(e):
+    if e.inaxes == fax:
+        pause()
+    elif not hard_paused:
+        resume()
+
+def reload_data():
+    global data
+    print("Reading %d data-frames for %s" % (length,zvar))
+    for i in range(len(f)):
+        data[i] = athena_read.bin(f[i],False,**kwargs)
+    
+
+argparser = ArgumentParser(description='plots the athena tab files specified')
+argparser.add_argument('-d', '--dir', help='the athena run directory containing tab/ or *.tab files', required=True)
+argparser.add_argument('-n', '--name', help='name of the problem being plotted') # primarily just used by the other gui
+argparser.add_argument('-f', '--fix', action='store_true', help='fixes the x and y axes of the animation based on the animation\'s first frame')
+args = argparser.parse_args()
+
+f = glob.glob(args.dir + '/bin/*.bin')
+f.sort()
+length = len(f)
+
+if length==0:
+    print("No bin data found")
+    sys.exit(0)
+
+data = list(range(len(f)))
+reload_data()
+
+
+#print('DEBUG: %s has %d files' % (fnames,len(f)))
+
+# global vars
+current_frame = 0
+is_playing = False
+hard_paused = True
+loop = False
+frame_sliding = False
+xlim = None
+ylim = None
+
+# plot settings
+left = 0.34
+bottom = 0.25
+top = 0.85
+
+# the time in seconds between frames
+delay = 100 / 1000
+
+# getting the variable names
+variables = athena_read.bin(f[0],True)
+print("bin variables detected:",variables)
+var_len = len(variables)
+
+# 0-based, change from animate2
+xcol=variables[0]
+ycol=variables[0]
+ixcol = 0
+iycol = 0
+
+# plotting configuration
+fig, ax = plt.subplots()
+fig.subplots_adjust(left=left, bottom=bottom, top=top) # old bottom was 0.34
+divider = make_axes_locatable(ax)
+cax = divider.append_axes('right', size='5%', pad=0.05)
+sax = ax.secondary_xaxis('bottom')
+# plt.rcParams['font.family'] = 'Arial'
+# pause on close otherwise we might freeze
+fig.canvas.mpl_connect('close_event', pause)
+# fig.set_size_inches(10, 10)
+
+plt.get_current_fig_manager().set_window_title(args.name if args.name else f[0].split('.')[0])
+
+rheight = var_len / 25
+rwidth = 0.045
+rdleft = 0.03
+rbot = (bottom + top - rheight) / 2
+
+rax = fig.add_axes([rdleft, rbot, rwidth, rheight])
+
+
+# use same axes to add text in order to make it easier to adjust
+# rax.text(-0.055, 0.05, 'X')
+rax.text(0.055, 0.05, 'Y')
+
+# rax = fig.add_axes([rdleft + 0.015, rbot, 0.25, rheight])
+radio2 = RadioButtons(rax, 
+                      tuple(variables), 
+                      radio_props={'color': ['#1f77b4' for _ in variables], 'edgecolor': ['black' for _ in variables]}
+                      )
+rax.axis('off')
+radio2.on_clicked(select_v)
+
+if True:
+
+    bwidth = 0.04
+    bheight = 0.05
+    bspace = 0.02
+    bstart = 0.05
+
+    brestart = Button(fig.add_axes([bstart, 0.125, bwidth, bheight]), '$\u21A9$')
+    brestart.on_clicked(restart)
+
+    bplay = Button(fig.add_axes([bstart + bwidth + bspace, 0.125, bwidth, bheight]), '$\u25B6$')
+    bplay.on_clicked(play)
+
+    bloop = Button(fig.add_axes([bstart + 2 * bwidth + 2 * bspace, 0.125, bwidth, bheight]), '$\u27F3$')
+    bloop.on_clicked(loopf)
+
+    # make slider nonlinear
+    delay_slider = Slider(
+        ax=fig.add_axes([0.18, 0.05, 0.65, 0.03]),
+        label='Delay (ms)',
+        valmin=1,
+        valmax=20,
+        valinit=5,
+    )
+
+    delay_slider.on_changed(update_delay)
+
+    fax = fig.add_axes([0.18, 0.95, 0.65, 0.03])
+    fslider = Slider(
+        ax=fax,
+        label='Frame',
+        valmin=1,
+        valmax=length,
+        valinit=1,
+        valstep=1
+    )
+
+    fslider.on_changed(update_fslider)
+
+
+    # in order to pause the animation when using the frame slider
+    fig.canvas.mpl_connect('motion_notify_event', mouse_moved)
+
+plt.show()
diff --git a/pythena_run.py b/pythena_run.py
index 725bc17..892b938 100755
--- a/pythena_run.py
+++ b/pythena_run.py
@@ -237,10 +237,13 @@ def run(self, odir_input, plot):
             #Popen(['python', 'plot1d.py', '-d', odir,         '--hst', '-n', info['problem'] + ' history'])
             cmd1 = 'python plot1d.py  -d %s/tab' % (odir)
             cmd2 = 'python plot1d.py  -d %s --hst' % (odir)
+            cmd3 = 'python plot2d.py  -d %s' % (odir)            
             print(cmd1)
             print(cmd2)
+            print(cmd3)
             Popen(cmd1.split())
             Popen(cmd2.split())
+            Popen(cmd3.split())
         else:
             w = DisplayWindow(cmd)
             self.windows.append(w)
diff --git a/qtrun.py b/qtrun.py
old mode 100644
new mode 100755
index 973ac63..d5c6312
--- a/qtrun.py
+++ b/qtrun.py
@@ -1,3 +1,5 @@
+#! /usr/bin/env python
+#
 import sys
 import os
 from PyQt5 import QtCore, QtWidgets
@@ -60,6 +62,7 @@ def initUI(self):
         self.createWidgetsFromGroups()
     
     # runs the file. takes the input file and runs it, piping the set paraments into the file as args
+    # it assumes the script can self-execute with the appropriate #! directive on the first line
     def run(self):
         print("running: " + self.input_file)
         contents = self.gather_data()
@@ -354,4 +357,4 @@ def parsefile(file):
         print('opening window')
         sys.exit(app.exec())
     except SystemExit:
-        print('closing window')
\ No newline at end of file
+        print('closing window')
diff --git a/version b/version
new file mode 100755
index 0000000..2447ac4
--- /dev/null
+++ b/version
@@ -0,0 +1,2 @@
+cd $(dirname $0)
+git rev-list --count HEAD