From 81576eaf412d49449c7679ad97c542c157dd256d Mon Sep 17 00:00:00 2001
From: Peter Teuben <teuben@gmail.com>
Date: Fri, 25 Aug 2023 09:57:27 -0400
Subject: [PATCH] return more meta-data with the numpy array

---
 athena_read.py | 32 +++++++++++++++++++++++---------
 1 file changed, 23 insertions(+), 9 deletions(-)

diff --git a/athena_read.py b/athena_read.py
index b2f2542..8afb1da 100644
--- a/athena_read.py
+++ b/athena_read.py
@@ -154,7 +154,10 @@ def hst(filename, raw=False):
                 check_nan(val)
     return data
 
-# Read .bin files and return dict ?
+# 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
@@ -171,7 +174,10 @@ def bin(filename, show_vars=False, **kwargs):
             print(line)
             raise RuntimeError('Unrecognized data file format.')
         next(f)
-        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=':
@@ -203,7 +209,6 @@ def bin(filename, show_vars=False, **kwargs):
         num_variables_base = len(variable_names_base)
         if show_vars:
             return variable_names_base
-        # print("variable_names_base: ", variable_names_base)
             
         if True:
             variable_name = kwargs['variable']
@@ -221,6 +226,9 @@ def bin(filename, show_vars=False, **kwargs):
         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
@@ -336,13 +344,19 @@ def bin(filename, show_vars=False, **kwargs):
     # Extract quantity without derivation
     quantity = quantities[variable_name]
 
-    # print("PJT",num_blocks_used,quantity.shape)
-
     if kwargs['output_file'] == None:
-        if num_blocks_used == 1:
-            return quantity[0]
-        print("Cannot merge blocks yet")
-        return quantity[0]
+        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: