Merge 2

osh5io.py

@@ -167,13 +167,66 @@ def read_raw(filename, path=None):
         #TODO: TIMESTAMP is not set in HDF5 file as of now (Aug 2019) so we make one up, check back when file format changes
         d['TIMESTAMP'] = timestamp
 
-        dtype = [(q, data[q].dtype) for q in quants]
+
+        dtype = [(q, data[q].dtype, (2,)) if q.lower() == 'tag' else (q, data[q].dtype) for q in quants]
+        # dtype = [(q, data[q].dtype) for q in quants]
+        print('dtype=',dtype)
+        print(data[dtype[0][0]].shape)
         r = PartData(data[dtype[0][0]].shape, dtype=dtype, attrs=d)
         for dt in dtype:
             r[dt[0]] = data[dt[0]]
 
     return r
 
+def read_h5_tracks(filename, path=None):
+    """
+    subroutine that reads a OSIRIS track file into an ordered list of NUMPY arrays.
+    Currently this subroutines supports the older (ver 1) track files
+
+    and the attributes that describe the data (e.g. title, units, scale).
+
+    tracks, quants = read_h5_tracks(filename)
+    print(tracks.shape())           # shows how many tracks are in the file
+    print(tracks[0].shape)          # this prints 2 numbers, the first number is the total number 
+                                    # of timesteps for the first track, and the second number is the number
+                                    # of particle quantities are stored in the track 
+    print(quants.index['x1'])       # this prints the array index for the quantity 'x1',
+                                    # possible track quantities include 't', 'x1','x2','x3','p1','p2','p3'
+                                    # 'E1', 'E2', 'E3', 'B1', 'B2', 'B3', and 'ene', for quasi 3D there is 
+                                    # an 'x4'
+    plt.plot(tracks[i][:,quants.index('p1')],tracks[i][:,quants.index('p2')])
+                                    # plots p1 vs p2 for the i-th particle
+
+    ## Example Python Code 
+
+    (a,b) = read_h5_tracks('FILENAME')
+    plt.figure()
+    for j in range(len(a)):
+        plt.plot(a[j][:,b.index('p1')],a[j][:b.index('p2')])
+    plt.show()
+
+    """
+    fname = filename if not path else path + '/' + filename
+    data_file = h5py.File(fname, 'r')
+    keys = list(data_file.keys())
+    num_tracks = len(keys)
+    quants = data_file.attrs['QUANTS']
+    quants_string = []
+    # here we convert the quant array from bytes to string 
+    for x in quants:
+        quants_string.append(x.decode())
+    # 
+
+    quants_string = []
+    particle_list = []
+
+    for j in range(0,num_tracks):
+        key = keys[j]
+        track = data_file[key]
+        particle_list.append(track)
+
+    return(particle_list, quants_string)
+
 
 def read_h5_openpmd(filename, path=None):
     """
@@ -523,7 +576,8 @@ def write_h5_openpmd(data, filename=None, path=None, dataset_name=None, overwrit
         local_offset[0]= np.float32(0.0)
 
     elif (number_axis_objects_we_need == 2):
-        local_axislabels=[b'x1', b'x2']
+        # local_axislabels=[b'x1', b'x2']
+        local_axislabels=np.array([np.string_('x1'),np.string('x2')])
         deltax[0] = deltax[0]/data_object.shape[0]
         deltax[1] = deltax[1]/data_object.shape[1]
         temp=deltax[0]
@@ -538,7 +592,9 @@ def write_h5_openpmd(data, filename=None, path=None, dataset_name=None, overwrit
         local_offset[1]= np.float32(0.0)
 
     else:
-        local_axislabels=[b'x1',b'x2',b'x3']
+        # local_axislabels=[b'x1',b'x2',b'x3']
+        # local_axislabels=['x1','x2','x3']
+        local_axislabels=np.array([np.string_('x1'),np.string_('x2'),np.string_('x3')])
         deltax[0] = deltax[0]/data_object.shape[0]
         deltax[1] = deltax[1]/data_object.shape[1]
         deltax[2] = deltax[2]/data_object.shape[2]
@@ -557,9 +613,10 @@ def write_h5_openpmd(data, filename=None, path=None, dataset_name=None, overwrit
 
 
 
-    datasetid.attrs['dataOrder'] = 'F'
-    datasetid.attrs['geometry'] = 'cartesian'
-    datasetid.attrs['geometryParameters'] =  'cartesian'
+    datasetid.attrs['dataOrder'] = np.string_('F')
+    datasetid.attrs['geometry'] = np.string_('cartesian')
+    datasetid.attrs['geometryParameters'] =  np.string_('cartesian')
+    print("axislabels=",local_axislabels)
     datasetid.attrs['axisLabels'] = local_axislabels
     datasetid.attrs['gridUnitSI'] = np.float64(length_to_si)
     datasetid.attrs['unitSI'] = np.float64(data_to_si)

osh5utils.py

@@ -593,7 +593,7 @@ def field_decompose(fldarr, ffted=True, idim=None, finalize=None, outquants=('L'
     if not finalize:
         finalize = __idle
     if np.issubdtype(fldarr[0].dtype, np.floating):
-        ftf, iftf = fftn, ifftn if norft else rfftn, irfftn
+        ftf, iftf = (fftn, partial(ifftn, result=np.real(result))) if norft else (rfftn, irfftn)
 
     def wrap_up(data):
         if idim:

osh5vis.py

@@ -68,7 +68,7 @@ def default_title(h5data, show_time=True, title=None, **kwargs):
 
 
 def tex(s):
-    return '$' + s + '$' if s else ''
+    return '$' + str(s) + '$' if s else ''
 
 
 def axis_format(name=None, unit=None):